MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
solvers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "linalg.hpp"
13 #include "../general/annotation.hpp"
14 #include "../general/forall.hpp"
15 #include "../general/globals.hpp"
16 #include "../fem/bilinearform.hpp"
17 #include <iostream>
18 #include <iomanip>
19 #include <algorithm>
20 #include <cmath>
21 #include <set>
22 
23 namespace mfem
24 {
25 
26 using namespace std;
27 
29  : Solver(0, true)
30 {
31  oper = NULL;
32  prec = NULL;
33  max_iter = 10;
34  rel_tol = abs_tol = 0.0;
35 #ifdef MFEM_USE_MPI
36  dot_prod_type = 0;
37 #endif
38 }
39 
40 #ifdef MFEM_USE_MPI
41 
43  : Solver(0, true)
44 {
45  oper = NULL;
46  prec = NULL;
47  max_iter = 10;
48  rel_tol = abs_tol = 0.0;
49  dot_prod_type = 1;
50  comm = comm_;
51 }
52 
53 #endif // MFEM_USE_MPI
54 
55 double IterativeSolver::Dot(const Vector &x, const Vector &y) const
56 {
57 #ifndef MFEM_USE_MPI
58  return (x * y);
59 #else
60  if (dot_prod_type == 0)
61  {
62  return (x * y);
63  }
64  else
65  {
66  return InnerProduct(comm, x, y);
67  }
68 #endif
69 }
70 
71 void IterativeSolver::SetPrintLevel(int print_lvl)
72 {
74  int print_level_ = print_lvl;
75 
76 #ifdef MFEM_USE_MPI
77  if (dot_prod_type != 0)
78  {
79  int rank;
80  MPI_Comm_rank(comm, &rank);
81  if (rank != 0) // Suppress output.
82  {
83  print_level_ = -1;
85  }
86  }
87 #endif
88 
89  print_level = print_level_;
90 }
91 
93 {
94  print_options = options;
95 
96  int derived_print_level = GuessLegacyPrintLevel(options);
97 
98 #ifdef MFEM_USE_MPI
99  if (dot_prod_type != 0)
100  {
101  int rank;
102  MPI_Comm_rank(comm, &rank);
103  if (rank != 0)
104  {
105  derived_print_level = -1;
107  }
108  }
109 #endif
110 
111  print_level = derived_print_level;
112 }
113 
115  int print_level_)
116 {
117 #ifdef MFEM_USE_MPI
118  int rank = 0;
119  if (comm != MPI_COMM_NULL)
120  {
121  MPI_Comm_rank(comm, &rank);
122  }
123 #endif
124 
125  switch (print_level_)
126  {
127  case -1:
128  return PrintLevel();
129  case 0:
130  return PrintLevel().Errors().Warnings();
131  case 1:
132  return PrintLevel().Errors().Warnings().Iterations();
133  case 2:
134  return PrintLevel().Errors().Warnings().Summary();
135  case 3:
136  return PrintLevel().Errors().Warnings().FirstAndLast();
137  default:
138 #ifdef MFEM_USE_MPI
139  if (rank == 0)
140 #endif
141  {
142  MFEM_WARNING("Unknown print level " << print_level_ <<
143  ". Defaulting to level 0.");
144  }
145  return PrintLevel().Errors().Warnings();
146  }
147 }
148 
150 {
151  if (print_options_.iterations)
152  {
153  return 1;
154  }
155  else if (print_options_.first_and_last)
156  {
157  return 3;
158  }
159  else if (print_options_.summary)
160  {
161  return 2;
162  }
163  else if (print_options_.errors && print_options_.warnings)
164  {
165  return 0;
166  }
167  else
168  {
169  return -1;
170  }
171 }
172 
174 {
175  prec = &pr;
176  prec->iterative_mode = false;
177 }
178 
180 {
181  oper = &op;
182  height = op.Height();
183  width = op.Width();
184  if (prec)
185  {
186  prec->SetOperator(*oper);
187  }
188 }
189 
190 void IterativeSolver::Monitor(int it, double norm, const Vector& r,
191  const Vector& x, bool final) const
192 {
193  if (monitor != nullptr)
194  {
195  monitor->MonitorResidual(it, norm, r, final);
196  monitor->MonitorSolution(it, norm, x, final);
197  }
198 }
199 
201  : damping(dmpng),
202  ess_tdof_list(nullptr),
203  oper(nullptr),
204  allow_updates(true)
205 { }
206 
208  const Array<int> &ess_tdofs,
209  const double dmpng)
210  :
211  Solver(a.FESpace()->GetTrueVSize()),
212  dinv(height),
213  damping(dmpng),
214  ess_tdof_list(&ess_tdofs),
215  residual(height),
216  allow_updates(false)
217 {
218  Vector &diag(residual);
219  a.AssembleDiagonal(diag);
220  // 'a' cannot be used for iterative_mode == true because its size may be
221  // different.
222  oper = nullptr;
223  Setup(diag);
224 }
225 
227  const Array<int> &ess_tdofs,
228  const double dmpng)
229  :
230  Solver(d.Size()),
231  dinv(height),
232  damping(dmpng),
233  ess_tdof_list(&ess_tdofs),
234  residual(height),
235  oper(NULL),
236  allow_updates(false)
237 {
238  Setup(d);
239 }
240 
242 {
243  if (!allow_updates)
244  {
245  // original behavior of this method
246  oper = &op; return;
247  }
248 
249  // Treat (Par)BilinearForm objects as a special case since their
250  // AssembleDiagonal method returns the true-dof diagonal whereas the form
251  // itself may act as an ldof operator. This is for compatibility with the
252  // constructor that takes a BilinearForm parameter.
253  const BilinearForm *blf = dynamic_cast<const BilinearForm *>(&op);
254  if (blf)
255  {
256  // 'a' cannot be used for iterative_mode == true because its size may be
257  // different.
258  oper = nullptr;
259  height = width = blf->FESpace()->GetTrueVSize();
260  }
261  else
262  {
263  oper = &op;
264  height = op.Height();
265  width = op.Width();
266  MFEM_ASSERT(height == width, "not a square matrix!");
267  // ess_tdof_list is only used with BilinearForm
268  ess_tdof_list = nullptr;
269  }
270  dinv.SetSize(height);
271  residual.SetSize(height);
272  Vector &diag(residual);
273  op.AssembleDiagonal(diag);
274  Setup(diag);
275 }
276 
278 {
279  residual.UseDevice(true);
280  const double delta = damping;
281  auto D = diag.Read();
282  auto DI = dinv.Write();
283  const bool use_abs_diag_ = use_abs_diag;
284  MFEM_FORALL(i, height,
285  {
286  if (D[i] == 0.0)
287  {
288  MFEM_ABORT_KERNEL("Zero diagonal entry in OperatorJacobiSmoother");
289  }
290  if (!use_abs_diag_) { DI[i] = delta / D[i]; }
291  else { DI[i] = delta / std::abs(D[i]); }
292  });
293  if (ess_tdof_list && ess_tdof_list->Size() > 0)
294  {
295  auto I = ess_tdof_list->Read();
296  MFEM_FORALL(i, ess_tdof_list->Size(), DI[I[i]] = delta; );
297  }
298 }
299 
301 {
302  // For empty MPI ranks, height may be 0:
303  // MFEM_VERIFY(Height() > 0, "The diagonal hasn't been computed.");
304  MFEM_ASSERT(x.Size() == Width(), "invalid input vector");
305  MFEM_ASSERT(y.Size() == Height(), "invalid output vector");
306 
307  if (iterative_mode)
308  {
309  MFEM_VERIFY(oper, "iterative_mode == true requires the forward operator");
310  oper->Mult(y, residual); // r = A y
311  subtract(x, residual, residual); // r = x - A y
312  }
313  else
314  {
315  residual = x;
316  y.UseDevice(true);
317  y = 0.0;
318  }
319  auto DI = dinv.Read();
320  auto R = residual.Read();
321  auto Y = y.ReadWrite();
322  MFEM_FORALL(i, height,
323  {
324  Y[i] += DI[i] * R[i];
325  });
326 }
327 
329  const Vector &d,
330  const Array<int>& ess_tdofs,
331  int order_, double max_eig_estimate_)
332  :
333  Solver(d.Size()),
334  order(order_),
335  max_eig_estimate(max_eig_estimate_),
336  N(d.Size()),
337  dinv(N),
338  diag(d),
339  coeffs(order),
340  ess_tdof_list(ess_tdofs),
341  residual(N),
342  oper(&oper_) { Setup(); }
343 
344 #ifdef MFEM_USE_MPI
346  const Vector &d,
347  const Array<int>& ess_tdofs,
348  int order_, MPI_Comm comm, int power_iterations, double power_tolerance)
349 #else
351  const Vector &d,
352  const Array<int>& ess_tdofs,
353  int order_, int power_iterations, double power_tolerance)
354 #endif
355  : Solver(d.Size()),
356  order(order_),
357  N(d.Size()),
358  dinv(N),
359  diag(d),
360  coeffs(order),
361  ess_tdof_list(ess_tdofs),
362  residual(N),
363  oper(&oper_)
364 {
365  OperatorJacobiSmoother invDiagOperator(diag, ess_tdofs, 1.0);
366  ProductOperator diagPrecond(&invDiagOperator, oper, false, false);
367 
368 #ifdef MFEM_USE_MPI
369  PowerMethod powerMethod(comm);
370 #else
371  PowerMethod powerMethod;
372 #endif
373  Vector ev(oper->Width());
374  max_eig_estimate = powerMethod.EstimateLargestEigenvalue(diagPrecond, ev,
375  power_iterations, power_tolerance);
376 
377  Setup();
378 }
379 
381  const Vector &d,
382  const Array<int>& ess_tdofs,
383  int order_, double max_eig_estimate_)
384  : OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, max_eig_estimate_) { }
385 
386 #ifdef MFEM_USE_MPI
388  const Vector &d,
389  const Array<int>& ess_tdofs,
390  int order_, MPI_Comm comm, int power_iterations, double power_tolerance)
391  : OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, comm,
392  power_iterations, power_tolerance) { }
393 #else
395  const Vector &d,
396  const Array<int>& ess_tdofs,
397  int order_, int power_iterations, double power_tolerance)
398  : OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, power_iterations,
399  power_tolerance) { }
400 #endif
401 
403 {
404  // Invert diagonal
405  residual.UseDevice(true);
406  auto D = diag.Read();
407  auto X = dinv.Write();
408  MFEM_FORALL(i, N, X[i] = 1.0 / D[i]; );
409  auto I = ess_tdof_list.Read();
410  MFEM_FORALL(i, ess_tdof_list.Size(), X[I[i]] = 1.0; );
411 
412  // Set up Chebyshev coefficients
413  // For reference, see e.g., Parallel multigrid smoothing: polynomial versus
414  // Gauss-Seidel by Adams et al.
415  double upper_bound = 1.2 * max_eig_estimate;
416  double lower_bound = 0.3 * max_eig_estimate;
417  double theta = 0.5 * (upper_bound + lower_bound);
418  double delta = 0.5 * (upper_bound - lower_bound);
419 
420  switch (order-1)
421  {
422  case 0:
423  {
424  coeffs[0] = 1.0 / theta;
425  break;
426  }
427  case 1:
428  {
429  double tmp_0 = 1.0/(pow(delta, 2) - 2*pow(theta, 2));
430  coeffs[0] = -4*theta*tmp_0;
431  coeffs[1] = 2*tmp_0;
432  break;
433  }
434  case 2:
435  {
436  double tmp_0 = 3*pow(delta, 2);
437  double tmp_1 = pow(theta, 2);
438  double tmp_2 = 1.0/(-4*pow(theta, 3) + theta*tmp_0);
439  coeffs[0] = tmp_2*(tmp_0 - 12*tmp_1);
440  coeffs[1] = 12/(tmp_0 - 4*tmp_1);
441  coeffs[2] = -4*tmp_2;
442  break;
443  }
444  case 3:
445  {
446  double tmp_0 = pow(delta, 2);
447  double tmp_1 = pow(theta, 2);
448  double tmp_2 = 8*tmp_0;
449  double tmp_3 = 1.0/(pow(delta, 4) + 8*pow(theta, 4) - tmp_1*tmp_2);
450  coeffs[0] = tmp_3*(32*pow(theta, 3) - 16*theta*tmp_0);
451  coeffs[1] = tmp_3*(-48*tmp_1 + tmp_2);
452  coeffs[2] = 32*theta*tmp_3;
453  coeffs[3] = -8*tmp_3;
454  break;
455  }
456  case 4:
457  {
458  double tmp_0 = 5*pow(delta, 4);
459  double tmp_1 = pow(theta, 4);
460  double tmp_2 = pow(theta, 2);
461  double tmp_3 = pow(delta, 2);
462  double tmp_4 = 60*tmp_3;
463  double tmp_5 = 20*tmp_3;
464  double tmp_6 = 1.0/(16*pow(theta, 5) - pow(theta, 3)*tmp_5 + theta*tmp_0);
465  double tmp_7 = 160*tmp_2;
466  double tmp_8 = 1.0/(tmp_0 + 16*tmp_1 - tmp_2*tmp_5);
467  coeffs[0] = tmp_6*(tmp_0 + 80*tmp_1 - tmp_2*tmp_4);
468  coeffs[1] = tmp_8*(tmp_4 - tmp_7);
469  coeffs[2] = tmp_6*(-tmp_5 + tmp_7);
470  coeffs[3] = -80*tmp_8;
471  coeffs[4] = 16*tmp_6;
472  break;
473  }
474  default:
475  MFEM_ABORT("Chebyshev smoother not implemented for order = " << order);
476  }
477 }
478 
480 {
481  if (iterative_mode)
482  {
483  MFEM_ABORT("Chebyshev smoother not implemented for iterative mode");
484  }
485 
486  if (!oper)
487  {
488  MFEM_ABORT("Chebyshev smoother requires operator");
489  }
490 
491  residual = x;
492  helperVector.SetSize(x.Size());
493 
494  y.UseDevice(true);
495  y = 0.0;
496 
497  for (int k = 0; k < order; ++k)
498  {
499  // Apply
500  if (k > 0)
501  {
502  oper->Mult(residual, helperVector);
503  residual = helperVector;
504  }
505 
506  // Scale residual by inverse diagonal
507  const int n = N;
508  auto Dinv = dinv.Read();
509  auto R = residual.ReadWrite();
510  MFEM_FORALL(i, n, R[i] *= Dinv[i]; );
511 
512  // Add weighted contribution to y
513  auto Y = y.ReadWrite();
514  auto C = coeffs.Read();
515  MFEM_FORALL(i, n, Y[i] += C[k] * R[i]; );
516  }
517 }
518 
520 {
521  r.SetSize(width);
522  z.SetSize(width);
523 }
524 
525 void SLISolver::Mult(const Vector &b, Vector &x) const
526 {
527  int i;
528 
529  // Optimized preconditioned SLI with fixed number of iterations and given
530  // initial guess
531  if (!rel_tol && iterative_mode && prec)
532  {
533  for (i = 0; i < max_iter; i++)
534  {
535  oper->Mult(x, r); // r = A x
536  subtract(b, r, r); // r = b - A x
537  prec->Mult(r, z); // z = B r
538  add(x, 1.0, z, x); // x = x + B (b - A x)
539  }
540  converged = true;
541  final_iter = i;
542  return;
543  }
544 
545  // Optimized preconditioned SLI with fixed number of iterations and zero
546  // initial guess
547  if (!rel_tol && !iterative_mode && prec)
548  {
549  prec->Mult(b, x); // x = B b (initial guess 0)
550  for (i = 1; i < max_iter; i++)
551  {
552  oper->Mult(x, r); // r = A x
553  subtract(b, r, r); // r = b - A x
554  prec->Mult(r, z); // z = B r
555  add(x, 1.0, z, x); // x = x + B (b - A x)
556  }
557  converged = true;
558  final_iter = i;
559  return;
560  }
561 
562  // General version of SLI with a relative tolerance, optional preconditioner
563  // and optional initial guess
564  double r0, nom, nom0, nomold = 1, cf;
565 
566  if (iterative_mode)
567  {
568  oper->Mult(x, r);
569  subtract(b, r, r); // r = b - A x
570  }
571  else
572  {
573  r = b;
574  x = 0.0;
575  }
576 
577  if (prec)
578  {
579  prec->Mult(r, z); // z = B r
580  nom0 = nom = sqrt(Dot(z, z));
581  }
582  else
583  {
584  nom0 = nom = sqrt(Dot(r, r));
585  }
586 
588  {
589  mfem::out << " Iteration : " << setw(3) << right << 0 << " ||Br|| = "
590  << nom << (print_options.first_and_last ? " ..." : "") << '\n';
591  }
592 
593  r0 = std::max(nom*rel_tol, abs_tol);
594  if (nom <= r0)
595  {
596  converged = true;
597  final_iter = 0;
598  final_norm = nom;
599  return;
600  }
601 
602  // start iteration
603  converged = false;
605  for (i = 1; true; )
606  {
607  if (prec) // x = x + B (b - A x)
608  {
609  add(x, 1.0, z, x);
610  }
611  else
612  {
613  add(x, 1.0, r, x);
614  }
615 
616  oper->Mult(x, r);
617  subtract(b, r, r); // r = b - A x
618 
619  if (prec)
620  {
621  prec->Mult(r, z); // z = B r
622  nom = sqrt(Dot(z, z));
623  }
624  else
625  {
626  nom = sqrt(Dot(r, r));
627  }
628 
629  cf = nom/nomold;
630  nomold = nom;
631 
632  bool done = false;
633  if (nom < r0)
634  {
635  converged = true;
636  final_iter = i;
637  done = true;
638  }
639 
640  if (++i > max_iter)
641  {
642  done = true;
643  }
644 
646  {
647  mfem::out << " Iteration : " << setw(3) << right << (i-1)
648  << " ||Br|| = " << setw(11) << left << nom
649  << "\tConv. rate: " << cf << '\n';
650  }
651 
652  if (done) { break; }
653  }
654 
656  {
657  const auto rf = pow (nom/nom0, 1.0/final_iter);
658  mfem::out << "SLI: Number of iterations: " << final_iter << '\n'
659  << "Conv. rate: " << cf << '\n'
660  << "Average reduction factor: "<< rf << '\n';
661  }
663  {
664  mfem::out << "SLI: No convergence!" << '\n';
665  }
666 
667  final_norm = nom;
668 }
669 
670 void SLI(const Operator &A, const Vector &b, Vector &x,
671  int print_iter, int max_num_iter,
672  double RTOLERANCE, double ATOLERANCE)
673 {
674  MFEM_PERF_FUNCTION;
675 
676  SLISolver sli;
677  sli.SetPrintLevel(print_iter);
678  sli.SetMaxIter(max_num_iter);
679  sli.SetRelTol(sqrt(RTOLERANCE));
680  sli.SetAbsTol(sqrt(ATOLERANCE));
681  sli.SetOperator(A);
682  sli.Mult(b, x);
683 }
684 
685 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
686  int print_iter, int max_num_iter,
687  double RTOLERANCE, double ATOLERANCE)
688 {
689  MFEM_PERF_FUNCTION;
690 
691  SLISolver sli;
692  sli.SetPrintLevel(print_iter);
693  sli.SetMaxIter(max_num_iter);
694  sli.SetRelTol(sqrt(RTOLERANCE));
695  sli.SetAbsTol(sqrt(ATOLERANCE));
696  sli.SetOperator(A);
697  sli.SetPreconditioner(B);
698  sli.Mult(b, x);
699 }
700 
701 
703 {
705 
706  r.SetSize(width, mt); r.UseDevice(true);
707  d.SetSize(width, mt); d.UseDevice(true);
708  z.SetSize(width, mt); z.UseDevice(true);
709 }
710 
711 void CGSolver::Mult(const Vector &b, Vector &x) const
712 {
713  int i;
714  double r0, den, nom, nom0, betanom, alpha, beta;
715 
716  x.UseDevice(true);
717  if (iterative_mode)
718  {
719  oper->Mult(x, r);
720  subtract(b, r, r); // r = b - A x
721  }
722  else
723  {
724  r = b;
725  x = 0.0;
726  }
727 
728  if (prec)
729  {
730  prec->Mult(r, z); // z = B r
731  d = z;
732  }
733  else
734  {
735  d = r;
736  }
737  nom0 = nom = Dot(d, r);
738  MFEM_ASSERT(IsFinite(nom), "nom = " << nom);
740  {
741  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
742  << nom << (print_options.first_and_last ? " ...\n" : "\n");
743  }
744  Monitor(0, nom, r, x);
745 
746  if (nom < 0.0)
747  {
749  {
750  mfem::out << "PCG: The preconditioner is not positive definite. (Br, r) = "
751  << nom << '\n';
752  }
753  converged = false;
754  final_iter = 0;
755  final_norm = nom;
756  return;
757  }
758  r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
759  if (nom <= r0)
760  {
761  converged = true;
762  final_iter = 0;
763  final_norm = sqrt(nom);
764  return;
765  }
766 
767  oper->Mult(d, z); // z = A d
768  den = Dot(z, d);
769  MFEM_ASSERT(IsFinite(den), "den = " << den);
770  if (den <= 0.0)
771  {
772  if (Dot(d, d) > 0.0 && print_options.warnings)
773  {
774  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
775  << den << '\n';
776  }
777  if (den == 0.0)
778  {
779  converged = false;
780  final_iter = 0;
781  final_norm = sqrt(nom);
782  return;
783  }
784  }
785 
786  // start iteration
787  converged = false;
789  for (i = 1; true; )
790  {
791  alpha = nom/den;
792  add(x, alpha, d, x); // x = x + alpha d
793  add(r, -alpha, z, r); // r = r - alpha A d
794 
795  if (prec)
796  {
797  prec->Mult(r, z); // z = B r
798  betanom = Dot(r, z);
799  }
800  else
801  {
802  betanom = Dot(r, r);
803  }
804  MFEM_ASSERT(IsFinite(betanom), "betanom = " << betanom);
805  if (betanom < 0.0)
806  {
808  {
809  mfem::out << "PCG: The preconditioner is not positive definite. (Br, r) = "
810  << betanom << '\n';
811  }
812  converged = false;
813  final_iter = i;
814  break;
815  }
816 
818  {
819  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
820  << betanom << '\n';
821  }
822 
823  Monitor(i, betanom, r, x);
824 
825  if (betanom <= r0)
826  {
827  converged = true;
828  final_iter = i;
829  break;
830  }
831 
832  if (++i > max_iter)
833  {
834  break;
835  }
836 
837  beta = betanom/nom;
838  if (prec)
839  {
840  add(z, beta, d, d); // d = z + beta d
841  }
842  else
843  {
844  add(r, beta, d, d);
845  }
846  oper->Mult(d, z); // z = A d
847  den = Dot(d, z);
848  MFEM_ASSERT(IsFinite(den), "den = " << den);
849  if (den <= 0.0)
850  {
851  if (Dot(d, d) > 0.0 && print_options.warnings)
852  {
853  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
854  << den << '\n';
855  }
856  if (den == 0.0)
857  {
858  final_iter = i;
859  break;
860  }
861  }
862  nom = betanom;
863  }
865  {
866  mfem::out << " Iteration : " << setw(3) << final_iter << " (B r, r) = "
867  << betanom << '\n';
868  }
870  {
871  mfem::out << "PCG: Number of iterations: " << final_iter << '\n';
872  }
875  {
876  const auto arf = pow (betanom/nom0, 0.5/final_iter);
877  mfem::out << "Average reduction factor = " << arf << '\n';
878  }
880  {
881  mfem::out << "PCG: No convergence!" << '\n';
882  }
883 
884  final_norm = sqrt(betanom);
885 
886  Monitor(final_iter, final_norm, r, x, true);
887 }
888 
889 void CG(const Operator &A, const Vector &b, Vector &x,
890  int print_iter, int max_num_iter,
891  double RTOLERANCE, double ATOLERANCE)
892 {
893  MFEM_PERF_FUNCTION;
894 
895  CGSolver cg;
896  cg.SetPrintLevel(print_iter);
897  cg.SetMaxIter(max_num_iter);
898  cg.SetRelTol(sqrt(RTOLERANCE));
899  cg.SetAbsTol(sqrt(ATOLERANCE));
900  cg.SetOperator(A);
901  cg.Mult(b, x);
902 }
903 
904 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
905  int print_iter, int max_num_iter,
906  double RTOLERANCE, double ATOLERANCE)
907 {
908  MFEM_PERF_FUNCTION;
909 
910  CGSolver pcg;
911  pcg.SetPrintLevel(print_iter);
912  pcg.SetMaxIter(max_num_iter);
913  pcg.SetRelTol(sqrt(RTOLERANCE));
914  pcg.SetAbsTol(sqrt(ATOLERANCE));
915  pcg.SetOperator(A);
916  pcg.SetPreconditioner(B);
917  pcg.Mult(b, x);
918 }
919 
920 
921 inline void GeneratePlaneRotation(double &dx, double &dy,
922  double &cs, double &sn)
923 {
924  if (dy == 0.0)
925  {
926  cs = 1.0;
927  sn = 0.0;
928  }
929  else if (fabs(dy) > fabs(dx))
930  {
931  double temp = dx / dy;
932  sn = 1.0 / sqrt( 1.0 + temp*temp );
933  cs = temp * sn;
934  }
935  else
936  {
937  double temp = dy / dx;
938  cs = 1.0 / sqrt( 1.0 + temp*temp );
939  sn = temp * cs;
940  }
941 }
942 
943 inline void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
944 {
945  double temp = cs * dx + sn * dy;
946  dy = -sn * dx + cs * dy;
947  dx = temp;
948 }
949 
950 inline void Update(Vector &x, int k, DenseMatrix &h, Vector &s,
951  Array<Vector*> &v)
952 {
953  Vector y(s);
954 
955  // Backsolve:
956  for (int i = k; i >= 0; i--)
957  {
958  y(i) /= h(i,i);
959  for (int j = i - 1; j >= 0; j--)
960  {
961  y(j) -= h(j,i) * y(i);
962  }
963  }
964 
965  for (int j = 0; j <= k; j++)
966  {
967  x.Add(y(j), *v[j]);
968  }
969 }
970 
971 void GMRESSolver::Mult(const Vector &b, Vector &x) const
972 {
973  // Generalized Minimum Residual method following the algorithm
974  // on p. 20 of the SIAM Templates book.
975 
976  int n = width;
977 
978  DenseMatrix H(m+1, m);
979  Vector s(m+1), cs(m+1), sn(m+1);
980  Vector r(n), w(n);
981  Array<Vector *> v;
982 
983  double resid;
984  int i, j, k;
985 
986  if (iterative_mode)
987  {
988  oper->Mult(x, r);
989  }
990  else
991  {
992  x = 0.0;
993  }
994 
995  if (prec)
996  {
997  if (iterative_mode)
998  {
999  subtract(b, r, w);
1000  prec->Mult(w, r); // r = M (b - A x)
1001  }
1002  else
1003  {
1004  prec->Mult(b, r);
1005  }
1006  }
1007  else
1008  {
1009  if (iterative_mode)
1010  {
1011  subtract(b, r, r);
1012  }
1013  else
1014  {
1015  r = b;
1016  }
1017  }
1018  double beta = Norm(r); // beta = ||r||
1019  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1020 
1021  final_norm = std::max(rel_tol*beta, abs_tol);
1022 
1023  if (beta <= final_norm)
1024  {
1025  final_norm = beta;
1026  final_iter = 0;
1027  converged = true;
1028  goto finish;
1029  }
1030 
1032  {
1033  mfem::out << " Pass : " << setw(2) << 1
1034  << " Iteration : " << setw(3) << 0
1035  << " ||B r|| = " << beta
1036  << (print_options.first_and_last ? " ...\n" : "\n");
1037  }
1038 
1039  Monitor(0, beta, r, x);
1040 
1041  v.SetSize(m+1, NULL);
1042 
1043  for (j = 1; j <= max_iter; )
1044  {
1045  if (v[0] == NULL) { v[0] = new Vector(n); }
1046  v[0]->Set(1.0/beta, r);
1047  s = 0.0; s(0) = beta;
1048 
1049  for (i = 0; i < m && j <= max_iter; i++, j++)
1050  {
1051  if (prec)
1052  {
1053  oper->Mult(*v[i], r);
1054  prec->Mult(r, w); // w = M A v[i]
1055  }
1056  else
1057  {
1058  oper->Mult(*v[i], w);
1059  }
1060 
1061  for (k = 0; k <= i; k++)
1062  {
1063  H(k,i) = Dot(w, *v[k]); // H(k,i) = w * v[k]
1064  w.Add(-H(k,i), *v[k]); // w -= H(k,i) * v[k]
1065  }
1066 
1067  H(i+1,i) = Norm(w); // H(i+1,i) = ||w||
1068  MFEM_ASSERT(IsFinite(H(i+1,i)), "Norm(w) = " << H(i+1,i));
1069  if (v[i+1] == NULL) { v[i+1] = new Vector(n); }
1070  v[i+1]->Set(1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
1071 
1072  for (k = 0; k < i; k++)
1073  {
1074  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
1075  }
1076 
1077  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1078  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1079  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
1080 
1081  resid = fabs(s(i+1));
1082  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1083 
1084  if (resid <= final_norm)
1085  {
1086  Update(x, i, H, s, v);
1087  final_norm = resid;
1088  final_iter = j;
1089  converged = true;
1090  goto finish;
1091  }
1092 
1094  {
1095  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1096  << " Iteration : " << setw(3) << j
1097  << " ||B r|| = " << resid << '\n';
1098  }
1099 
1100  Monitor(j, resid, r, x);
1101  }
1102 
1103  if (print_options.iterations && j <= max_iter)
1104  {
1105  mfem::out << "Restarting..." << '\n';
1106  }
1107 
1108  Update(x, i-1, H, s, v);
1109 
1110  oper->Mult(x, r);
1111  if (prec)
1112  {
1113  subtract(b, r, w);
1114  prec->Mult(w, r); // r = M (b - A x)
1115  }
1116  else
1117  {
1118  subtract(b, r, r);
1119  }
1120  beta = Norm(r); // beta = ||r||
1121  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1122  if (beta <= final_norm)
1123  {
1124  final_norm = beta;
1125  final_iter = j;
1126  converged = true;
1127  goto finish;
1128  }
1129  }
1130 
1131  final_norm = beta;
1132  final_iter = max_iter;
1133  converged = false;
1134 
1135 finish:
1137  {
1138  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1139  << " Iteration : " << setw(3) << final_iter
1140  << " ||B r|| = " << resid << '\n';
1141  }
1143  {
1144  mfem::out << "GMRES: Number of iterations: " << final_iter << '\n';
1145  }
1147  {
1148  mfem::out << "GMRES: No convergence!\n";
1149  }
1150 
1151  Monitor(final_iter, final_norm, r, x, true);
1152 
1153  for (i = 0; i < v.Size(); i++)
1154  {
1155  delete v[i];
1156  }
1157 }
1158 
1159 void FGMRESSolver::Mult(const Vector &b, Vector &x) const
1160 {
1161  DenseMatrix H(m+1,m);
1162  Vector s(m+1), cs(m+1), sn(m+1);
1163  Vector r(b.Size());
1164 
1165  int i, j, k;
1166 
1167 
1168  if (iterative_mode)
1169  {
1170  oper->Mult(x, r);
1171  subtract(b,r,r);
1172  }
1173  else
1174  {
1175  x = 0.;
1176  r = b;
1177  }
1178  double beta = Norm(r); // beta = ||r||
1179  // We need to preallocate this to report the correct result in the case of
1180  // no convergence.
1181  double resid;
1182  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1183 
1184  final_norm = std::max(rel_tol*beta, abs_tol);
1185 
1186  if (beta <= final_norm)
1187  {
1188  final_norm = beta;
1189  final_iter = 0;
1190  converged = true;
1191  return;
1192  }
1193 
1195  {
1196  mfem::out << " Pass : " << setw(2) << 1
1197  << " Iteration : " << setw(3) << 0
1198  << " || r || = " << beta
1199  << (print_options.first_and_last ? " ...\n" : "\n");
1200  }
1201 
1202  Monitor(0, beta, r, x);
1203 
1204  Array<Vector*> v(m+1);
1205  Array<Vector*> z(m+1);
1206  for (i= 0; i<=m; i++)
1207  {
1208  v[i] = NULL;
1209  z[i] = NULL;
1210  }
1211 
1212  j = 1;
1213  while (j <= max_iter)
1214  {
1215  if (v[0] == NULL) { v[0] = new Vector(b.Size()); }
1216  (*v[0]) = 0.0;
1217  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
1218  s = 0.0; s(0) = beta;
1219 
1220  for (i = 0; i < m && j <= max_iter; i++, j++)
1221  {
1222 
1223  if (z[i] == NULL) { z[i] = new Vector(b.Size()); }
1224  (*z[i]) = 0.0;
1225 
1226  if (prec)
1227  {
1228  prec->Mult(*v[i], *z[i]);
1229  }
1230  else
1231  {
1232  (*z[i]) = (*v[i]);
1233  }
1234  oper->Mult(*z[i], r);
1235 
1236  for (k = 0; k <= i; k++)
1237  {
1238  H(k,i) = Dot( r, *v[k]); // H(k,i) = r * v[k]
1239  r.Add(-H(k,i), (*v[k])); // r -= H(k,i) * v[k]
1240  }
1241 
1242  H(i+1,i) = Norm(r); // H(i+1,i) = ||r||
1243  if (v[i+1] == NULL) { v[i+1] = new Vector(b.Size()); }
1244  (*v[i+1]) = 0.0;
1245  v[i+1] -> Add (1.0/H(i+1,i), r); // v[i+1] = r / H(i+1,i)
1246 
1247  for (k = 0; k < i; k++)
1248  {
1249  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
1250  }
1251 
1252  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1253  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1254  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
1255 
1256  resid = fabs(s(i+1));
1257  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1259  resid <= final_norm))
1260  {
1261  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1262  << " Iteration : " << setw(3) << j
1263  << " || r || = " << resid << endl;
1264  }
1265  Monitor(j, resid, r, x, resid <= final_norm);
1266 
1267  if (resid <= final_norm)
1268  {
1269  Update(x, i, H, s, z);
1270  final_norm = resid;
1271  final_iter = j;
1272  converged = true;
1273 
1274  if (print_options.summary)
1275  {
1276  mfem::out << "FGMRES: Number of iterations: " << final_iter << '\n';
1277  }
1278 
1279  for (i= 0; i<=m; i++)
1280  {
1281  if (v[i]) { delete v[i]; }
1282  if (z[i]) { delete z[i]; }
1283  }
1284  return;
1285  }
1286  }
1287 
1289  {
1290  mfem::out << "Restarting..." << endl;
1291  }
1292 
1293  Update(x, i-1, H, s, z);
1294 
1295  oper->Mult(x, r);
1296  subtract(b,r,r);
1297  beta = Norm(r);
1298  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1299  if (beta <= final_norm)
1300  {
1301  final_norm = beta;
1302  final_iter = j;
1303  converged = true;
1304 
1305  break;
1306  }
1307  }
1308 
1309  // Clean buffers up
1310  for (i = 0; i <= m; i++)
1311  {
1312  if (v[i]) { delete v[i]; }
1313  if (z[i]) { delete z[i]; }
1314  }
1315  converged = false;
1316 
1317  // Note: j is off by one when we arrive here
1319  {
1320  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1321  << " Iteration : " << setw(3) << j-1
1322  << " || r || = " << resid << endl;
1323  }
1325  {
1326  mfem::out << "FGMRES: Number of iterations: " << j-1 << '\n';
1327  }
1329  {
1330  mfem::out << "FGMRES: No convergence!\n";
1331  }
1332 }
1333 
1334 
1335 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
1336  int &max_iter, int m, double &tol, double atol, int printit)
1337 {
1338  MFEM_PERF_FUNCTION;
1339 
1340  GMRESSolver gmres;
1341  gmres.SetPrintLevel(printit);
1342  gmres.SetMaxIter(max_iter);
1343  gmres.SetKDim(m);
1344  gmres.SetRelTol(sqrt(tol));
1345  gmres.SetAbsTol(sqrt(atol));
1346  gmres.SetOperator(A);
1347  gmres.SetPreconditioner(M);
1348  gmres.Mult(b, x);
1349  max_iter = gmres.GetNumIterations();
1350  tol = gmres.GetFinalNorm()*gmres.GetFinalNorm();
1351  return gmres.GetConverged();
1352 }
1353 
1354 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
1355  int print_iter, int max_num_iter, int m, double rtol, double atol)
1356 {
1357  GMRES(A, x, b, B, max_num_iter, m, rtol, atol, print_iter);
1358 }
1359 
1360 
1362 {
1363  p.SetSize(width);
1364  phat.SetSize(width);
1365  s.SetSize(width);
1366  shat.SetSize(width);
1367  t.SetSize(width);
1368  v.SetSize(width);
1369  r.SetSize(width);
1370  rtilde.SetSize(width);
1371 }
1372 
1373 void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
1374 {
1375  // BiConjugate Gradient Stabilized method following the algorithm
1376  // on p. 27 of the SIAM Templates book.
1377 
1378  int i;
1379  double resid, tol_goal;
1380  double rho_1, rho_2=1.0, alpha=1.0, beta, omega=1.0;
1381 
1382  if (iterative_mode)
1383  {
1384  oper->Mult(x, r);
1385  subtract(b, r, r); // r = b - A x
1386  }
1387  else
1388  {
1389  x = 0.0;
1390  r = b;
1391  }
1392  rtilde = r;
1393 
1394  resid = Norm(r);
1395  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1397  {
1398  mfem::out << " Iteration : " << setw(3) << 0
1399  << " ||r|| = " << resid << (print_options.first_and_last ? " ...\n" : "\n");
1400  }
1401 
1402  Monitor(0, resid, r, x);
1403 
1404  tol_goal = std::max(resid*rel_tol, abs_tol);
1405 
1406  if (resid <= tol_goal)
1407  {
1408  final_norm = resid;
1409  final_iter = 0;
1410  converged = true;
1411  return;
1412  }
1413 
1414  for (i = 1; i <= max_iter; i++)
1415  {
1416  rho_1 = Dot(rtilde, r);
1417  if (rho_1 == 0)
1418  {
1420  {
1421  mfem::out << " Iteration : " << setw(3) << i
1422  << " ||r|| = " << resid << '\n';
1423  }
1424 
1425  Monitor(i, resid, r, x);
1426 
1427  final_norm = resid;
1428  final_iter = i;
1429  converged = false;
1431  {
1432  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1433  }
1434  if (print_options.warnings)
1435  {
1436  mfem::out << "BiCGStab: No convergence!\n";
1437  }
1438  return;
1439  }
1440  if (i == 1)
1441  {
1442  p = r;
1443  }
1444  else
1445  {
1446  beta = (rho_1/rho_2) * (alpha/omega);
1447  add(p, -omega, v, p); // p = p - omega * v
1448  add(r, beta, p, p); // p = r + beta * p
1449  }
1450  if (prec)
1451  {
1452  prec->Mult(p, phat); // phat = M^{-1} * p
1453  }
1454  else
1455  {
1456  phat = p;
1457  }
1458  oper->Mult(phat, v); // v = A * phat
1459  alpha = rho_1 / Dot(rtilde, v);
1460  add(r, -alpha, v, s); // s = r - alpha * v
1461  resid = Norm(s);
1462  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1463  if (resid < tol_goal)
1464  {
1465  x.Add(alpha, phat); // x = x + alpha * phat
1467  {
1468  mfem::out << " Iteration : " << setw(3) << i
1469  << " ||s|| = " << resid << '\n';
1470  }
1471  final_norm = resid;
1472  final_iter = i;
1473  converged = true;
1475  {
1476  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1477  }
1478  return;
1479  }
1481  {
1482  mfem::out << " Iteration : " << setw(3) << i
1483  << " ||s|| = " << resid;
1484  }
1485  Monitor(i, resid, r, x);
1486  if (prec)
1487  {
1488  prec->Mult(s, shat); // shat = M^{-1} * s
1489  }
1490  else
1491  {
1492  shat = s;
1493  }
1494  oper->Mult(shat, t); // t = A * shat
1495  omega = Dot(t, s) / Dot(t, t);
1496  x.Add(alpha, phat); // x += alpha * phat
1497  x.Add(omega, shat); // x += omega * shat
1498  add(s, -omega, t, r); // r = s - omega * t
1499 
1500  rho_2 = rho_1;
1501  resid = Norm(r);
1502  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1504  {
1505  mfem::out << " ||r|| = " << resid << '\n';
1506  }
1507  Monitor(i, resid, r, x);
1508  if (resid < tol_goal)
1509  {
1510  final_norm = resid;
1511  final_iter = i;
1512  converged = true;
1514  {
1515  mfem::out << " Iteration : " << setw(3) << i
1516  << " ||r|| = " << resid << '\n';
1517  }
1519  {
1520  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1521  }
1522  return;
1523  }
1524  if (omega == 0)
1525  {
1526  final_norm = resid;
1527  final_iter = i;
1528  converged = false;
1530  {
1531  mfem::out << " Iteration : " << setw(3) << i
1532  << " ||r|| = " << resid << '\n';
1533  }
1535  {
1536  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1537  }
1538  if (print_options.warnings)
1539  {
1540  mfem::out << "BiCGStab: No convergence!\n";
1541  }
1542  return;
1543  }
1544  }
1545 
1546  final_norm = resid;
1547  final_iter = max_iter;
1548  converged = false;
1549 
1551  {
1552  mfem::out << " Iteration : " << setw(3) << final_iter
1553  << " ||r|| = " << resid << '\n';
1554  }
1556  {
1557  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1558  }
1559  if (print_options.warnings)
1560  {
1561  mfem::out << "BiCGStab: No convergence!\n";
1562  }
1563 }
1564 
1565 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
1566  int &max_iter, double &tol, double atol, int printit)
1567 {
1568  BiCGSTABSolver bicgstab;
1569  bicgstab.SetPrintLevel(printit);
1570  bicgstab.SetMaxIter(max_iter);
1571  bicgstab.SetRelTol(sqrt(tol));
1572  bicgstab.SetAbsTol(sqrt(atol));
1573  bicgstab.SetOperator(A);
1574  bicgstab.SetPreconditioner(M);
1575  bicgstab.Mult(b, x);
1576  max_iter = bicgstab.GetNumIterations();
1577  tol = bicgstab.GetFinalNorm()*bicgstab.GetFinalNorm();
1578  return bicgstab.GetConverged();
1579 }
1580 
1581 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
1582  int print_iter, int max_num_iter, double rtol, double atol)
1583 {
1584  BiCGSTAB(A, x, b, B, max_num_iter, rtol, atol, print_iter);
1585 }
1586 
1587 
1589 {
1591  v0.SetSize(width);
1592  v1.SetSize(width);
1593  w0.SetSize(width);
1594  w1.SetSize(width);
1595  q.SetSize(width);
1596  if (prec)
1597  {
1598  u1.SetSize(width);
1599  }
1600 }
1601 
1602 void MINRESSolver::Mult(const Vector &b, Vector &x) const
1603 {
1604  // Based on the MINRES algorithm on p. 86, Fig. 6.9 in
1605  // "Iterative Krylov Methods for Large Linear Systems",
1606  // by Henk A. van der Vorst, 2003.
1607  // Extended to support an SPD preconditioner.
1608 
1609  int it;
1610  double beta, eta, gamma0, gamma1, sigma0, sigma1;
1611  double alpha, delta, rho1, rho2, rho3, norm_goal;
1612  Vector *z = (prec) ? &u1 : &v1;
1613 
1614  converged = true;
1615 
1616  if (!iterative_mode)
1617  {
1618  v1 = b;
1619  x = 0.;
1620  }
1621  else
1622  {
1623  oper->Mult(x, v1);
1624  subtract(b, v1, v1);
1625  }
1626 
1627  if (prec)
1628  {
1629  prec->Mult(v1, u1);
1630  }
1631  eta = beta = sqrt(Dot(*z, v1));
1632  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1633  gamma0 = gamma1 = 1.;
1634  sigma0 = sigma1 = 0.;
1635 
1636  norm_goal = std::max(rel_tol*eta, abs_tol);
1637 
1638  if (eta <= norm_goal)
1639  {
1640  it = 0;
1641  goto loop_end;
1642  }
1643 
1645  {
1646  mfem::out << "MINRES: iteration " << setw(3) << 0 << ": ||r||_B = "
1647  << eta << (print_options.first_and_last ? " ..." : "") << '\n';
1648  }
1649  Monitor(0, eta, *z, x);
1650 
1651  for (it = 1; it <= max_iter; it++)
1652  {
1653  v1 /= beta;
1654  if (prec)
1655  {
1656  u1 /= beta;
1657  }
1658  oper->Mult(*z, q);
1659  alpha = Dot(*z, q);
1660  MFEM_ASSERT(IsFinite(alpha), "alpha = " << alpha);
1661  if (it > 1) // (v0 == 0) for (it == 1)
1662  {
1663  q.Add(-beta, v0);
1664  }
1665  add(q, -alpha, v1, v0);
1666 
1667  delta = gamma1*alpha - gamma0*sigma1*beta;
1668  rho3 = sigma0*beta;
1669  rho2 = sigma1*alpha + gamma0*gamma1*beta;
1670  if (!prec)
1671  {
1672  beta = Norm(v0);
1673  }
1674  else
1675  {
1676  prec->Mult(v0, q);
1677  beta = sqrt(Dot(v0, q));
1678  }
1679  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1680  rho1 = hypot(delta, beta);
1681 
1682  if (it == 1)
1683  {
1684  w0.Set(1./rho1, *z); // (w0 == 0) and (w1 == 0)
1685  }
1686  else if (it == 2)
1687  {
1688  add(1./rho1, *z, -rho2/rho1, w1, w0); // (w0 == 0)
1689  }
1690  else
1691  {
1692  add(-rho3/rho1, w0, -rho2/rho1, w1, w0);
1693  w0.Add(1./rho1, *z);
1694  }
1695 
1696  gamma0 = gamma1;
1697  gamma1 = delta/rho1;
1698 
1699  x.Add(gamma1*eta, w0);
1700 
1701  sigma0 = sigma1;
1702  sigma1 = beta/rho1;
1703 
1704  eta = -sigma1*eta;
1705  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1706 
1707  if (fabs(eta) <= norm_goal)
1708  {
1709  goto loop_end;
1710  }
1711 
1713  {
1714  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1715  << fabs(eta) << '\n';
1716  }
1717  Monitor(it, fabs(eta), *z, x);
1718 
1719  if (prec)
1720  {
1721  Swap(u1, q);
1722  }
1723  Swap(v0, v1);
1724  Swap(w0, w1);
1725  }
1726  converged = false;
1727  it--;
1728 
1729 loop_end:
1730  final_iter = it;
1731  final_norm = fabs(eta);
1732 
1734  {
1735  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1736  << fabs(eta) << '\n';
1737  }
1738 
1740  {
1741  mfem::out << "MINRES: Number of iterations: " << setw(3) << final_iter << '\n';
1742  }
1743 
1744  Monitor(final_iter, final_norm, *z, x, true);
1745 
1746  // if (print_options.iteration_details || (!converged && print_options.errors))
1747  // {
1748  // oper->Mult(x, v1);
1749  // subtract(b, v1, v1);
1750  // if (prec)
1751  // {
1752  // prec->Mult(v1, u1);
1753  // }
1754  // eta = sqrt(Dot(*z, v1));
1755  // mfem::out << "MINRES: iteration " << setw(3) << it << '\n'
1756  // << " ||r||_B = " << eta << " (re-computed)" << '\n';
1757  // }
1758 
1759  if (!converged && (print_options.warnings))
1760  {
1761  mfem::out << "MINRES: No convergence!\n";
1762  }
1763 }
1764 
1765 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it,
1766  int max_it, double rtol, double atol)
1767 {
1768  MFEM_PERF_FUNCTION;
1769 
1770  MINRESSolver minres;
1771  minres.SetPrintLevel(print_it);
1772  minres.SetMaxIter(max_it);
1773  minres.SetRelTol(sqrt(rtol));
1774  minres.SetAbsTol(sqrt(atol));
1775  minres.SetOperator(A);
1776  minres.Mult(b, x);
1777 }
1778 
1779 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
1780  int print_it, int max_it, double rtol, double atol)
1781 {
1782  MINRESSolver minres;
1783  minres.SetPrintLevel(print_it);
1784  minres.SetMaxIter(max_it);
1785  minres.SetRelTol(sqrt(rtol));
1786  minres.SetAbsTol(sqrt(atol));
1787  minres.SetOperator(A);
1788  minres.SetPreconditioner(B);
1789  minres.Mult(b, x);
1790 }
1791 
1792 
1794 {
1795  oper = &op;
1796  height = op.Height();
1797  width = op.Width();
1798  MFEM_ASSERT(height == width, "square Operator is required.");
1799 
1800  r.SetSize(width);
1801  c.SetSize(width);
1802 }
1803 
1804 void NewtonSolver::Mult(const Vector &b, Vector &x) const
1805 {
1806  MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
1807  MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver).");
1808 
1809  int it;
1810  double norm0, norm, norm_goal;
1811  const bool have_b = (b.Size() == Height());
1812 
1813  if (!iterative_mode)
1814  {
1815  x = 0.0;
1816  }
1817 
1818  ProcessNewState(x);
1819 
1820  oper->Mult(x, r);
1821  if (have_b)
1822  {
1823  r -= b;
1824  }
1825 
1826  norm0 = norm = Norm(r);
1828  {
1829  mfem::out << "Newton iteration " << setw(2) << 0
1830  << " : ||r|| = " << norm << "...\n";
1831  }
1832  norm_goal = std::max(rel_tol*norm, abs_tol);
1833 
1834  prec->iterative_mode = false;
1835 
1836  // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b]
1837  for (it = 0; true; it++)
1838  {
1839  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
1841  {
1842  mfem::out << "Newton iteration " << setw(2) << it
1843  << " : ||r|| = " << norm;
1844  if (it > 0)
1845  {
1846  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
1847  }
1848  mfem::out << '\n';
1849  }
1850  Monitor(it, norm, r, x);
1851 
1852  if (norm <= norm_goal)
1853  {
1854  converged = true;
1855  break;
1856  }
1857 
1858  if (it >= max_iter)
1859  {
1860  converged = false;
1861  break;
1862  }
1863 
1864  grad = &oper->GetGradient(x);
1865  prec->SetOperator(*grad);
1866 
1867  if (lin_rtol_type)
1868  {
1869  AdaptiveLinRtolPreSolve(x, it, norm);
1870  }
1871 
1872  prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b]
1873 
1874  if (lin_rtol_type)
1875  {
1876  AdaptiveLinRtolPostSolve(c, r, it, norm);
1877  }
1878 
1879  const double c_scale = ComputeScalingFactor(x, b);
1880  if (c_scale == 0.0)
1881  {
1882  converged = false;
1883  break;
1884  }
1885  add(x, -c_scale, c, x);
1886 
1887  ProcessNewState(x);
1888 
1889  oper->Mult(x, r);
1890  if (have_b)
1891  {
1892  r -= b;
1893  }
1894  norm = Norm(r);
1895  }
1896 
1897  final_iter = it;
1898  final_norm = norm;
1899 
1902  {
1903  mfem::out << "Newton: Number of iterations: " << final_iter << '\n'
1904  << " ||r|| = " << final_norm << '\n';
1905  }
1907  {
1908  mfem::out << "Newton: No convergence!\n";
1909  }
1910 }
1911 
1913  const double rtol0,
1914  const double rtol_max,
1915  const double alpha_,
1916  const double gamma_)
1917 {
1918  lin_rtol_type = type;
1919  lin_rtol0 = rtol0;
1920  lin_rtol_max = rtol_max;
1921  this->alpha = alpha_;
1922  this->gamma = gamma_;
1923 }
1924 
1926  const int it,
1927  const double fnorm) const
1928 {
1929  // Assume that when adaptive linear solver relative tolerance is activated,
1930  // we are working with an iterative solver.
1931  auto iterative_solver = static_cast<IterativeSolver *>(prec);
1932  // Adaptive linear solver relative tolerance
1933  double eta;
1934  // Safeguard threshold
1935  double sg_threshold = 0.1;
1936 
1937  if (it == 0)
1938  {
1939  eta = lin_rtol0;
1940  }
1941  else
1942  {
1943  if (lin_rtol_type == 1)
1944  {
1945  // eta = gamma * abs(||F(x1)|| - ||F(x0) + DF(x0) s0||) / ||F(x0)||
1946  eta = gamma * abs(fnorm - lnorm_last) / fnorm_last;
1947  }
1948  else if (lin_rtol_type == 2)
1949  {
1950  // eta = gamma * (||F(x1)|| / ||F(x0)||)^alpha
1951  eta = gamma * pow(fnorm / fnorm_last, alpha);
1952  }
1953  else
1954  {
1955  MFEM_ABORT("Unknown adaptive linear solver rtol version");
1956  }
1957 
1958  // Safeguard rtol from "oversolving" ?!
1959  const double sg_eta = gamma * pow(eta_last, alpha);
1960  if (sg_eta > sg_threshold) { eta = std::max(eta, sg_eta); }
1961  }
1962 
1963  eta = std::min(eta, lin_rtol_max);
1964  iterative_solver->SetRelTol(eta);
1965  eta_last = eta;
1967  {
1968  mfem::out << "Eisenstat-Walker rtol = " << eta << "\n";
1969  }
1970 }
1971 
1973  const Vector &b,
1974  const int it,
1975  const double fnorm) const
1976 {
1977  fnorm_last = fnorm;
1978 
1979  // If version 1 is chosen, the true linear residual norm has to be computed
1980  // and in most cases we can only retrieve the preconditioned linear residual
1981  // norm.
1982  if (lin_rtol_type == 1)
1983  {
1984  // lnorm_last = ||F(x0) + DF(x0) s0||
1985  Vector linres(x.Size());
1986  grad->Mult(x, linres);
1987  linres -= b;
1988  lnorm_last = Norm(linres);
1989  }
1990 }
1991 
1992 void LBFGSSolver::Mult(const Vector &b, Vector &x) const
1993 {
1994  MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");
1995 
1996  // Quadrature points that are checked for negative Jacobians etc.
1997  Vector sk, rk, yk, rho, alpha;
1998 
1999  // r - r_{k+1}, c - descent direction
2000  sk.SetSize(width); // x_{k+1}-x_k
2001  rk.SetSize(width); // nabla(f(x_{k}))
2002  yk.SetSize(width); // r_{k+1}-r_{k}
2003  rho.SetSize(m); // 1/(dot(yk,sk)
2004  alpha.SetSize(m); // rhok*sk'*c
2005  int last_saved_id = -1;
2006 
2007  int it;
2008  double norm0, norm, norm_goal;
2009  const bool have_b = (b.Size() == Height());
2010 
2011  if (!iterative_mode)
2012  {
2013  x = 0.0;
2014  }
2015 
2016  ProcessNewState(x);
2017 
2018  // r = F(x)-b
2019  oper->Mult(x, r);
2020  if (have_b) { r -= b; }
2021 
2022  c = r; // initial descent direction
2023 
2024  norm0 = norm = Norm(r);
2026  {
2027  mfem::out << "LBFGS iteration " << setw(2) << 0
2028  << " : ||r|| = " << norm << "...\n";
2029  }
2030  norm_goal = std::max(rel_tol*norm, abs_tol);
2031  for (it = 0; true; it++)
2032  {
2033  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
2035  {
2036  mfem::out << "LBFGS iteration " << it
2037  << " : ||r|| = " << norm;
2038  if (it > 0)
2039  {
2040  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
2041  }
2042  mfem::out << '\n';
2043  }
2044 
2045  if (norm <= norm_goal)
2046  {
2047  converged = true;
2048  break;
2049  }
2050 
2051  if (it >= max_iter)
2052  {
2053  converged = false;
2054  break;
2055  }
2056 
2057  rk = r;
2058  const double c_scale = ComputeScalingFactor(x, b);
2059  if (c_scale == 0.0)
2060  {
2061  converged = false;
2062  break;
2063  }
2064  add(x, -c_scale, c, x); // x_{k+1} = x_k - c_scale*c
2065 
2066  ProcessNewState(x);
2067 
2068  oper->Mult(x, r);
2069  if (have_b)
2070  {
2071  r -= b;
2072  }
2073 
2074  // LBFGS - construct descent direction
2075  subtract(r, rk, yk); // yk = r_{k+1} - r_{k}
2076  sk = c; sk *= -c_scale; //sk = x_{k+1} - x_{k} = -c_scale*c
2077  const double gamma = Dot(sk, yk)/Dot(yk, yk);
2078 
2079  // Save last m vectors
2080  last_saved_id = (last_saved_id == m-1) ? 0 : last_saved_id+1;
2081  *skArray[last_saved_id] = sk;
2082  *ykArray[last_saved_id] = yk;
2083 
2084  c = r;
2085  for (int i = last_saved_id; i > -1; i--)
2086  {
2087  rho(i) = 1.0/Dot((*skArray[i]),(*ykArray[i]));
2088  alpha(i) = rho(i)*Dot((*skArray[i]),c);
2089  add(c, -alpha(i), (*ykArray[i]), c);
2090  }
2091  if (it > m-1)
2092  {
2093  for (int i = m-1; i > last_saved_id; i--)
2094  {
2095  rho(i) = 1./Dot((*skArray[i]), (*ykArray[i]));
2096  alpha(i) = rho(i)*Dot((*skArray[i]),c);
2097  add(c, -alpha(i), (*ykArray[i]), c);
2098  }
2099  }
2100 
2101  c *= gamma; // scale search direction
2102  if (it > m-1)
2103  {
2104  for (int i = last_saved_id+1; i < m ; i++)
2105  {
2106  double betai = rho(i)*Dot((*ykArray[i]), c);
2107  add(c, alpha(i)-betai, (*skArray[i]), c);
2108  }
2109  }
2110  for (int i = 0; i < last_saved_id+1 ; i++)
2111  {
2112  double betai = rho(i)*Dot((*ykArray[i]), c);
2113  add(c, alpha(i)-betai, (*skArray[i]), c);
2114  }
2115 
2116  norm = Norm(r);
2117  }
2118 
2119  final_iter = it;
2120  final_norm = norm;
2121 
2124  {
2125  mfem::out << "LBFGS: Number of iterations: " << final_iter << '\n'
2126  << " ||r|| = " << final_norm << '\n';
2127  }
2129  {
2130  mfem::out << "LBFGS: No convergence!\n";
2131  }
2132 }
2133 
2134 int aGMRES(const Operator &A, Vector &x, const Vector &b,
2135  const Operator &M, int &max_iter,
2136  int m_max, int m_min, int m_step, double cf,
2137  double &tol, double &atol, int printit)
2138 {
2139  int n = A.Width();
2140 
2141  int m = m_max;
2142 
2143  DenseMatrix H(m+1,m);
2144  Vector s(m+1), cs(m+1), sn(m+1);
2145  Vector w(n), av(n);
2146 
2147  double r1, resid;
2148  int i, j, k;
2149 
2150  M.Mult(b,w);
2151  double normb = w.Norml2(); // normb = ||M b||
2152  if (normb == 0.0)
2153  {
2154  normb = 1;
2155  }
2156 
2157  Vector r(n);
2158  A.Mult(x, r);
2159  subtract(b,r,w);
2160  M.Mult(w, r); // r = M (b - A x)
2161  double beta = r.Norml2(); // beta = ||r||
2162 
2163  resid = beta / normb;
2164 
2165  if (resid * resid <= tol)
2166  {
2167  tol = resid * resid;
2168  max_iter = 0;
2169  return 0;
2170  }
2171 
2172  if (printit)
2173  {
2174  mfem::out << " Pass : " << setw(2) << 1
2175  << " Iteration : " << setw(3) << 0
2176  << " (r, r) = " << beta*beta << '\n';
2177  }
2178 
2179  tol *= (normb*normb);
2180  tol = (atol > tol) ? atol : tol;
2181 
2182  m = m_max;
2183  Array<Vector *> v(m+1);
2184  for (i= 0; i<=m; i++)
2185  {
2186  v[i] = new Vector(n);
2187  (*v[i]) = 0.0;
2188  }
2189 
2190  j = 1;
2191  while (j <= max_iter)
2192  {
2193  (*v[0]) = 0.0;
2194  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
2195  s = 0.0; s(0) = beta;
2196 
2197  r1 = beta;
2198 
2199  for (i = 0; i < m && j <= max_iter; i++)
2200  {
2201  A.Mult((*v[i]),av);
2202  M.Mult(av,w); // w = M A v[i]
2203 
2204  for (k = 0; k <= i; k++)
2205  {
2206  H(k,i) = w * (*v[k]); // H(k,i) = w * v[k]
2207  w.Add(-H(k,i), (*v[k])); // w -= H(k,i) * v[k]
2208  }
2209 
2210  H(i+1,i) = w.Norml2(); // H(i+1,i) = ||w||
2211  (*v[i+1]) = 0.0;
2212  v[i+1] -> Add (1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
2213 
2214  for (k = 0; k < i; k++)
2215  {
2216  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
2217  }
2218 
2219  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
2220  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
2221  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
2222 
2223  resid = fabs(s(i+1));
2224  if (printit)
2225  {
2226  mfem::out << " Pass : " << setw(2) << j
2227  << " Iteration : " << setw(3) << i+1
2228  << " (r, r) = " << resid*resid << '\n';
2229  }
2230 
2231  if ( resid*resid < tol)
2232  {
2233  Update(x, i, H, s, v);
2234  tol = resid * resid;
2235  max_iter = j;
2236  for (i= 0; i<=m; i++)
2237  {
2238  delete v[i];
2239  }
2240  return 0;
2241  }
2242  }
2243 
2244  if (printit)
2245  {
2246  mfem::out << "Restarting..." << '\n';
2247  }
2248 
2249  Update(x, i-1, H, s, v);
2250 
2251  A.Mult(x, r);
2252  subtract(b,r,w);
2253  M.Mult(w, r); // r = M (b - A x)
2254  beta = r.Norml2(); // beta = ||r||
2255  if ( resid*resid < tol)
2256  {
2257  tol = resid * resid;
2258  max_iter = j;
2259  for (i= 0; i<=m; i++)
2260  {
2261  delete v[i];
2262  }
2263  return 0;
2264  }
2265 
2266  if (beta/r1 > cf)
2267  {
2268  if (m - m_step >= m_min)
2269  {
2270  m -= m_step;
2271  }
2272  else
2273  {
2274  m = m_max;
2275  }
2276  }
2277 
2278  j++;
2279  }
2280 
2281  tol = resid * resid;
2282  for (i= 0; i<=m; i++)
2283  {
2284  delete v[i];
2285  }
2286  return 1;
2287 }
2288 
2290  const Operator *C_,
2291  const Operator *D_)
2292  : C(C_), D(D_), c_e(NULL), d_lo(NULL), d_hi(NULL), x_lo(NULL), x_hi(NULL),
2293  input_size(insize)
2294 {
2295  if (C) { MFEM_ASSERT(C->Width() == input_size, "Wrong width of C."); }
2296  if (D) { MFEM_ASSERT(D->Width() == input_size, "Wrong width of D."); }
2297 }
2298 
2300 {
2301  MFEM_ASSERT(C, "The C operator is unspecified -- can't set constraints.");
2302  MFEM_ASSERT(c.Size() == C->Height(), "Wrong size of the constraint.");
2303 
2304  c_e = &c;
2305 }
2306 
2308  const Vector &dh)
2309 {
2310  MFEM_ASSERT(D, "The D operator is unspecified -- can't set constraints.");
2311  MFEM_ASSERT(dl.Size() == D->Height() && dh.Size() == D->Height(),
2312  "Wrong size of the constraint.");
2313 
2314  d_lo = &dl; d_hi = &dh;
2315 }
2316 
2318 {
2319  MFEM_ASSERT(xl.Size() == input_size && xh.Size() == input_size,
2320  "Wrong size of the constraint.");
2321 
2322  x_lo = &xl; x_hi = &xh;
2323 }
2324 
2326 {
2327  int m = 0;
2328  if (C) { m += C->Height(); }
2329  if (D) { m += D->Height(); }
2330  return m;
2331 }
2332 
2334 {
2335  if (print_options.warnings)
2336  {
2337  MFEM_WARNING("Objective functional is ignored as SLBQP always minimizes"
2338  "the l2 norm of (x - x_target).");
2339  }
2340  MFEM_ASSERT(prob.GetC(), "Linear constraint is not set.");
2341  MFEM_ASSERT(prob.GetC()->Height() == 1, "Solver expects scalar constraint.");
2342 
2343  problem = &prob;
2344 }
2345 
2346 void SLBQPOptimizer::SetBounds(const Vector &lo_, const Vector &hi_)
2347 {
2348  lo.SetDataAndSize(lo_.GetData(), lo_.Size());
2349  hi.SetDataAndSize(hi_.GetData(), hi_.Size());
2350 }
2351 
2353 {
2354  w.SetDataAndSize(w_.GetData(), w_.Size());
2355  a = a_;
2356 }
2357 
2358 inline void SLBQPOptimizer::print_iteration(int it, double r, double l) const
2359 {
2361  {
2362  mfem::out << "SLBQP iteration " << it << ": residual = " << r
2363  << ", lambda = " << l << '\n';
2364  }
2365 }
2366 
2367 void SLBQPOptimizer::Mult(const Vector& xt, Vector& x) const
2368 {
2369  // Based on code provided by Denis Ridzal, dridzal@sandia.gov.
2370  // Algorithm adapted from Dai and Fletcher, "New Algorithms for
2371  // Singly Linearly Constrained Quadratic Programs Subject to Lower
2372  // and Upper Bounds", Numerical Analysis Report NA/216, 2003.
2373 
2374  // Set some algorithm-specific constants and temporaries.
2375  int nclip = 0;
2376  double l = 0;
2377  double llow = 0;
2378  double lupp = 0;
2379  double lnew = 0;
2380  double dl = 2;
2381  double r = 0;
2382  double rlow = 0;
2383  double rupp = 0;
2384  double s = 0;
2385 
2386  const double smin = 0.1;
2387 
2388  const double tol = max(abs_tol, rel_tol*a);
2389 
2390  // *** Start bracketing phase of SLBQP ***
2392  {
2393  mfem::out << "SLBQP bracketing phase" << '\n';
2394  }
2395 
2396  // Solve QP with fixed Lagrange multiplier
2397  r = solve(l,xt,x,nclip);
2398  print_iteration(nclip, r, l);
2399 
2400 
2401  // If x=xt was already within bounds and satisfies the linear
2402  // constraint, then we already have the solution.
2403  if (fabs(r) <= tol)
2404  {
2405  converged = true;
2406  goto slbqp_done;
2407  }
2408 
2409  if (r < 0)
2410  {
2411  llow = l; rlow = r; l = l + dl;
2412 
2413  // Solve QP with fixed Lagrange multiplier
2414  r = solve(l,xt,x,nclip);
2415  print_iteration(nclip, r, l);
2416 
2417  while ((r < 0) && (nclip < max_iter))
2418  {
2419  llow = l;
2420  s = rlow/r - 1.0;
2421  if (s < smin) { s = smin; }
2422  dl = dl + dl/s;
2423  l = l + dl;
2424 
2425  // Solve QP with fixed Lagrange multiplier
2426  r = solve(l,xt,x,nclip);
2427  print_iteration(nclip, r, l);
2428  }
2429 
2430  lupp = l; rupp = r;
2431  }
2432  else
2433  {
2434  lupp = l; rupp = r; l = l - dl;
2435 
2436  // Solve QP with fixed Lagrange multiplier
2437  r = solve(l,xt,x,nclip);
2438  print_iteration(nclip, r, l);
2439 
2440  while ((r > 0) && (nclip < max_iter))
2441  {
2442  lupp = l;
2443  s = rupp/r - 1.0;
2444  if (s < smin) { s = smin; }
2445  dl = dl + dl/s;
2446  l = l - dl;
2447 
2448  // Solve QP with fixed Lagrange multiplier
2449  r = solve(l,xt,x,nclip);
2450  print_iteration(nclip, r, l);
2451  }
2452 
2453  llow = l; rlow = r;
2454  }
2455 
2456  // *** Stop bracketing phase of SLBQP ***
2457 
2458 
2459  // *** Start secant phase of SLBQP ***
2461  {
2462  mfem::out << "SLBQP secant phase" << '\n';
2463  }
2464 
2465  s = 1.0 - rlow/rupp; dl = dl/s; l = lupp - dl;
2466 
2467  // Solve QP with fixed Lagrange multiplier
2468  r = solve(l,xt,x,nclip);
2469  print_iteration(nclip, r, l);
2470 
2471  while ( (fabs(r) > tol) && (nclip < max_iter) )
2472  {
2473  if (r > 0)
2474  {
2475  if (s <= 2.0)
2476  {
2477  lupp = l; rupp = r; s = 1.0 - rlow/rupp;
2478  dl = (lupp - llow)/s; l = lupp - dl;
2479  }
2480  else
2481  {
2482  s = rupp/r - 1.0;
2483  if (s < smin) { s = smin; }
2484  dl = (lupp - l)/s;
2485  lnew = 0.75*llow + 0.25*l;
2486  if (lnew < l-dl) { lnew = l-dl; }
2487  lupp = l; rupp = r; l = lnew;
2488  s = (lupp - llow)/(lupp - l);
2489  }
2490 
2491  }
2492  else
2493  {
2494  if (s >= 2.0)
2495  {
2496  llow = l; rlow = r; s = 1.0 - rlow/rupp;
2497  dl = (lupp - llow)/s; l = lupp - dl;
2498  }
2499  else
2500  {
2501  s = rlow/r - 1.0;
2502  if (s < smin) { s = smin; }
2503  dl = (l - llow)/s;
2504  lnew = 0.75*lupp + 0.25*l;
2505  if (lnew < l+dl) { lnew = l+dl; }
2506  llow = l; rlow = r; l = lnew;
2507  s = (lupp - llow)/(lupp - l);
2508  }
2509  }
2510 
2511  // Solve QP with fixed Lagrange multiplier
2512  r = solve(l,xt,x,nclip);
2513  print_iteration(nclip, r, l);
2514  }
2515 
2516  // *** Stop secant phase of SLBQP ***
2517  converged = (fabs(r) <= tol);
2518 
2519 slbqp_done:
2520 
2521  final_iter = nclip;
2522  final_norm = r;
2523 
2526  {
2527  mfem::out << "SLBQP: Number of iterations: " << final_iter << '\n'
2528  << " lambda = " << l << '\n'
2529  << " ||r|| = " << final_norm << '\n';
2530  }
2532  {
2533  mfem::out << "SLBQP: No convergence!" << '\n';
2534  }
2535 }
2536 
2537 struct WeightMinHeap
2538 {
2539  const std::vector<double> &w;
2540  std::vector<size_t> c;
2541  std::vector<int> loc;
2542 
2543  WeightMinHeap(const std::vector<double> &w_) : w(w_)
2544  {
2545  c.reserve(w.size());
2546  loc.resize(w.size());
2547  for (size_t i=0; i<w.size(); ++i) { push(i); }
2548  }
2549 
2550  size_t percolate_up(size_t pos, double val)
2551  {
2552  for (; pos > 0 && w[c[(pos-1)/2]] > val; pos = (pos-1)/2)
2553  {
2554  c[pos] = c[(pos-1)/2];
2555  loc[c[(pos-1)/2]] = pos;
2556  }
2557  return pos;
2558  }
2559 
2560  size_t percolate_down(size_t pos, double val)
2561  {
2562  while (2*pos+1 < c.size())
2563  {
2564  size_t left = 2*pos+1;
2565  size_t right = left+1;
2566  size_t tgt;
2567  if (right < c.size() && w[c[right]] < w[c[left]]) { tgt = right; }
2568  else { tgt = left; }
2569  if (w[c[tgt]] < val)
2570  {
2571  c[pos] = c[tgt];
2572  loc[c[tgt]] = pos;
2573  pos = tgt;
2574  }
2575  else
2576  {
2577  break;
2578  }
2579  }
2580  return pos;
2581  }
2582 
2583  void push(size_t i)
2584  {
2585  double val = w[i];
2586  c.push_back(0);
2587  size_t pos = c.size()-1;
2588  pos = percolate_up(pos, val);
2589  c[pos] = i;
2590  loc[i] = pos;
2591  }
2592 
2593  int pop()
2594  {
2595  size_t i = c[0];
2596  size_t j = c.back();
2597  c.pop_back();
2598  // Mark as removed
2599  loc[i] = -1;
2600  if (c.empty()) { return i; }
2601  double val = w[j];
2602  size_t pos = 0;
2603  pos = percolate_down(pos, val);
2604  c[pos] = j;
2605  loc[j] = pos;
2606  return i;
2607  }
2608 
2609  void update(size_t i)
2610  {
2611  size_t pos = loc[i];
2612  double val = w[i];
2613  pos = percolate_up(pos, val);
2614  pos = percolate_down(pos, val);
2615  c[pos] = i;
2616  loc[i] = pos;
2617  }
2618 
2619  bool picked(size_t i)
2620  {
2621  return loc[i] < 0;
2622  }
2623 };
2624 
2626 {
2627  int n = C.Width();
2628  // Scale rows by reciprocal of diagonal and take absolute value
2629  Vector D;
2630  C.GetDiag(D);
2631  int *I = C.GetI();
2632  int *J = C.GetJ();
2633  double *V = C.GetData();
2634  for (int i=0; i<n; ++i)
2635  {
2636  for (int j=I[i]; j<I[i+1]; ++j)
2637  {
2638  V[j] = abs(V[j]/D[i]);
2639  }
2640  }
2641 
2642  std::vector<double> w(n, 0.0);
2643  for (int k=0; k<n; ++k)
2644  {
2645  // Find all neighbors i of k
2646  for (int ii=I[k]; ii<I[k+1]; ++ii)
2647  {
2648  int i = J[ii];
2649  // Find value of (i,k)
2650  double C_ik = 0.0;
2651  for (int kk=I[i]; kk<I[i+1]; ++kk)
2652  {
2653  if (J[kk] == k)
2654  {
2655  C_ik = V[kk];
2656  break;
2657  }
2658  }
2659  for (int jj=I[k]; jj<I[k+1]; ++jj)
2660  {
2661  int j = J[jj];
2662  if (j == k) { continue; }
2663  double C_kj = V[jj];
2664  bool ij_exists = false;
2665  for (int jj2=I[i]; jj2<I[i+1]; ++jj2)
2666  {
2667  if (J[jj2] == j)
2668  {
2669  ij_exists = true;
2670  break;
2671  }
2672  }
2673  if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2674  }
2675  }
2676  w[k] = sqrt(w[k]);
2677  }
2678 
2679  WeightMinHeap w_heap(w);
2680 
2681  // Compute ordering
2682  p.SetSize(n);
2683  for (int ii=0; ii<n; ++ii)
2684  {
2685  int pi = w_heap.pop();
2686  p[ii] = pi;
2687  w[pi] = -1;
2688  for (int kk=I[pi]; kk<I[pi+1]; ++kk)
2689  {
2690  int k = J[kk];
2691  if (w_heap.picked(k)) { continue; }
2692  // Recompute weight
2693  w[k] = 0.0;
2694  // Find all neighbors i of k
2695  for (int ii2=I[k]; ii2<I[k+1]; ++ii2)
2696  {
2697  int i = J[ii2];
2698  if (w_heap.picked(i)) { continue; }
2699  // Find value of (i,k)
2700  double C_ik = 0.0;
2701  for (int kk2=I[i]; kk2<I[i+1]; ++kk2)
2702  {
2703  if (J[kk2] == k)
2704  {
2705  C_ik = V[kk2];
2706  break;
2707  }
2708  }
2709  for (int jj=I[k]; jj<I[k+1]; ++jj)
2710  {
2711  int j = J[jj];
2712  if (j == k || w_heap.picked(j)) { continue; }
2713  double C_kj = V[jj];
2714  bool ij_exists = false;
2715  for (int jj2=I[i]; jj2<I[i+1]; ++jj2)
2716  {
2717  if (J[jj2] == j)
2718  {
2719  ij_exists = true;
2720  break;
2721  }
2722  }
2723  if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2724  }
2725  }
2726  w[k] = sqrt(w[k]);
2727  w_heap.update(k);
2728  }
2729  }
2730 }
2731 
2732 BlockILU::BlockILU(int block_size_,
2733  Reordering reordering_,
2734  int k_fill_)
2735  : Solver(0),
2736  block_size(block_size_),
2737  k_fill(k_fill_),
2738  reordering(reordering_)
2739 { }
2740 
2742  int block_size_,
2743  Reordering reordering_,
2744  int k_fill_)
2745  : BlockILU(block_size_, reordering_, k_fill_)
2746 {
2747  SetOperator(op);
2748 }
2749 
2751 {
2752  const SparseMatrix *A = NULL;
2753 #ifdef MFEM_USE_MPI
2754  const HypreParMatrix *A_par = dynamic_cast<const HypreParMatrix *>(&op);
2755  SparseMatrix A_par_diag;
2756  if (A_par != NULL)
2757  {
2758  A_par->GetDiag(A_par_diag);
2759  A = &A_par_diag;
2760  }
2761 #endif
2762  if (A == NULL)
2763  {
2764  A = dynamic_cast<const SparseMatrix *>(&op);
2765  if (A == NULL)
2766  {
2767  MFEM_ABORT("BlockILU must be created with a SparseMatrix or HypreParMatrix");
2768  }
2769  }
2770  height = op.Height();
2771  width = op.Width();
2772  MFEM_ASSERT(A->Finalized(), "Matrix must be finalized.");
2773  CreateBlockPattern(*A);
2774  Factorize();
2775 }
2776 
2777 void BlockILU::CreateBlockPattern(const SparseMatrix &A)
2778 {
2779  MFEM_VERIFY(k_fill == 0, "Only block ILU(0) is currently supported.");
2780  if (A.Height() % block_size != 0)
2781  {
2782  MFEM_ABORT("BlockILU: block size must evenly divide the matrix size");
2783  }
2784 
2785  int nrows = A.Height();
2786  const int *I = A.GetI();
2787  const int *J = A.GetJ();
2788  const double *V = A.GetData();
2789  int nnz = 0;
2790  int nblockrows = nrows / block_size;
2791 
2792  std::vector<std::set<int>> unique_block_cols(nblockrows);
2793 
2794  for (int iblock = 0; iblock < nblockrows; ++iblock)
2795  {
2796  for (int bi = 0; bi < block_size; ++bi)
2797  {
2798  int i = iblock * block_size + bi;
2799  for (int k = I[i]; k < I[i + 1]; ++k)
2800  {
2801  unique_block_cols[iblock].insert(J[k] / block_size);
2802  }
2803  }
2804  nnz += unique_block_cols[iblock].size();
2805  }
2806 
2807  if (reordering != Reordering::NONE)
2808  {
2809  SparseMatrix C(nblockrows, nblockrows);
2810  for (int iblock = 0; iblock < nblockrows; ++iblock)
2811  {
2812  for (int jblock : unique_block_cols[iblock])
2813  {
2814  for (int bi = 0; bi < block_size; ++bi)
2815  {
2816  int i = iblock * block_size + bi;
2817  for (int k = I[i]; k < I[i + 1]; ++k)
2818  {
2819  int j = J[k];
2820  if (j >= jblock * block_size && j < (jblock + 1) * block_size)
2821  {
2822  C.Add(iblock, jblock, V[k]*V[k]);
2823  }
2824  }
2825  }
2826  }
2827  }
2828  C.Finalize(false);
2829  double *CV = C.GetData();
2830  for (int i=0; i<C.NumNonZeroElems(); ++i)
2831  {
2832  CV[i] = sqrt(CV[i]);
2833  }
2834 
2835  switch (reordering)
2836  {
2839  break;
2840  default:
2841  MFEM_ABORT("BlockILU: unknown reordering")
2842  }
2843  }
2844  else
2845  {
2846  // No reordering: permutation is identity
2847  P.SetSize(nblockrows);
2848  for (int i=0; i<nblockrows; ++i)
2849  {
2850  P[i] = i;
2851  }
2852  }
2853 
2854  // Compute inverse permutation
2855  Pinv.SetSize(nblockrows);
2856  for (int i=0; i<nblockrows; ++i)
2857  {
2858  Pinv[P[i]] = i;
2859  }
2860 
2861  // Permute columns
2862  std::vector<std::vector<int>> unique_block_cols_perminv(nblockrows);
2863  for (int i=0; i<nblockrows; ++i)
2864  {
2865  std::vector<int> &cols = unique_block_cols_perminv[i];
2866  for (int j : unique_block_cols[P[i]])
2867  {
2868  cols.push_back(Pinv[j]);
2869  }
2870  std::sort(cols.begin(), cols.end());
2871  }
2872 
2873  ID.SetSize(nblockrows);
2874  IB.SetSize(nblockrows + 1);
2875  IB[0] = 0;
2876  JB.SetSize(nnz);
2877  AB.SetSize(block_size, block_size, nnz);
2878  DB.SetSize(block_size, block_size, nblockrows);
2879  AB = 0.0;
2880  DB = 0.0;
2881  ipiv.SetSize(block_size*nblockrows);
2882  int counter = 0;
2883 
2884  for (int iblock = 0; iblock < nblockrows; ++iblock)
2885  {
2886  int iblock_perm = P[iblock];
2887  for (int jblock : unique_block_cols_perminv[iblock])
2888  {
2889  int jblock_perm = P[jblock];
2890  if (iblock == jblock)
2891  {
2892  ID[iblock] = counter;
2893  }
2894  JB[counter] = jblock;
2895  for (int bi = 0; bi < block_size; ++bi)
2896  {
2897  int i = iblock_perm*block_size + bi;
2898  for (int k = I[i]; k < I[i + 1]; ++k)
2899  {
2900  int j = J[k];
2901  if (j >= jblock_perm*block_size && j < (jblock_perm + 1)*block_size)
2902  {
2903  int bj = j - jblock_perm*block_size;
2904  double val = V[k];
2905  AB(bi, bj, counter) = val;
2906  // Extract the diagonal
2907  if (iblock == jblock)
2908  {
2909  DB(bi, bj, iblock) = val;
2910  }
2911  }
2912  }
2913  }
2914  ++counter;
2915  }
2916  IB[iblock + 1] = counter;
2917  }
2918 }
2919 
2920 void BlockILU::Factorize()
2921 {
2922  int nblockrows = Height()/block_size;
2923 
2924  // Precompute LU factorization of diagonal blocks
2925  for (int i=0; i<nblockrows; ++i)
2926  {
2927  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2928  factorization.Factor(block_size);
2929  }
2930 
2931  // Note: we use UseExternalData to extract submatrices from the tensor AB
2932  // instead of the DenseTensor call operator, because the call operator does
2933  // not allow for two simultaneous submatrix views into the same tensor
2934  DenseMatrix A_ik, A_ij, A_kj;
2935  // Loop over block rows (starting with second block row)
2936  for (int i=1; i<nblockrows; ++i)
2937  {
2938  // Find all nonzeros to the left of the diagonal in row i
2939  for (int kk=IB[i]; kk<IB[i+1]; ++kk)
2940  {
2941  int k = JB[kk];
2942  // Make sure we're still to the left of the diagonal
2943  if (k == i) { break; }
2944  if (k > i)
2945  {
2946  MFEM_ABORT("Matrix must be sorted with nonzero diagonal");
2947  }
2948  LUFactors A_kk_inv(DB.GetData(k), &ipiv[k*block_size]);
2949  A_ik.UseExternalData(&AB(0,0,kk), block_size, block_size);
2950  // A_ik = A_ik * A_kk^{-1}
2951  A_kk_inv.RightSolve(block_size, block_size, A_ik.GetData());
2952  // Modify everything to the right of k in row i
2953  for (int jj=kk+1; jj<IB[i+1]; ++jj)
2954  {
2955  int j = JB[jj];
2956  if (j <= k) { continue; } // Superfluous because JB is sorted?
2957  A_ij.UseExternalData(&AB(0,0,jj), block_size, block_size);
2958  for (int ll=IB[k]; ll<IB[k+1]; ++ll)
2959  {
2960  int l = JB[ll];
2961  if (l == j)
2962  {
2963  A_kj.UseExternalData(&AB(0,0,ll), block_size, block_size);
2964  // A_ij = A_ij - A_ik*A_kj;
2965  AddMult_a(-1.0, A_ik, A_kj, A_ij);
2966  // If we need to, update diagonal factorization
2967  if (j == i)
2968  {
2969  DB(i) = A_ij;
2970  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2971  factorization.Factor(block_size);
2972  }
2973  break;
2974  }
2975  }
2976  }
2977  }
2978  }
2979 }
2980 
2981 void BlockILU::Mult(const Vector &b, Vector &x) const
2982 {
2983  MFEM_ASSERT(height > 0, "BlockILU(0) preconditioner is not constructed");
2984  int nblockrows = Height()/block_size;
2985  y.SetSize(Height());
2986 
2987  DenseMatrix B;
2988  Vector yi, yj, xi, xj;
2989  Vector tmp(block_size);
2990  // Forward substitute to solve Ly = b
2991  // Implicitly, L has identity on the diagonal
2992  y = 0.0;
2993  for (int i=0; i<nblockrows; ++i)
2994  {
2995  yi.SetDataAndSize(&y[i*block_size], block_size);
2996  for (int ib=0; ib<block_size; ++ib)
2997  {
2998  yi[ib] = b[ib + P[i]*block_size];
2999  }
3000  for (int k=IB[i]; k<ID[i]; ++k)
3001  {
3002  int j = JB[k];
3003  const DenseMatrix &L_ij = AB(k);
3004  yj.SetDataAndSize(&y[j*block_size], block_size);
3005  // y_i = y_i - L_ij*y_j
3006  L_ij.AddMult_a(-1.0, yj, yi);
3007  }
3008  }
3009  // Backward substitution to solve Ux = y
3010  for (int i=nblockrows-1; i >= 0; --i)
3011  {
3012  xi.SetDataAndSize(&x[P[i]*block_size], block_size);
3013  for (int ib=0; ib<block_size; ++ib)
3014  {
3015  xi[ib] = y[ib + i*block_size];
3016  }
3017  for (int k=ID[i]+1; k<IB[i+1]; ++k)
3018  {
3019  int j = JB[k];
3020  const DenseMatrix &U_ij = AB(k);
3021  xj.SetDataAndSize(&x[P[j]*block_size], block_size);
3022  // x_i = x_i - U_ij*x_j
3023  U_ij.AddMult_a(-1.0, xj, xi);
3024  }
3025  LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
3026  // x_i = D_ii^{-1} x_i
3027  A_ii_inv.Solve(block_size, 1, xi);
3028  }
3029 }
3030 
3031 
3033  int it, double norm, const Vector &r, bool final)
3034 {
3035  if (!ess_dofs_list) { return; }
3036 
3037  double bc_norm_squared = 0.0;
3038  r.HostRead();
3040  for (int i = 0; i < ess_dofs_list->Size(); i++)
3041  {
3042  const double r_entry = r((*ess_dofs_list)[i]);
3043  bc_norm_squared += r_entry*r_entry;
3044  }
3045  bool print = true;
3046 #ifdef MFEM_USE_MPI
3047  MPI_Comm comm = iter_solver->GetComm();
3048  if (comm != MPI_COMM_NULL)
3049  {
3050  double glob_bc_norm_squared = 0.0;
3051  MPI_Reduce(&bc_norm_squared, &glob_bc_norm_squared, 1, MPI_DOUBLE,
3052  MPI_SUM, 0, comm);
3053  bc_norm_squared = glob_bc_norm_squared;
3054  int rank;
3055  MPI_Comm_rank(comm, &rank);
3056  print = (rank == 0);
3057  }
3058 #endif
3059  if ((it == 0 || final || bc_norm_squared > 0.0) && print)
3060  {
3061  mfem::out << " ResidualBCMonitor : b.c. residual norm = "
3062  << sqrt(bc_norm_squared) << endl;
3063  }
3064 }
3065 
3066 
3067 #ifdef MFEM_USE_SUITESPARSE
3068 
3070 {
3071  mat = NULL;
3072  Numeric = NULL;
3073  AI = AJ = NULL;
3074  if (!use_long_ints)
3075  {
3076  umfpack_di_defaults(Control);
3077  }
3078  else
3079  {
3080  umfpack_dl_defaults(Control);
3081  }
3082 }
3083 
3085 {
3086  void *Symbolic;
3087 
3088  if (Numeric)
3089  {
3090  if (!use_long_ints)
3091  {
3092  umfpack_di_free_numeric(&Numeric);
3093  }
3094  else
3095  {
3096  umfpack_dl_free_numeric(&Numeric);
3097  }
3098  }
3099 
3100  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
3101  MFEM_VERIFY(mat, "not a SparseMatrix");
3102 
3103  // UMFPack requires that the column-indices in mat corresponding to each
3104  // row be sorted.
3105  // Generally, this will modify the ordering of the entries of mat.
3107 
3108  height = mat->Height();
3109  width = mat->Width();
3110  MFEM_VERIFY(width == height, "not a square matrix");
3111 
3112  const int * Ap = mat->HostReadI();
3113  const int * Ai = mat->HostReadJ();
3114  const double * Ax = mat->HostReadData();
3115 
3116  if (!use_long_ints)
3117  {
3118  int status = umfpack_di_symbolic(width, width, Ap, Ai, Ax, &Symbolic,
3119  Control, Info);
3120  if (status < 0)
3121  {
3122  umfpack_di_report_info(Control, Info);
3123  umfpack_di_report_status(Control, status);
3124  mfem_error("UMFPackSolver::SetOperator :"
3125  " umfpack_di_symbolic() failed!");
3126  }
3127 
3128  status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric,
3129  Control, Info);
3130  if (status < 0)
3131  {
3132  umfpack_di_report_info(Control, Info);
3133  umfpack_di_report_status(Control, status);
3134  mfem_error("UMFPackSolver::SetOperator :"
3135  " umfpack_di_numeric() failed!");
3136  }
3137  umfpack_di_free_symbolic(&Symbolic);
3138  }
3139  else
3140  {
3141  SuiteSparse_long status;
3142 
3143  delete [] AJ;
3144  delete [] AI;
3145  AI = new SuiteSparse_long[width + 1];
3146  AJ = new SuiteSparse_long[Ap[width]];
3147  for (int i = 0; i <= width; i++)
3148  {
3149  AI[i] = (SuiteSparse_long)(Ap[i]);
3150  }
3151  for (int i = 0; i < Ap[width]; i++)
3152  {
3153  AJ[i] = (SuiteSparse_long)(Ai[i]);
3154  }
3155 
3156  status = umfpack_dl_symbolic(width, width, AI, AJ, Ax, &Symbolic,
3157  Control, Info);
3158  if (status < 0)
3159  {
3160  umfpack_dl_report_info(Control, Info);
3161  umfpack_dl_report_status(Control, status);
3162  mfem_error("UMFPackSolver::SetOperator :"
3163  " umfpack_dl_symbolic() failed!");
3164  }
3165 
3166  status = umfpack_dl_numeric(AI, AJ, Ax, Symbolic, &Numeric,
3167  Control, Info);
3168  if (status < 0)
3169  {
3170  umfpack_dl_report_info(Control, Info);
3171  umfpack_dl_report_status(Control, status);
3172  mfem_error("UMFPackSolver::SetOperator :"
3173  " umfpack_dl_numeric() failed!");
3174  }
3175  umfpack_dl_free_symbolic(&Symbolic);
3176  }
3177 }
3178 
3179 void UMFPackSolver::Mult(const Vector &b, Vector &x) const
3180 {
3181  if (mat == NULL)
3182  mfem_error("UMFPackSolver::Mult : matrix is not set!"
3183  " Call SetOperator first!");
3184  b.HostRead();
3185  x.HostReadWrite();
3186  if (!use_long_ints)
3187  {
3188  int status =
3189  umfpack_di_solve(UMFPACK_At, mat->HostReadI(), mat->HostReadJ(),
3190  mat->HostReadData(), x, b, Numeric, Control, Info);
3191  umfpack_di_report_info(Control, Info);
3192  if (status < 0)
3193  {
3194  umfpack_di_report_status(Control, status);
3195  mfem_error("UMFPackSolver::Mult : umfpack_di_solve() failed!");
3196  }
3197  }
3198  else
3199  {
3200  SuiteSparse_long status =
3201  umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->HostReadData(), x, b,
3202  Numeric, Control, Info);
3203  umfpack_dl_report_info(Control, Info);
3204  if (status < 0)
3205  {
3206  umfpack_dl_report_status(Control, status);
3207  mfem_error("UMFPackSolver::Mult : umfpack_dl_solve() failed!");
3208  }
3209  }
3210 }
3211 
3213 {
3214  if (mat == NULL)
3215  mfem_error("UMFPackSolver::MultTranspose : matrix is not set!"
3216  " Call SetOperator first!");
3217  b.HostRead();
3218  x.HostReadWrite();
3219  if (!use_long_ints)
3220  {
3221  int status =
3222  umfpack_di_solve(UMFPACK_A, mat->HostReadI(), mat->HostReadJ(),
3223  mat->HostReadData(), x, b, Numeric, Control, Info);
3224  umfpack_di_report_info(Control, Info);
3225  if (status < 0)
3226  {
3227  umfpack_di_report_status(Control, status);
3228  mfem_error("UMFPackSolver::MultTranspose :"
3229  " umfpack_di_solve() failed!");
3230  }
3231  }
3232  else
3233  {
3234  SuiteSparse_long status =
3235  umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->HostReadData(), x, b,
3236  Numeric, Control, Info);
3237  umfpack_dl_report_info(Control, Info);
3238  if (status < 0)
3239  {
3240  umfpack_dl_report_status(Control, status);
3241  mfem_error("UMFPackSolver::MultTranspose :"
3242  " umfpack_dl_solve() failed!");
3243  }
3244  }
3245 }
3246 
3248 {
3249  delete [] AJ;
3250  delete [] AI;
3251  if (Numeric)
3252  {
3253  if (!use_long_ints)
3254  {
3255  umfpack_di_free_numeric(&Numeric);
3256  }
3257  else
3258  {
3259  umfpack_dl_free_numeric(&Numeric);
3260  }
3261  }
3262 }
3263 
3265 {
3266  klu_defaults(&Common);
3267 }
3268 
3270 {
3271  if (Numeric)
3272  {
3273  MFEM_ASSERT(Symbolic != 0,
3274  "Had Numeric pointer in KLU, but not Symbolic");
3275  klu_free_symbolic(&Symbolic, &Common);
3276  Symbolic = 0;
3277  klu_free_numeric(&Numeric, &Common);
3278  Numeric = 0;
3279  }
3280 
3281  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
3282  MFEM_VERIFY(mat != NULL, "not a SparseMatrix");
3283 
3284  // KLU requires that the column-indices in mat corresponding to each row be
3285  // sorted. Generally, this will modify the ordering of the entries of mat.
3287 
3288  height = mat->Height();
3289  width = mat->Width();
3290  MFEM_VERIFY(width == height, "not a square matrix");
3291 
3292  int * Ap = mat->GetI();
3293  int * Ai = mat->GetJ();
3294  double * Ax = mat->GetData();
3295 
3296  Symbolic = klu_analyze( height, Ap, Ai, &Common);
3297  Numeric = klu_factor(Ap, Ai, Ax, Symbolic, &Common);
3298 }
3299 
3300 void KLUSolver::Mult(const Vector &b, Vector &x) const
3301 {
3302  MFEM_VERIFY(mat != NULL,
3303  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3304 
3305  int n = mat->Height();
3306  int numRhs = 1;
3307  // Copy B into X, so we can pass it in and overwrite it.
3308  x = b;
3309  // Solve the transpose, since KLU thinks the matrix is compressed column
3310  // format.
3311  klu_tsolve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
3312 }
3313 
3314 void KLUSolver::MultTranspose(const Vector &b, Vector &x) const
3315 {
3316  MFEM_VERIFY(mat != NULL,
3317  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3318 
3319  int n = mat->Height();
3320  int numRhs = 1;
3321  // Copy B into X, so we can pass it in and overwrite it.
3322  x = b;
3323  // Solve the regular matrix, not the transpose, since KLU thinks the matrix
3324  // is compressed column format.
3325  klu_solve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
3326 }
3327 
3329 {
3330  klu_free_symbolic (&Symbolic, &Common) ;
3331  klu_free_numeric (&Numeric, &Common) ;
3332  Symbolic = 0;
3333  Numeric = 0;
3334 }
3335 
3336 #endif // MFEM_USE_SUITESPARSE
3337 
3339  const SparseMatrix &block_dof_)
3340  : Solver(A.NumRows()), block_dof(const_cast<SparseMatrix&>(block_dof_)),
3341  block_solvers(new DenseMatrixInverse[block_dof.NumRows()])
3342 {
3343  DenseMatrix sub_A;
3344  for (int i = 0; i < block_dof.NumRows(); ++i)
3345  {
3346  local_dofs.MakeRef(block_dof.GetRowColumns(i), block_dof.RowSize(i));
3347  sub_A.SetSize(local_dofs.Size());
3348  A.GetSubMatrix(local_dofs, local_dofs, sub_A);
3349  block_solvers[i].SetOperator(sub_A);
3350  }
3351 }
3352 
3353 void DirectSubBlockSolver::Mult(const Vector &x, Vector &y) const
3354 {
3355  y.SetSize(x.Size());
3356  y = 0.0;
3357 
3358  for (int i = 0; i < block_dof.NumRows(); ++i)
3359  {
3360  local_dofs.MakeRef(block_dof.GetRowColumns(i), block_dof.RowSize(i));
3361  x.GetSubVector(local_dofs, sub_rhs);
3362  sub_sol.SetSize(local_dofs.Size());
3363  block_solvers[i].Mult(sub_rhs, sub_sol);
3364  y.AddElementVector(local_dofs, sub_sol);
3365  }
3366 }
3367 
3368 void ProductSolver::Mult(const Vector & x, Vector & y) const
3369 {
3370  y.SetSize(x.Size());
3371  y = 0.0;
3372  S0->Mult(x, y);
3373 
3374  Vector z(x.Size());
3375  z = 0.0;
3376  A->Mult(y, z);
3377  add(-1.0, z, 1.0, x, z); // z = (I - A * S0) x
3378 
3379  Vector S1z(x.Size());
3380  S1z = 0.0;
3381  S1->Mult(z, S1z);
3382  y += S1z;
3383 }
3384 
3385 void ProductSolver::MultTranspose(const Vector & x, Vector & y) const
3386 {
3387  y.SetSize(x.Size());
3388  y = 0.0;
3389  S1->MultTranspose(x, y);
3390 
3391  Vector z(x.Size());
3392  z = 0.0;
3393  A->MultTranspose(y, z);
3394  add(-1.0, z, 1.0, x, z); // z = (I - A^T * S1^T) x
3395 
3396  Vector S0Tz(x.Size());
3397  S0Tz = 0.0;
3398  S0->MultTranspose(z, S0Tz);
3399  y += S0Tz;
3400 }
3401 
3402 #ifdef MFEM_USE_MPI
3404  HypreParMatrix *aux_map,
3405  bool op_is_symmetric,
3406  bool own_aux_map)
3407  : Solver(op.NumRows()), aux_map_(aux_map, own_aux_map)
3408 {
3409  aux_system_.Reset(RAP(&op, aux_map));
3410  aux_system_.As<HypreParMatrix>()->EliminateZeroRows();
3411  aux_smoother_.Reset(new HypreSmoother(*aux_system_.As<HypreParMatrix>()));
3412  aux_smoother_.As<HypreSmoother>()->SetOperatorSymmetry(op_is_symmetric);
3413 }
3414 
3415 void AuxSpaceSmoother::Mult(const Vector &x, Vector &y, bool transpose) const
3416 {
3417  Vector aux_rhs(aux_map_->NumCols());
3418  aux_map_->MultTranspose(x, aux_rhs);
3419 
3420  Vector aux_sol(aux_rhs.Size());
3421  if (transpose)
3422  {
3423  aux_smoother_->MultTranspose(aux_rhs, aux_sol);
3424  }
3425  else
3426  {
3427  aux_smoother_->Mult(aux_rhs, aux_sol);
3428  }
3429 
3430  y.SetSize(aux_map_->NumRows());
3431  aux_map_->Mult(aux_sol, y);
3432 }
3433 #endif // MFEM_USE_MPI
3434 
3435 }
const Array< int > * ess_dofs_list
Not owned.
Definition: solvers.hpp:1011
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:943
PowerMethod helper class to estimate the largest eigenvalue of an operator using the iterative power ...
Definition: operator.hpp:953
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:336
Conjugate gradient method.
Definition: solvers.hpp:465
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.cpp:3212
double Info[UMFPACK_INFO]
Definition: solvers.hpp:1037
void Mult(const Vector &b, Vector &x) const
Solve the system LUx = b, where L and U are the block ILU factors.
Definition: solvers.cpp:2981
Chebyshev accelerated smoothing with given vector, no matrix necessary.
Definition: solvers.hpp:351
SparseMatrix * mat
Definition: solvers.hpp:1029
int GetNumIterations() const
Definition: solvers.hpp:246
const Operator * oper
Definition: solvers.hpp:120
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:264
Array< Vector * > ykArray
Definition: solvers.hpp:702
SuiteSparse_long * AI
Definition: solvers.hpp:1031
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:2673
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3353
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
void Monitor(int it, double norm, const Vector &r, const Vector &x, bool final=false) const
Definition: solvers.cpp:190
virtual void ProcessNewState(const Vector &x) const
This method can be overloaded in derived classes to perform computations that need knowledge of the n...
Definition: solvers.hpp:676
void MinimumDiscardedFillOrdering(SparseMatrix &C, Array< int > &p)
Definition: solvers.cpp:2625
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:383
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
Definition: operator.hpp:98
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition: array.hpp:308
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Definition: solvers.hpp:91
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
bool first_and_last
Information about the first and last iteration will be printed to mfem::out.
Definition: solvers.hpp:94
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:792
OperatorJacobiSmoother(const double damping=1.0)
Default constructor: the diagonal will be computed by subsequent calls to SetOperator() using the Ope...
Definition: solvers.cpp:200
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:214
class IterativeSolver * iter_solver
The last IterativeSolver to which this monitor was attached.
Definition: solvers.hpp:40
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:534
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:300
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1159
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
int * GetI()
Return the array I.
Definition: sparsemat.hpp:209
bool errors
If a fatal problem has been detected the failure will be reported to mfem::err.
Definition: solvers.hpp:82
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a tdof Vector.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.cpp:3314
void SetOperator(const Operator &op)
Definition: solvers.cpp:2750
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:655
MINRES method.
Definition: solvers.hpp:575
const Operator * GetC() const
Definition: solvers.hpp:808
bool GetConverged() const
Definition: solvers.hpp:247
Reordering
The reordering method used by the BlockILU factorization.
Definition: solvers.hpp:923
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
int print_level
(DEPRECATED) Legacy print level definition, which is left for compatibility with custom iterative sol...
Definition: solvers.hpp:130
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:449
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2748
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:300
double Norm(const Vector &x) const
Definition: solvers.hpp:171
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
const Vector * d_hi
Definition: solvers.hpp:790
double GetFinalNorm() const
Definition: solvers.hpp:248
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:219
int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, double &tol, double atol, int printit)
BiCGSTAB method. (tolerances are squared)
Definition: solvers.cpp:1565
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:525
bool IsFinite(const double &val)
Definition: vector.hpp:490
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
const Vector * x_lo
Definition: solvers.hpp:790
Operator * grad
Definition: solvers.hpp:617
const Vector * x_hi
Definition: solvers.hpp:790
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1916
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:971
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1602
General product operator: x -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:743
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:588
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:859
int GetNumConstraints() const
Definition: solvers.cpp:2325
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:86
Data type sparse matrix.
Definition: sparsemat.hpp:46
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it, int max_it, double rtol, double atol)
MINRES method without preconditioner. (tolerances are squared)
Definition: solvers.cpp:1765
virtual ~UMFPackSolver()
Definition: solvers.cpp:3247
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:446
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:529
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
const double * HostReadData() const
Definition: sparsemat.hpp:264
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
double delta
Definition: lissajous.cpp:43
const int * HostReadJ() const
Definition: sparsemat.hpp:248
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d, const Array< int > &ess_tdof_list, int order, double max_eig_estimate)
Definition: solvers.cpp:328
void SetLinearConstraint(const Vector &w_, double a_)
Definition: solvers.cpp:2352
DirectSubBlockSolver(const SparseMatrix &A, const SparseMatrix &block_dof)
Definition: solvers.cpp:3338
SuiteSparse_long * AJ
Definition: solvers.hpp:1031
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:510
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:889
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:2358
double * GetData(int k)
Definition: densemat.hpp:886
void Setup(const Vector &diag)
Definition: solvers.cpp:277
Parallel smoothers in hypre.
Definition: hypre.hpp:896
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:566
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:838
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1480
void Solve(int m, int n, double *X) const
Definition: densemat.cpp:2999
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:3269
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1588
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:55
Array< Vector * > skArray
Definition: solvers.hpp:702
prob_type prob
Definition: ex25.cpp:153
klu_symbolic * Symbolic
Definition: solvers.hpp:1070
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition: vector.cpp:626
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:442
BlockILU(int block_size_, Reordering reordering_=Reordering::MINIMUM_DISCARDED_FILL, int k_fill_=0)
Definition: solvers.cpp:2732
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1804
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:521
void AdaptiveLinRtolPostSolve(const Vector &x, const Vector &b, const int it, const double fnorm) const
Method for the adaptive linear solver rtol invoked after the linear solve.
Definition: solvers.cpp:1972
Stationary linear iteration: x &lt;- x + B (b - A x)
Definition: solvers.hpp:433
double omega
Definition: ex25.cpp:141
void SetInequalityConstraint(const Vector &dl, const Vector &dh)
Definition: solvers.cpp:2307
void SetEqualityConstraint(const Vector &c)
Definition: solvers.cpp:2299
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
void SetAbsTol(double atol)
Definition: solvers.hpp:199
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.cpp:3385
PrintLevel print_options
Output behavior for the iterative solver.
Definition: solvers.hpp:137
void SetRelTol(double rtol)
Definition: solvers.hpp:198
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1036
Abstract base class for iterative solver.
Definition: solvers.hpp:66
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1992
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:638
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1373
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3300
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:447
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
void SetAdaptiveLinRtol(const int type=2, const double rtol0=0.5, const double rtol_max=0.9, const double alpha=0.5 *(1.0+sqrt(5.0)), const double gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
Definition: solvers.cpp:1912
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
int max_iter
Limit for the number of iterations the solver is allowed to do.
Definition: solvers.hpp:151
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
GMRES method.
Definition: solvers.hpp:497
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
Definition: solvers.cpp:950
static int GuessLegacyPrintLevel(PrintLevel)
Use some heuristics to guess a legacy print level corresponding to the given PrintLevel.
Definition: solvers.cpp:149
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:267
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.cpp:2333
double a
Definition: lissajous.cpp:41
void UpdateVectors()
Definition: solvers.cpp:519
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:426
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:454
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:252
klu_common Common
Definition: solvers.hpp:1091
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:574
double rel_tol
Relative tolerance.
Definition: solvers.hpp:154
virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
This method can be overloaded in derived classes to implement line search algorithms.
Definition: solvers.hpp:671
OptimizationProblem(int insize, const Operator *C_, const Operator *D_)
In parallel, insize is the number of the local true dofs.
Definition: solvers.cpp:2289
void SetSolutionBounds(const Vector &xl, const Vector &xh)
Definition: solvers.cpp:2317
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition: operator.hpp:107
AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map, bool op_is_symmetric=true, bool own_aux_map=false)
Definition: solvers.cpp:3403
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:479
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
void SetOperator(const Operator &op)
Recompute the diagonal using the method AssembleDiagonal of the given new Operator, op.
Definition: solvers.cpp:241
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:921
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:55
void AdaptiveLinRtolPreSolve(const Vector &x, const int it, const double fnorm) const
Method for the adaptive linear solver rtol invoked before the linear solve.
Definition: solvers.cpp:1925
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:813
PrintLevel FromLegacyPrintLevel(int)
Convert a legacy print level integer to a PrintLevel object.
Definition: solvers.cpp:114
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:558
virtual void Mult(const Vector &xt, Vector &x) const
Definition: solvers.cpp:2367
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:864
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
const double alpha
Definition: ex15.cpp:369
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:247
Vector data type.
Definition: vector.hpp:60
const int * HostReadI() const
Definition: sparsemat.hpp:232
void MonitorResidual(int it, double norm, const Vector &r, bool final) override
Monitor the residual vector r.
Definition: solvers.cpp:3032
const Vector * d_lo
Definition: solvers.hpp:790
virtual void MonitorSolution(int it, double norm, const Vector &x, bool final)
Monitor the solution vector x.
Definition: solvers.hpp:54
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
int aGMRES(const Operator &A, Vector &x, const Vector &b, const Operator &M, int &max_iter, int m_max, int m_min, int m_step, double cf, double &tol, double &atol, int printit)
Definition: solvers.cpp:2134
virtual ~KLUSolver()
Definition: solvers.cpp:3328
RefCoord s[3]
void UpdateVectors()
Definition: solvers.cpp:702
BiCGSTAB method.
Definition: solvers.hpp:544
bool warnings
If a non-fatal problem has been detected some context-specific information will be reported to mfem::...
Definition: solvers.hpp:85
bool iterations
Detailed information about each iteration will be reported to mfem::out.
Definition: solvers.hpp:88
SparseMatrix * mat
Definition: solvers.hpp:1069
Base class for solvers.
Definition: operator.hpp:651
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
const OptimizationProblem * problem
Definition: solvers.hpp:823
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
double EstimateLargestEigenvalue(Operator &opr, Vector &v0, int numSteps=10, double tolerance=1e-8, int seed=12345)
Returns an estimate of the largest eigenvalue of the operator opr using the iterative power method...
Definition: operator.cpp:632
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition: solvers.cpp:2346
void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
Definition: densemat.cpp:2009
const Operator * D
Definition: solvers.hpp:789
const Vector * c_e
Definition: solvers.hpp:790
IterativeSolverMonitor * monitor
Definition: solvers.hpp:122
klu_numeric * Numeric
Definition: solvers.hpp:1071
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3368
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:438
virtual void MonitorResidual(int it, double norm, const Vector &r, bool final)
Monitor the residual vector r.
Definition: solvers.hpp:48
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:458
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3179
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:447
Settings for the output behavior of the IterativeSolver.
Definition: solvers.hpp:78
const Operator * C
Not owned, some can remain unused (NULL).
Definition: solvers.hpp:789
double abs_tol
Absolute tolerance.
Definition: solvers.hpp:157
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
void SLI(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Stationary linear iteration. (tolerances are squared)
Definition: solvers.cpp:670
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, double &tol, double atol, int printit)
GMRES method. (tolerances are squared)
Definition: solvers.cpp:1335
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1793
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3084