MFEM  v4.5.1
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 == 0.0 && 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 == 0.0 && !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 << std::endl;
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  v0.UseDevice(true);
1602  v1.UseDevice(true);
1603  w0.UseDevice(true);
1604  w1.UseDevice(true);
1605  q.UseDevice(true);
1606  u1.UseDevice(true);
1607 }
1608 
1609 void MINRESSolver::Mult(const Vector &b, Vector &x) const
1610 {
1611  // Based on the MINRES algorithm on p. 86, Fig. 6.9 in
1612  // "Iterative Krylov Methods for Large Linear Systems",
1613  // by Henk A. van der Vorst, 2003.
1614  // Extended to support an SPD preconditioner.
1615 
1616  b.UseDevice(true);
1617  x.UseDevice(true);
1618 
1619  int it;
1620  double beta, eta, gamma0, gamma1, sigma0, sigma1;
1621  double alpha, delta, rho1, rho2, rho3, norm_goal;
1622  Vector *z = (prec) ? &u1 : &v1;
1623 
1624  converged = true;
1625 
1626  if (!iterative_mode)
1627  {
1628  v1 = b;
1629  x = 0.;
1630  }
1631  else
1632  {
1633  oper->Mult(x, v1);
1634  subtract(b, v1, v1);
1635  }
1636 
1637  if (prec)
1638  {
1639  prec->Mult(v1, u1);
1640  }
1641  eta = beta = sqrt(Dot(*z, v1));
1642  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1643  gamma0 = gamma1 = 1.;
1644  sigma0 = sigma1 = 0.;
1645 
1646  norm_goal = std::max(rel_tol*eta, abs_tol);
1647 
1648  if (eta <= norm_goal)
1649  {
1650  it = 0;
1651  goto loop_end;
1652  }
1653 
1655  {
1656  mfem::out << "MINRES: iteration " << setw(3) << 0 << ": ||r||_B = "
1657  << eta << (print_options.first_and_last ? " ..." : "") << '\n';
1658  }
1659  Monitor(0, eta, *z, x);
1660 
1661  for (it = 1; it <= max_iter; it++)
1662  {
1663  v1 /= beta;
1664  if (prec)
1665  {
1666  u1 /= beta;
1667  }
1668  oper->Mult(*z, q);
1669  alpha = Dot(*z, q);
1670  MFEM_ASSERT(IsFinite(alpha), "alpha = " << alpha);
1671  if (it > 1) // (v0 == 0) for (it == 1)
1672  {
1673  q.Add(-beta, v0);
1674  }
1675  add(q, -alpha, v1, v0);
1676 
1677  delta = gamma1*alpha - gamma0*sigma1*beta;
1678  rho3 = sigma0*beta;
1679  rho2 = sigma1*alpha + gamma0*gamma1*beta;
1680  if (!prec)
1681  {
1682  beta = Norm(v0);
1683  }
1684  else
1685  {
1686  prec->Mult(v0, q);
1687  beta = sqrt(Dot(v0, q));
1688  }
1689  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1690  rho1 = hypot(delta, beta);
1691 
1692  if (it == 1)
1693  {
1694  w0.Set(1./rho1, *z); // (w0 == 0) and (w1 == 0)
1695  }
1696  else if (it == 2)
1697  {
1698  add(1./rho1, *z, -rho2/rho1, w1, w0); // (w0 == 0)
1699  }
1700  else
1701  {
1702  add(-rho3/rho1, w0, -rho2/rho1, w1, w0);
1703  w0.Add(1./rho1, *z);
1704  }
1705 
1706  gamma0 = gamma1;
1707  gamma1 = delta/rho1;
1708 
1709  x.Add(gamma1*eta, w0);
1710 
1711  sigma0 = sigma1;
1712  sigma1 = beta/rho1;
1713 
1714  eta = -sigma1*eta;
1715  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1716 
1717  if (fabs(eta) <= norm_goal)
1718  {
1719  goto loop_end;
1720  }
1721 
1723  {
1724  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1725  << fabs(eta) << '\n';
1726  }
1727  Monitor(it, fabs(eta), *z, x);
1728 
1729  if (prec)
1730  {
1731  Swap(u1, q);
1732  }
1733  Swap(v0, v1);
1734  Swap(w0, w1);
1735  }
1736  converged = false;
1737  it--;
1738 
1739 loop_end:
1740  final_iter = it;
1741  final_norm = fabs(eta);
1742 
1744  {
1745  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1746  << fabs(eta) << '\n';
1747  }
1748 
1750  {
1751  mfem::out << "MINRES: Number of iterations: " << setw(3) << final_iter << '\n';
1752  }
1753 
1754  Monitor(final_iter, final_norm, *z, x, true);
1755 
1756  // if (print_options.iteration_details || (!converged && print_options.errors))
1757  // {
1758  // oper->Mult(x, v1);
1759  // subtract(b, v1, v1);
1760  // if (prec)
1761  // {
1762  // prec->Mult(v1, u1);
1763  // }
1764  // eta = sqrt(Dot(*z, v1));
1765  // mfem::out << "MINRES: iteration " << setw(3) << it << '\n'
1766  // << " ||r||_B = " << eta << " (re-computed)" << '\n';
1767  // }
1768 
1769  if (!converged && (print_options.warnings))
1770  {
1771  mfem::out << "MINRES: No convergence!\n";
1772  }
1773 }
1774 
1775 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it,
1776  int max_it, double rtol, double atol)
1777 {
1778  MFEM_PERF_FUNCTION;
1779 
1780  MINRESSolver minres;
1781  minres.SetPrintLevel(print_it);
1782  minres.SetMaxIter(max_it);
1783  minres.SetRelTol(sqrt(rtol));
1784  minres.SetAbsTol(sqrt(atol));
1785  minres.SetOperator(A);
1786  minres.Mult(b, x);
1787 }
1788 
1789 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
1790  int print_it, int max_it, double rtol, double atol)
1791 {
1792  MINRESSolver minres;
1793  minres.SetPrintLevel(print_it);
1794  minres.SetMaxIter(max_it);
1795  minres.SetRelTol(sqrt(rtol));
1796  minres.SetAbsTol(sqrt(atol));
1797  minres.SetOperator(A);
1798  minres.SetPreconditioner(B);
1799  minres.Mult(b, x);
1800 }
1801 
1802 
1804 {
1805  oper = &op;
1806  height = op.Height();
1807  width = op.Width();
1808  MFEM_ASSERT(height == width, "square Operator is required.");
1809 
1810  r.SetSize(width);
1811  c.SetSize(width);
1812 }
1813 
1814 void NewtonSolver::Mult(const Vector &b, Vector &x) const
1815 {
1816  MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
1817  MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver).");
1818 
1819  int it;
1820  double norm0, norm, norm_goal;
1821  const bool have_b = (b.Size() == Height());
1822 
1823  if (!iterative_mode)
1824  {
1825  x = 0.0;
1826  }
1827 
1828  ProcessNewState(x);
1829 
1830  oper->Mult(x, r);
1831  if (have_b)
1832  {
1833  r -= b;
1834  }
1835 
1836  norm0 = norm = Norm(r);
1838  {
1839  mfem::out << "Newton iteration " << setw(2) << 0
1840  << " : ||r|| = " << norm << "...\n";
1841  }
1842  norm_goal = std::max(rel_tol*norm, abs_tol);
1843 
1844  prec->iterative_mode = false;
1845 
1846  // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b]
1847  for (it = 0; true; it++)
1848  {
1849  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
1851  {
1852  mfem::out << "Newton iteration " << setw(2) << it
1853  << " : ||r|| = " << norm;
1854  if (it > 0)
1855  {
1856  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
1857  }
1858  mfem::out << '\n';
1859  }
1860  Monitor(it, norm, r, x);
1861 
1862  if (norm <= norm_goal)
1863  {
1864  converged = true;
1865  break;
1866  }
1867 
1868  if (it >= max_iter)
1869  {
1870  converged = false;
1871  break;
1872  }
1873 
1874  grad = &oper->GetGradient(x);
1875  prec->SetOperator(*grad);
1876 
1877  if (lin_rtol_type)
1878  {
1879  AdaptiveLinRtolPreSolve(x, it, norm);
1880  }
1881 
1882  prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b]
1883 
1884  if (lin_rtol_type)
1885  {
1886  AdaptiveLinRtolPostSolve(c, r, it, norm);
1887  }
1888 
1889  const double c_scale = ComputeScalingFactor(x, b);
1890  if (c_scale == 0.0)
1891  {
1892  converged = false;
1893  break;
1894  }
1895  add(x, -c_scale, c, x);
1896 
1897  ProcessNewState(x);
1898 
1899  oper->Mult(x, r);
1900  if (have_b)
1901  {
1902  r -= b;
1903  }
1904  norm = Norm(r);
1905  }
1906 
1907  final_iter = it;
1908  final_norm = norm;
1909 
1912  {
1913  mfem::out << "Newton: Number of iterations: " << final_iter << '\n'
1914  << " ||r|| = " << final_norm << '\n';
1915  }
1917  {
1918  mfem::out << "Newton: No convergence!\n";
1919  }
1920 }
1921 
1923  const double rtol0,
1924  const double rtol_max,
1925  const double alpha_,
1926  const double gamma_)
1927 {
1928  lin_rtol_type = type;
1929  lin_rtol0 = rtol0;
1930  lin_rtol_max = rtol_max;
1931  this->alpha = alpha_;
1932  this->gamma = gamma_;
1933 }
1934 
1936  const int it,
1937  const double fnorm) const
1938 {
1939  // Assume that when adaptive linear solver relative tolerance is activated,
1940  // we are working with an iterative solver.
1941  auto iterative_solver = static_cast<IterativeSolver *>(prec);
1942  // Adaptive linear solver relative tolerance
1943  double eta;
1944  // Safeguard threshold
1945  double sg_threshold = 0.1;
1946 
1947  if (it == 0)
1948  {
1949  eta = lin_rtol0;
1950  }
1951  else
1952  {
1953  if (lin_rtol_type == 1)
1954  {
1955  // eta = gamma * abs(||F(x1)|| - ||F(x0) + DF(x0) s0||) / ||F(x0)||
1956  eta = gamma * abs(fnorm - lnorm_last) / fnorm_last;
1957  }
1958  else if (lin_rtol_type == 2)
1959  {
1960  // eta = gamma * (||F(x1)|| / ||F(x0)||)^alpha
1961  eta = gamma * pow(fnorm / fnorm_last, alpha);
1962  }
1963  else
1964  {
1965  MFEM_ABORT("Unknown adaptive linear solver rtol version");
1966  }
1967 
1968  // Safeguard rtol from "oversolving" ?!
1969  const double sg_eta = gamma * pow(eta_last, alpha);
1970  if (sg_eta > sg_threshold) { eta = std::max(eta, sg_eta); }
1971  }
1972 
1973  eta = std::min(eta, lin_rtol_max);
1974  iterative_solver->SetRelTol(eta);
1975  eta_last = eta;
1977  {
1978  mfem::out << "Eisenstat-Walker rtol = " << eta << "\n";
1979  }
1980 }
1981 
1983  const Vector &b,
1984  const int it,
1985  const double fnorm) const
1986 {
1987  fnorm_last = fnorm;
1988 
1989  // If version 1 is chosen, the true linear residual norm has to be computed
1990  // and in most cases we can only retrieve the preconditioned linear residual
1991  // norm.
1992  if (lin_rtol_type == 1)
1993  {
1994  // lnorm_last = ||F(x0) + DF(x0) s0||
1995  Vector linres(x.Size());
1996  grad->Mult(x, linres);
1997  linres -= b;
1998  lnorm_last = Norm(linres);
1999  }
2000 }
2001 
2002 void LBFGSSolver::Mult(const Vector &b, Vector &x) const
2003 {
2004  MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");
2005 
2006  // Quadrature points that are checked for negative Jacobians etc.
2007  Vector sk, rk, yk, rho, alpha;
2008 
2009  // r - r_{k+1}, c - descent direction
2010  sk.SetSize(width); // x_{k+1}-x_k
2011  rk.SetSize(width); // nabla(f(x_{k}))
2012  yk.SetSize(width); // r_{k+1}-r_{k}
2013  rho.SetSize(m); // 1/(dot(yk,sk)
2014  alpha.SetSize(m); // rhok*sk'*c
2015  int last_saved_id = -1;
2016 
2017  int it;
2018  double norm0, norm, norm_goal;
2019  const bool have_b = (b.Size() == Height());
2020 
2021  if (!iterative_mode)
2022  {
2023  x = 0.0;
2024  }
2025 
2026  ProcessNewState(x);
2027 
2028  // r = F(x)-b
2029  oper->Mult(x, r);
2030  if (have_b) { r -= b; }
2031 
2032  c = r; // initial descent direction
2033 
2034  norm0 = norm = Norm(r);
2036  {
2037  mfem::out << "LBFGS iteration " << setw(2) << 0
2038  << " : ||r|| = " << norm << "...\n";
2039  }
2040  norm_goal = std::max(rel_tol*norm, abs_tol);
2041  for (it = 0; true; it++)
2042  {
2043  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
2045  {
2046  mfem::out << "LBFGS iteration " << it
2047  << " : ||r|| = " << norm;
2048  if (it > 0)
2049  {
2050  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
2051  }
2052  mfem::out << '\n';
2053  }
2054 
2055  if (norm <= norm_goal)
2056  {
2057  converged = true;
2058  break;
2059  }
2060 
2061  if (it >= max_iter)
2062  {
2063  converged = false;
2064  break;
2065  }
2066 
2067  rk = r;
2068  const double c_scale = ComputeScalingFactor(x, b);
2069  if (c_scale == 0.0)
2070  {
2071  converged = false;
2072  break;
2073  }
2074  add(x, -c_scale, c, x); // x_{k+1} = x_k - c_scale*c
2075 
2076  ProcessNewState(x);
2077 
2078  oper->Mult(x, r);
2079  if (have_b)
2080  {
2081  r -= b;
2082  }
2083 
2084  // LBFGS - construct descent direction
2085  subtract(r, rk, yk); // yk = r_{k+1} - r_{k}
2086  sk = c; sk *= -c_scale; //sk = x_{k+1} - x_{k} = -c_scale*c
2087  const double gamma = Dot(sk, yk)/Dot(yk, yk);
2088 
2089  // Save last m vectors
2090  last_saved_id = (last_saved_id == m-1) ? 0 : last_saved_id+1;
2091  *skArray[last_saved_id] = sk;
2092  *ykArray[last_saved_id] = yk;
2093 
2094  c = r;
2095  for (int i = last_saved_id; i > -1; i--)
2096  {
2097  rho(i) = 1.0/Dot((*skArray[i]),(*ykArray[i]));
2098  alpha(i) = rho(i)*Dot((*skArray[i]),c);
2099  add(c, -alpha(i), (*ykArray[i]), c);
2100  }
2101  if (it > m-1)
2102  {
2103  for (int i = m-1; i > last_saved_id; i--)
2104  {
2105  rho(i) = 1./Dot((*skArray[i]), (*ykArray[i]));
2106  alpha(i) = rho(i)*Dot((*skArray[i]),c);
2107  add(c, -alpha(i), (*ykArray[i]), c);
2108  }
2109  }
2110 
2111  c *= gamma; // scale search direction
2112  if (it > m-1)
2113  {
2114  for (int i = last_saved_id+1; i < m ; i++)
2115  {
2116  double betai = rho(i)*Dot((*ykArray[i]), c);
2117  add(c, alpha(i)-betai, (*skArray[i]), c);
2118  }
2119  }
2120  for (int i = 0; i < last_saved_id+1 ; i++)
2121  {
2122  double betai = rho(i)*Dot((*ykArray[i]), c);
2123  add(c, alpha(i)-betai, (*skArray[i]), c);
2124  }
2125 
2126  norm = Norm(r);
2127  }
2128 
2129  final_iter = it;
2130  final_norm = norm;
2131 
2134  {
2135  mfem::out << "LBFGS: Number of iterations: " << final_iter << '\n'
2136  << " ||r|| = " << final_norm << '\n';
2137  }
2139  {
2140  mfem::out << "LBFGS: No convergence!\n";
2141  }
2142 }
2143 
2144 int aGMRES(const Operator &A, Vector &x, const Vector &b,
2145  const Operator &M, int &max_iter,
2146  int m_max, int m_min, int m_step, double cf,
2147  double &tol, double &atol, int printit)
2148 {
2149  int n = A.Width();
2150 
2151  int m = m_max;
2152 
2153  DenseMatrix H(m+1,m);
2154  Vector s(m+1), cs(m+1), sn(m+1);
2155  Vector w(n), av(n);
2156 
2157  double r1, resid;
2158  int i, j, k;
2159 
2160  M.Mult(b,w);
2161  double normb = w.Norml2(); // normb = ||M b||
2162  if (normb == 0.0)
2163  {
2164  normb = 1;
2165  }
2166 
2167  Vector r(n);
2168  A.Mult(x, r);
2169  subtract(b,r,w);
2170  M.Mult(w, r); // r = M (b - A x)
2171  double beta = r.Norml2(); // beta = ||r||
2172 
2173  resid = beta / normb;
2174 
2175  if (resid * resid <= tol)
2176  {
2177  tol = resid * resid;
2178  max_iter = 0;
2179  return 0;
2180  }
2181 
2182  if (printit)
2183  {
2184  mfem::out << " Pass : " << setw(2) << 1
2185  << " Iteration : " << setw(3) << 0
2186  << " (r, r) = " << beta*beta << '\n';
2187  }
2188 
2189  tol *= (normb*normb);
2190  tol = (atol > tol) ? atol : tol;
2191 
2192  m = m_max;
2193  Array<Vector *> v(m+1);
2194  for (i= 0; i<=m; i++)
2195  {
2196  v[i] = new Vector(n);
2197  (*v[i]) = 0.0;
2198  }
2199 
2200  j = 1;
2201  while (j <= max_iter)
2202  {
2203  (*v[0]) = 0.0;
2204  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
2205  s = 0.0; s(0) = beta;
2206 
2207  r1 = beta;
2208 
2209  for (i = 0; i < m && j <= max_iter; i++)
2210  {
2211  A.Mult((*v[i]),av);
2212  M.Mult(av,w); // w = M A v[i]
2213 
2214  for (k = 0; k <= i; k++)
2215  {
2216  H(k,i) = w * (*v[k]); // H(k,i) = w * v[k]
2217  w.Add(-H(k,i), (*v[k])); // w -= H(k,i) * v[k]
2218  }
2219 
2220  H(i+1,i) = w.Norml2(); // H(i+1,i) = ||w||
2221  (*v[i+1]) = 0.0;
2222  v[i+1] -> Add (1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
2223 
2224  for (k = 0; k < i; k++)
2225  {
2226  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
2227  }
2228 
2229  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
2230  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
2231  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
2232 
2233  resid = fabs(s(i+1));
2234  if (printit)
2235  {
2236  mfem::out << " Pass : " << setw(2) << j
2237  << " Iteration : " << setw(3) << i+1
2238  << " (r, r) = " << resid*resid << '\n';
2239  }
2240 
2241  if ( resid*resid < tol)
2242  {
2243  Update(x, i, H, s, v);
2244  tol = resid * resid;
2245  max_iter = j;
2246  for (i= 0; i<=m; i++)
2247  {
2248  delete v[i];
2249  }
2250  return 0;
2251  }
2252  }
2253 
2254  if (printit)
2255  {
2256  mfem::out << "Restarting..." << '\n';
2257  }
2258 
2259  Update(x, i-1, H, s, v);
2260 
2261  A.Mult(x, r);
2262  subtract(b,r,w);
2263  M.Mult(w, r); // r = M (b - A x)
2264  beta = r.Norml2(); // beta = ||r||
2265  if ( resid*resid < tol)
2266  {
2267  tol = resid * resid;
2268  max_iter = j;
2269  for (i= 0; i<=m; i++)
2270  {
2271  delete v[i];
2272  }
2273  return 0;
2274  }
2275 
2276  if (beta/r1 > cf)
2277  {
2278  if (m - m_step >= m_min)
2279  {
2280  m -= m_step;
2281  }
2282  else
2283  {
2284  m = m_max;
2285  }
2286  }
2287 
2288  j++;
2289  }
2290 
2291  tol = resid * resid;
2292  for (i= 0; i<=m; i++)
2293  {
2294  delete v[i];
2295  }
2296  return 1;
2297 }
2298 
2300  const Operator *C_,
2301  const Operator *D_)
2302  : C(C_), D(D_), c_e(NULL), d_lo(NULL), d_hi(NULL), x_lo(NULL), x_hi(NULL),
2303  input_size(insize)
2304 {
2305  if (C) { MFEM_ASSERT(C->Width() == input_size, "Wrong width of C."); }
2306  if (D) { MFEM_ASSERT(D->Width() == input_size, "Wrong width of D."); }
2307 }
2308 
2310 {
2311  MFEM_ASSERT(C, "The C operator is unspecified -- can't set constraints.");
2312  MFEM_ASSERT(c.Size() == C->Height(), "Wrong size of the constraint.");
2313 
2314  c_e = &c;
2315 }
2316 
2318  const Vector &dh)
2319 {
2320  MFEM_ASSERT(D, "The D operator is unspecified -- can't set constraints.");
2321  MFEM_ASSERT(dl.Size() == D->Height() && dh.Size() == D->Height(),
2322  "Wrong size of the constraint.");
2323 
2324  d_lo = &dl; d_hi = &dh;
2325 }
2326 
2328 {
2329  MFEM_ASSERT(xl.Size() == input_size && xh.Size() == input_size,
2330  "Wrong size of the constraint.");
2331 
2332  x_lo = &xl; x_hi = &xh;
2333 }
2334 
2336 {
2337  int m = 0;
2338  if (C) { m += C->Height(); }
2339  if (D) { m += D->Height(); }
2340  return m;
2341 }
2342 
2344 {
2345  if (print_options.warnings)
2346  {
2347  MFEM_WARNING("Objective functional is ignored as SLBQP always minimizes"
2348  "the l2 norm of (x - x_target).");
2349  }
2350  MFEM_ASSERT(prob.GetC(), "Linear constraint is not set.");
2351  MFEM_ASSERT(prob.GetC()->Height() == 1, "Solver expects scalar constraint.");
2352 
2353  problem = &prob;
2354 }
2355 
2356 void SLBQPOptimizer::SetBounds(const Vector &lo_, const Vector &hi_)
2357 {
2358  lo.SetDataAndSize(lo_.GetData(), lo_.Size());
2359  hi.SetDataAndSize(hi_.GetData(), hi_.Size());
2360 }
2361 
2363 {
2364  w.SetDataAndSize(w_.GetData(), w_.Size());
2365  a = a_;
2366 }
2367 
2368 inline void SLBQPOptimizer::print_iteration(int it, double r, double l) const
2369 {
2371  {
2372  mfem::out << "SLBQP iteration " << it << ": residual = " << r
2373  << ", lambda = " << l << '\n';
2374  }
2375 }
2376 
2377 void SLBQPOptimizer::Mult(const Vector& xt, Vector& x) const
2378 {
2379  // Based on code provided by Denis Ridzal, dridzal@sandia.gov.
2380  // Algorithm adapted from Dai and Fletcher, "New Algorithms for
2381  // Singly Linearly Constrained Quadratic Programs Subject to Lower
2382  // and Upper Bounds", Numerical Analysis Report NA/216, 2003.
2383 
2384  // Set some algorithm-specific constants and temporaries.
2385  int nclip = 0;
2386  double l = 0;
2387  double llow = 0;
2388  double lupp = 0;
2389  double lnew = 0;
2390  double dl = 2;
2391  double r = 0;
2392  double rlow = 0;
2393  double rupp = 0;
2394  double s = 0;
2395 
2396  const double smin = 0.1;
2397 
2398  const double tol = max(abs_tol, rel_tol*a);
2399 
2400  // *** Start bracketing phase of SLBQP ***
2402  {
2403  mfem::out << "SLBQP bracketing phase" << '\n';
2404  }
2405 
2406  // Solve QP with fixed Lagrange multiplier
2407  r = solve(l,xt,x,nclip);
2408  print_iteration(nclip, r, l);
2409 
2410 
2411  // If x=xt was already within bounds and satisfies the linear
2412  // constraint, then we already have the solution.
2413  if (fabs(r) <= tol)
2414  {
2415  converged = true;
2416  goto slbqp_done;
2417  }
2418 
2419  if (r < 0)
2420  {
2421  llow = l; rlow = r; l = l + dl;
2422 
2423  // Solve QP with fixed Lagrange multiplier
2424  r = solve(l,xt,x,nclip);
2425  print_iteration(nclip, r, l);
2426 
2427  while ((r < 0) && (nclip < max_iter))
2428  {
2429  llow = l;
2430  s = rlow/r - 1.0;
2431  if (s < smin) { s = smin; }
2432  dl = dl + dl/s;
2433  l = l + dl;
2434 
2435  // Solve QP with fixed Lagrange multiplier
2436  r = solve(l,xt,x,nclip);
2437  print_iteration(nclip, r, l);
2438  }
2439 
2440  lupp = l; rupp = r;
2441  }
2442  else
2443  {
2444  lupp = l; rupp = r; l = l - dl;
2445 
2446  // Solve QP with fixed Lagrange multiplier
2447  r = solve(l,xt,x,nclip);
2448  print_iteration(nclip, r, l);
2449 
2450  while ((r > 0) && (nclip < max_iter))
2451  {
2452  lupp = l;
2453  s = rupp/r - 1.0;
2454  if (s < smin) { s = smin; }
2455  dl = dl + dl/s;
2456  l = l - dl;
2457 
2458  // Solve QP with fixed Lagrange multiplier
2459  r = solve(l,xt,x,nclip);
2460  print_iteration(nclip, r, l);
2461  }
2462 
2463  llow = l; rlow = r;
2464  }
2465 
2466  // *** Stop bracketing phase of SLBQP ***
2467 
2468 
2469  // *** Start secant phase of SLBQP ***
2471  {
2472  mfem::out << "SLBQP secant phase" << '\n';
2473  }
2474 
2475  s = 1.0 - rlow/rupp; dl = dl/s; l = lupp - dl;
2476 
2477  // Solve QP with fixed Lagrange multiplier
2478  r = solve(l,xt,x,nclip);
2479  print_iteration(nclip, r, l);
2480 
2481  while ( (fabs(r) > tol) && (nclip < max_iter) )
2482  {
2483  if (r > 0)
2484  {
2485  if (s <= 2.0)
2486  {
2487  lupp = l; rupp = r; s = 1.0 - rlow/rupp;
2488  dl = (lupp - llow)/s; l = lupp - dl;
2489  }
2490  else
2491  {
2492  s = rupp/r - 1.0;
2493  if (s < smin) { s = smin; }
2494  dl = (lupp - l)/s;
2495  lnew = 0.75*llow + 0.25*l;
2496  if (lnew < l-dl) { lnew = l-dl; }
2497  lupp = l; rupp = r; l = lnew;
2498  s = (lupp - llow)/(lupp - l);
2499  }
2500 
2501  }
2502  else
2503  {
2504  if (s >= 2.0)
2505  {
2506  llow = l; rlow = r; s = 1.0 - rlow/rupp;
2507  dl = (lupp - llow)/s; l = lupp - dl;
2508  }
2509  else
2510  {
2511  s = rlow/r - 1.0;
2512  if (s < smin) { s = smin; }
2513  dl = (l - llow)/s;
2514  lnew = 0.75*lupp + 0.25*l;
2515  if (lnew < l+dl) { lnew = l+dl; }
2516  llow = l; rlow = r; l = lnew;
2517  s = (lupp - llow)/(lupp - l);
2518  }
2519  }
2520 
2521  // Solve QP with fixed Lagrange multiplier
2522  r = solve(l,xt,x,nclip);
2523  print_iteration(nclip, r, l);
2524  }
2525 
2526  // *** Stop secant phase of SLBQP ***
2527  converged = (fabs(r) <= tol);
2528 
2529 slbqp_done:
2530 
2531  final_iter = nclip;
2532  final_norm = r;
2533 
2536  {
2537  mfem::out << "SLBQP: Number of iterations: " << final_iter << '\n'
2538  << " lambda = " << l << '\n'
2539  << " ||r|| = " << final_norm << '\n';
2540  }
2542  {
2543  mfem::out << "SLBQP: No convergence!" << '\n';
2544  }
2545 }
2546 
2547 struct WeightMinHeap
2548 {
2549  const std::vector<double> &w;
2550  std::vector<size_t> c;
2551  std::vector<int> loc;
2552 
2553  WeightMinHeap(const std::vector<double> &w_) : w(w_)
2554  {
2555  c.reserve(w.size());
2556  loc.resize(w.size());
2557  for (size_t i=0; i<w.size(); ++i) { push(i); }
2558  }
2559 
2560  size_t percolate_up(size_t pos, double val)
2561  {
2562  for (; pos > 0 && w[c[(pos-1)/2]] > val; pos = (pos-1)/2)
2563  {
2564  c[pos] = c[(pos-1)/2];
2565  loc[c[(pos-1)/2]] = pos;
2566  }
2567  return pos;
2568  }
2569 
2570  size_t percolate_down(size_t pos, double val)
2571  {
2572  while (2*pos+1 < c.size())
2573  {
2574  size_t left = 2*pos+1;
2575  size_t right = left+1;
2576  size_t tgt;
2577  if (right < c.size() && w[c[right]] < w[c[left]]) { tgt = right; }
2578  else { tgt = left; }
2579  if (w[c[tgt]] < val)
2580  {
2581  c[pos] = c[tgt];
2582  loc[c[tgt]] = pos;
2583  pos = tgt;
2584  }
2585  else
2586  {
2587  break;
2588  }
2589  }
2590  return pos;
2591  }
2592 
2593  void push(size_t i)
2594  {
2595  double val = w[i];
2596  c.push_back(0);
2597  size_t pos = c.size()-1;
2598  pos = percolate_up(pos, val);
2599  c[pos] = i;
2600  loc[i] = pos;
2601  }
2602 
2603  int pop()
2604  {
2605  size_t i = c[0];
2606  size_t j = c.back();
2607  c.pop_back();
2608  // Mark as removed
2609  loc[i] = -1;
2610  if (c.empty()) { return i; }
2611  double val = w[j];
2612  size_t pos = 0;
2613  pos = percolate_down(pos, val);
2614  c[pos] = j;
2615  loc[j] = pos;
2616  return i;
2617  }
2618 
2619  void update(size_t i)
2620  {
2621  size_t pos = loc[i];
2622  double val = w[i];
2623  pos = percolate_up(pos, val);
2624  pos = percolate_down(pos, val);
2625  c[pos] = i;
2626  loc[i] = pos;
2627  }
2628 
2629  bool picked(size_t i)
2630  {
2631  return loc[i] < 0;
2632  }
2633 };
2634 
2636 {
2637  int n = C.Width();
2638  // Scale rows by reciprocal of diagonal and take absolute value
2639  Vector D;
2640  C.GetDiag(D);
2641  int *I = C.GetI();
2642  int *J = C.GetJ();
2643  double *V = C.GetData();
2644  for (int i=0; i<n; ++i)
2645  {
2646  for (int j=I[i]; j<I[i+1]; ++j)
2647  {
2648  V[j] = abs(V[j]/D[i]);
2649  }
2650  }
2651 
2652  std::vector<double> w(n, 0.0);
2653  for (int k=0; k<n; ++k)
2654  {
2655  // Find all neighbors i of k
2656  for (int ii=I[k]; ii<I[k+1]; ++ii)
2657  {
2658  int i = J[ii];
2659  // Find value of (i,k)
2660  double C_ik = 0.0;
2661  for (int kk=I[i]; kk<I[i+1]; ++kk)
2662  {
2663  if (J[kk] == k)
2664  {
2665  C_ik = V[kk];
2666  break;
2667  }
2668  }
2669  for (int jj=I[k]; jj<I[k+1]; ++jj)
2670  {
2671  int j = J[jj];
2672  if (j == k) { continue; }
2673  double C_kj = V[jj];
2674  bool ij_exists = false;
2675  for (int jj2=I[i]; jj2<I[i+1]; ++jj2)
2676  {
2677  if (J[jj2] == j)
2678  {
2679  ij_exists = true;
2680  break;
2681  }
2682  }
2683  if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2684  }
2685  }
2686  w[k] = sqrt(w[k]);
2687  }
2688 
2689  WeightMinHeap w_heap(w);
2690 
2691  // Compute ordering
2692  p.SetSize(n);
2693  for (int ii=0; ii<n; ++ii)
2694  {
2695  int pi = w_heap.pop();
2696  p[ii] = pi;
2697  w[pi] = -1;
2698  for (int kk=I[pi]; kk<I[pi+1]; ++kk)
2699  {
2700  int k = J[kk];
2701  if (w_heap.picked(k)) { continue; }
2702  // Recompute weight
2703  w[k] = 0.0;
2704  // Find all neighbors i of k
2705  for (int ii2=I[k]; ii2<I[k+1]; ++ii2)
2706  {
2707  int i = J[ii2];
2708  if (w_heap.picked(i)) { continue; }
2709  // Find value of (i,k)
2710  double C_ik = 0.0;
2711  for (int kk2=I[i]; kk2<I[i+1]; ++kk2)
2712  {
2713  if (J[kk2] == k)
2714  {
2715  C_ik = V[kk2];
2716  break;
2717  }
2718  }
2719  for (int jj=I[k]; jj<I[k+1]; ++jj)
2720  {
2721  int j = J[jj];
2722  if (j == k || w_heap.picked(j)) { continue; }
2723  double C_kj = V[jj];
2724  bool ij_exists = false;
2725  for (int jj2=I[i]; jj2<I[i+1]; ++jj2)
2726  {
2727  if (J[jj2] == j)
2728  {
2729  ij_exists = true;
2730  break;
2731  }
2732  }
2733  if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2734  }
2735  }
2736  w[k] = sqrt(w[k]);
2737  w_heap.update(k);
2738  }
2739  }
2740 }
2741 
2742 BlockILU::BlockILU(int block_size_,
2743  Reordering reordering_,
2744  int k_fill_)
2745  : Solver(0),
2746  block_size(block_size_),
2747  k_fill(k_fill_),
2748  reordering(reordering_)
2749 { }
2750 
2752  int block_size_,
2753  Reordering reordering_,
2754  int k_fill_)
2755  : BlockILU(block_size_, reordering_, k_fill_)
2756 {
2757  SetOperator(op);
2758 }
2759 
2761 {
2762  const SparseMatrix *A = NULL;
2763 #ifdef MFEM_USE_MPI
2764  const HypreParMatrix *A_par = dynamic_cast<const HypreParMatrix *>(&op);
2765  SparseMatrix A_par_diag;
2766  if (A_par != NULL)
2767  {
2768  A_par->GetDiag(A_par_diag);
2769  A = &A_par_diag;
2770  }
2771 #endif
2772  if (A == NULL)
2773  {
2774  A = dynamic_cast<const SparseMatrix *>(&op);
2775  if (A == NULL)
2776  {
2777  MFEM_ABORT("BlockILU must be created with a SparseMatrix or HypreParMatrix");
2778  }
2779  }
2780  height = op.Height();
2781  width = op.Width();
2782  MFEM_ASSERT(A->Finalized(), "Matrix must be finalized.");
2783  CreateBlockPattern(*A);
2784  Factorize();
2785 }
2786 
2787 void BlockILU::CreateBlockPattern(const SparseMatrix &A)
2788 {
2789  MFEM_VERIFY(k_fill == 0, "Only block ILU(0) is currently supported.");
2790  if (A.Height() % block_size != 0)
2791  {
2792  MFEM_ABORT("BlockILU: block size must evenly divide the matrix size");
2793  }
2794 
2795  int nrows = A.Height();
2796  const int *I = A.GetI();
2797  const int *J = A.GetJ();
2798  const double *V = A.GetData();
2799  int nnz = 0;
2800  int nblockrows = nrows / block_size;
2801 
2802  std::vector<std::set<int>> unique_block_cols(nblockrows);
2803 
2804  for (int iblock = 0; iblock < nblockrows; ++iblock)
2805  {
2806  for (int bi = 0; bi < block_size; ++bi)
2807  {
2808  int i = iblock * block_size + bi;
2809  for (int k = I[i]; k < I[i + 1]; ++k)
2810  {
2811  unique_block_cols[iblock].insert(J[k] / block_size);
2812  }
2813  }
2814  nnz += unique_block_cols[iblock].size();
2815  }
2816 
2817  if (reordering != Reordering::NONE)
2818  {
2819  SparseMatrix C(nblockrows, nblockrows);
2820  for (int iblock = 0; iblock < nblockrows; ++iblock)
2821  {
2822  for (int jblock : unique_block_cols[iblock])
2823  {
2824  for (int bi = 0; bi < block_size; ++bi)
2825  {
2826  int i = iblock * block_size + bi;
2827  for (int k = I[i]; k < I[i + 1]; ++k)
2828  {
2829  int j = J[k];
2830  if (j >= jblock * block_size && j < (jblock + 1) * block_size)
2831  {
2832  C.Add(iblock, jblock, V[k]*V[k]);
2833  }
2834  }
2835  }
2836  }
2837  }
2838  C.Finalize(false);
2839  double *CV = C.GetData();
2840  for (int i=0; i<C.NumNonZeroElems(); ++i)
2841  {
2842  CV[i] = sqrt(CV[i]);
2843  }
2844 
2845  switch (reordering)
2846  {
2849  break;
2850  default:
2851  MFEM_ABORT("BlockILU: unknown reordering")
2852  }
2853  }
2854  else
2855  {
2856  // No reordering: permutation is identity
2857  P.SetSize(nblockrows);
2858  for (int i=0; i<nblockrows; ++i)
2859  {
2860  P[i] = i;
2861  }
2862  }
2863 
2864  // Compute inverse permutation
2865  Pinv.SetSize(nblockrows);
2866  for (int i=0; i<nblockrows; ++i)
2867  {
2868  Pinv[P[i]] = i;
2869  }
2870 
2871  // Permute columns
2872  std::vector<std::vector<int>> unique_block_cols_perminv(nblockrows);
2873  for (int i=0; i<nblockrows; ++i)
2874  {
2875  std::vector<int> &cols = unique_block_cols_perminv[i];
2876  for (int j : unique_block_cols[P[i]])
2877  {
2878  cols.push_back(Pinv[j]);
2879  }
2880  std::sort(cols.begin(), cols.end());
2881  }
2882 
2883  ID.SetSize(nblockrows);
2884  IB.SetSize(nblockrows + 1);
2885  IB[0] = 0;
2886  JB.SetSize(nnz);
2887  AB.SetSize(block_size, block_size, nnz);
2888  DB.SetSize(block_size, block_size, nblockrows);
2889  AB = 0.0;
2890  DB = 0.0;
2891  ipiv.SetSize(block_size*nblockrows);
2892  int counter = 0;
2893 
2894  for (int iblock = 0; iblock < nblockrows; ++iblock)
2895  {
2896  int iblock_perm = P[iblock];
2897  for (int jblock : unique_block_cols_perminv[iblock])
2898  {
2899  int jblock_perm = P[jblock];
2900  if (iblock == jblock)
2901  {
2902  ID[iblock] = counter;
2903  }
2904  JB[counter] = jblock;
2905  for (int bi = 0; bi < block_size; ++bi)
2906  {
2907  int i = iblock_perm*block_size + bi;
2908  for (int k = I[i]; k < I[i + 1]; ++k)
2909  {
2910  int j = J[k];
2911  if (j >= jblock_perm*block_size && j < (jblock_perm + 1)*block_size)
2912  {
2913  int bj = j - jblock_perm*block_size;
2914  double val = V[k];
2915  AB(bi, bj, counter) = val;
2916  // Extract the diagonal
2917  if (iblock == jblock)
2918  {
2919  DB(bi, bj, iblock) = val;
2920  }
2921  }
2922  }
2923  }
2924  ++counter;
2925  }
2926  IB[iblock + 1] = counter;
2927  }
2928 }
2929 
2930 void BlockILU::Factorize()
2931 {
2932  int nblockrows = Height()/block_size;
2933 
2934  // Precompute LU factorization of diagonal blocks
2935  for (int i=0; i<nblockrows; ++i)
2936  {
2937  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2938  factorization.Factor(block_size);
2939  }
2940 
2941  // Note: we use UseExternalData to extract submatrices from the tensor AB
2942  // instead of the DenseTensor call operator, because the call operator does
2943  // not allow for two simultaneous submatrix views into the same tensor
2944  DenseMatrix A_ik, A_ij, A_kj;
2945  // Loop over block rows (starting with second block row)
2946  for (int i=1; i<nblockrows; ++i)
2947  {
2948  // Find all nonzeros to the left of the diagonal in row i
2949  for (int kk=IB[i]; kk<IB[i+1]; ++kk)
2950  {
2951  int k = JB[kk];
2952  // Make sure we're still to the left of the diagonal
2953  if (k == i) { break; }
2954  if (k > i)
2955  {
2956  MFEM_ABORT("Matrix must be sorted with nonzero diagonal");
2957  }
2958  LUFactors A_kk_inv(DB.GetData(k), &ipiv[k*block_size]);
2959  A_ik.UseExternalData(&AB(0,0,kk), block_size, block_size);
2960  // A_ik = A_ik * A_kk^{-1}
2961  A_kk_inv.RightSolve(block_size, block_size, A_ik.GetData());
2962  // Modify everything to the right of k in row i
2963  for (int jj=kk+1; jj<IB[i+1]; ++jj)
2964  {
2965  int j = JB[jj];
2966  if (j <= k) { continue; } // Superfluous because JB is sorted?
2967  A_ij.UseExternalData(&AB(0,0,jj), block_size, block_size);
2968  for (int ll=IB[k]; ll<IB[k+1]; ++ll)
2969  {
2970  int l = JB[ll];
2971  if (l == j)
2972  {
2973  A_kj.UseExternalData(&AB(0,0,ll), block_size, block_size);
2974  // A_ij = A_ij - A_ik*A_kj;
2975  AddMult_a(-1.0, A_ik, A_kj, A_ij);
2976  // If we need to, update diagonal factorization
2977  if (j == i)
2978  {
2979  DB(i) = A_ij;
2980  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2981  factorization.Factor(block_size);
2982  }
2983  break;
2984  }
2985  }
2986  }
2987  }
2988  }
2989 }
2990 
2991 void BlockILU::Mult(const Vector &b, Vector &x) const
2992 {
2993  MFEM_ASSERT(height > 0, "BlockILU(0) preconditioner is not constructed");
2994  int nblockrows = Height()/block_size;
2995  y.SetSize(Height());
2996 
2997  DenseMatrix B;
2998  Vector yi, yj, xi, xj;
2999  Vector tmp(block_size);
3000  // Forward substitute to solve Ly = b
3001  // Implicitly, L has identity on the diagonal
3002  y = 0.0;
3003  for (int i=0; i<nblockrows; ++i)
3004  {
3005  yi.SetDataAndSize(&y[i*block_size], block_size);
3006  for (int ib=0; ib<block_size; ++ib)
3007  {
3008  yi[ib] = b[ib + P[i]*block_size];
3009  }
3010  for (int k=IB[i]; k<ID[i]; ++k)
3011  {
3012  int j = JB[k];
3013  const DenseMatrix &L_ij = AB(k);
3014  yj.SetDataAndSize(&y[j*block_size], block_size);
3015  // y_i = y_i - L_ij*y_j
3016  L_ij.AddMult_a(-1.0, yj, yi);
3017  }
3018  }
3019  // Backward substitution to solve Ux = y
3020  for (int i=nblockrows-1; i >= 0; --i)
3021  {
3022  xi.SetDataAndSize(&x[P[i]*block_size], block_size);
3023  for (int ib=0; ib<block_size; ++ib)
3024  {
3025  xi[ib] = y[ib + i*block_size];
3026  }
3027  for (int k=ID[i]+1; k<IB[i+1]; ++k)
3028  {
3029  int j = JB[k];
3030  const DenseMatrix &U_ij = AB(k);
3031  xj.SetDataAndSize(&x[P[j]*block_size], block_size);
3032  // x_i = x_i - U_ij*x_j
3033  U_ij.AddMult_a(-1.0, xj, xi);
3034  }
3035  LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
3036  // x_i = D_ii^{-1} x_i
3037  A_ii_inv.Solve(block_size, 1, xi);
3038  }
3039 }
3040 
3041 
3043  int it, double norm, const Vector &r, bool final)
3044 {
3045  if (!ess_dofs_list) { return; }
3046 
3047  double bc_norm_squared = 0.0;
3048  r.HostRead();
3050  for (int i = 0; i < ess_dofs_list->Size(); i++)
3051  {
3052  const double r_entry = r((*ess_dofs_list)[i]);
3053  bc_norm_squared += r_entry*r_entry;
3054  }
3055  bool print = true;
3056 #ifdef MFEM_USE_MPI
3057  MPI_Comm comm = iter_solver->GetComm();
3058  if (comm != MPI_COMM_NULL)
3059  {
3060  double glob_bc_norm_squared = 0.0;
3061  MPI_Reduce(&bc_norm_squared, &glob_bc_norm_squared, 1, MPI_DOUBLE,
3062  MPI_SUM, 0, comm);
3063  bc_norm_squared = glob_bc_norm_squared;
3064  int rank;
3065  MPI_Comm_rank(comm, &rank);
3066  print = (rank == 0);
3067  }
3068 #endif
3069  if ((it == 0 || final || bc_norm_squared > 0.0) && print)
3070  {
3071  mfem::out << " ResidualBCMonitor : b.c. residual norm = "
3072  << sqrt(bc_norm_squared) << endl;
3073  }
3074 }
3075 
3076 
3077 #ifdef MFEM_USE_SUITESPARSE
3078 
3080 {
3081  mat = NULL;
3082  Numeric = NULL;
3083  AI = AJ = NULL;
3084  if (!use_long_ints)
3085  {
3086  umfpack_di_defaults(Control);
3087  }
3088  else
3089  {
3090  umfpack_dl_defaults(Control);
3091  }
3092 }
3093 
3095 {
3096  void *Symbolic;
3097 
3098  if (Numeric)
3099  {
3100  if (!use_long_ints)
3101  {
3102  umfpack_di_free_numeric(&Numeric);
3103  }
3104  else
3105  {
3106  umfpack_dl_free_numeric(&Numeric);
3107  }
3108  }
3109 
3110  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
3111  MFEM_VERIFY(mat, "not a SparseMatrix");
3112 
3113  // UMFPack requires that the column-indices in mat corresponding to each
3114  // row be sorted.
3115  // Generally, this will modify the ordering of the entries of mat.
3117 
3118  height = mat->Height();
3119  width = mat->Width();
3120  MFEM_VERIFY(width == height, "not a square matrix");
3121 
3122  const int * Ap = mat->HostReadI();
3123  const int * Ai = mat->HostReadJ();
3124  const double * Ax = mat->HostReadData();
3125 
3126  if (!use_long_ints)
3127  {
3128  int status = umfpack_di_symbolic(width, width, Ap, Ai, Ax, &Symbolic,
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_symbolic() failed!");
3136  }
3137 
3138  status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric,
3139  Control, Info);
3140  if (status < 0)
3141  {
3142  umfpack_di_report_info(Control, Info);
3143  umfpack_di_report_status(Control, status);
3144  mfem_error("UMFPackSolver::SetOperator :"
3145  " umfpack_di_numeric() failed!");
3146  }
3147  umfpack_di_free_symbolic(&Symbolic);
3148  }
3149  else
3150  {
3151  SuiteSparse_long status;
3152 
3153  delete [] AJ;
3154  delete [] AI;
3155  AI = new SuiteSparse_long[width + 1];
3156  AJ = new SuiteSparse_long[Ap[width]];
3157  for (int i = 0; i <= width; i++)
3158  {
3159  AI[i] = (SuiteSparse_long)(Ap[i]);
3160  }
3161  for (int i = 0; i < Ap[width]; i++)
3162  {
3163  AJ[i] = (SuiteSparse_long)(Ai[i]);
3164  }
3165 
3166  status = umfpack_dl_symbolic(width, width, AI, AJ, Ax, &Symbolic,
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_symbolic() failed!");
3174  }
3175 
3176  status = umfpack_dl_numeric(AI, AJ, Ax, Symbolic, &Numeric,
3177  Control, Info);
3178  if (status < 0)
3179  {
3180  umfpack_dl_report_info(Control, Info);
3181  umfpack_dl_report_status(Control, status);
3182  mfem_error("UMFPackSolver::SetOperator :"
3183  " umfpack_dl_numeric() failed!");
3184  }
3185  umfpack_dl_free_symbolic(&Symbolic);
3186  }
3187 }
3188 
3189 void UMFPackSolver::Mult(const Vector &b, Vector &x) const
3190 {
3191  if (mat == NULL)
3192  mfem_error("UMFPackSolver::Mult : matrix is not set!"
3193  " Call SetOperator first!");
3194  b.HostRead();
3195  x.HostReadWrite();
3196  if (!use_long_ints)
3197  {
3198  int status =
3199  umfpack_di_solve(UMFPACK_At, mat->HostReadI(), mat->HostReadJ(),
3200  mat->HostReadData(), x, b, Numeric, Control, Info);
3201  umfpack_di_report_info(Control, Info);
3202  if (status < 0)
3203  {
3204  umfpack_di_report_status(Control, status);
3205  mfem_error("UMFPackSolver::Mult : umfpack_di_solve() failed!");
3206  }
3207  }
3208  else
3209  {
3210  SuiteSparse_long status =
3211  umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->HostReadData(), x, b,
3212  Numeric, Control, Info);
3213  umfpack_dl_report_info(Control, Info);
3214  if (status < 0)
3215  {
3216  umfpack_dl_report_status(Control, status);
3217  mfem_error("UMFPackSolver::Mult : umfpack_dl_solve() failed!");
3218  }
3219  }
3220 }
3221 
3223 {
3224  if (mat == NULL)
3225  mfem_error("UMFPackSolver::MultTranspose : matrix is not set!"
3226  " Call SetOperator first!");
3227  b.HostRead();
3228  x.HostReadWrite();
3229  if (!use_long_ints)
3230  {
3231  int status =
3232  umfpack_di_solve(UMFPACK_A, mat->HostReadI(), mat->HostReadJ(),
3233  mat->HostReadData(), x, b, Numeric, Control, Info);
3234  umfpack_di_report_info(Control, Info);
3235  if (status < 0)
3236  {
3237  umfpack_di_report_status(Control, status);
3238  mfem_error("UMFPackSolver::MultTranspose :"
3239  " umfpack_di_solve() failed!");
3240  }
3241  }
3242  else
3243  {
3244  SuiteSparse_long status =
3245  umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->HostReadData(), x, b,
3246  Numeric, Control, Info);
3247  umfpack_dl_report_info(Control, Info);
3248  if (status < 0)
3249  {
3250  umfpack_dl_report_status(Control, status);
3251  mfem_error("UMFPackSolver::MultTranspose :"
3252  " umfpack_dl_solve() failed!");
3253  }
3254  }
3255 }
3256 
3258 {
3259  delete [] AJ;
3260  delete [] AI;
3261  if (Numeric)
3262  {
3263  if (!use_long_ints)
3264  {
3265  umfpack_di_free_numeric(&Numeric);
3266  }
3267  else
3268  {
3269  umfpack_dl_free_numeric(&Numeric);
3270  }
3271  }
3272 }
3273 
3275 {
3276  klu_defaults(&Common);
3277 }
3278 
3280 {
3281  if (Numeric)
3282  {
3283  MFEM_ASSERT(Symbolic != 0,
3284  "Had Numeric pointer in KLU, but not Symbolic");
3285  klu_free_symbolic(&Symbolic, &Common);
3286  Symbolic = 0;
3287  klu_free_numeric(&Numeric, &Common);
3288  Numeric = 0;
3289  }
3290 
3291  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
3292  MFEM_VERIFY(mat != NULL, "not a SparseMatrix");
3293 
3294  // KLU requires that the column-indices in mat corresponding to each row be
3295  // sorted. Generally, this will modify the ordering of the entries of mat.
3297 
3298  height = mat->Height();
3299  width = mat->Width();
3300  MFEM_VERIFY(width == height, "not a square matrix");
3301 
3302  int * Ap = mat->GetI();
3303  int * Ai = mat->GetJ();
3304  double * Ax = mat->GetData();
3305 
3306  Symbolic = klu_analyze( height, Ap, Ai, &Common);
3307  Numeric = klu_factor(Ap, Ai, Ax, Symbolic, &Common);
3308 }
3309 
3310 void KLUSolver::Mult(const Vector &b, Vector &x) const
3311 {
3312  MFEM_VERIFY(mat != NULL,
3313  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3314 
3315  int n = mat->Height();
3316  int numRhs = 1;
3317  // Copy B into X, so we can pass it in and overwrite it.
3318  x = b;
3319  // Solve the transpose, since KLU thinks the matrix is compressed column
3320  // format.
3321  klu_tsolve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
3322 }
3323 
3324 void KLUSolver::MultTranspose(const Vector &b, Vector &x) const
3325 {
3326  MFEM_VERIFY(mat != NULL,
3327  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3328 
3329  int n = mat->Height();
3330  int numRhs = 1;
3331  // Copy B into X, so we can pass it in and overwrite it.
3332  x = b;
3333  // Solve the regular matrix, not the transpose, since KLU thinks the matrix
3334  // is compressed column format.
3335  klu_solve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
3336 }
3337 
3339 {
3340  klu_free_symbolic (&Symbolic, &Common) ;
3341  klu_free_numeric (&Numeric, &Common) ;
3342  Symbolic = 0;
3343  Numeric = 0;
3344 }
3345 
3346 #endif // MFEM_USE_SUITESPARSE
3347 
3349  const SparseMatrix &block_dof_)
3350  : Solver(A.NumRows()), block_dof(const_cast<SparseMatrix&>(block_dof_)),
3351  block_solvers(new DenseMatrixInverse[block_dof.NumRows()])
3352 {
3353  DenseMatrix sub_A;
3354  for (int i = 0; i < block_dof.NumRows(); ++i)
3355  {
3356  local_dofs.MakeRef(block_dof.GetRowColumns(i), block_dof.RowSize(i));
3357  sub_A.SetSize(local_dofs.Size());
3358  A.GetSubMatrix(local_dofs, local_dofs, sub_A);
3359  block_solvers[i].SetOperator(sub_A);
3360  }
3361 }
3362 
3363 void DirectSubBlockSolver::Mult(const Vector &x, Vector &y) const
3364 {
3365  y.SetSize(x.Size());
3366  y = 0.0;
3367 
3368  for (int i = 0; i < block_dof.NumRows(); ++i)
3369  {
3370  local_dofs.MakeRef(block_dof.GetRowColumns(i), block_dof.RowSize(i));
3371  x.GetSubVector(local_dofs, sub_rhs);
3372  sub_sol.SetSize(local_dofs.Size());
3373  block_solvers[i].Mult(sub_rhs, sub_sol);
3374  y.AddElementVector(local_dofs, sub_sol);
3375  }
3376 }
3377 
3378 void ProductSolver::Mult(const Vector & x, Vector & y) const
3379 {
3380  y.SetSize(x.Size());
3381  y = 0.0;
3382  S0->Mult(x, y);
3383 
3384  Vector z(x.Size());
3385  z = 0.0;
3386  A->Mult(y, z);
3387  add(-1.0, z, 1.0, x, z); // z = (I - A * S0) x
3388 
3389  Vector S1z(x.Size());
3390  S1z = 0.0;
3391  S1->Mult(z, S1z);
3392  y += S1z;
3393 }
3394 
3395 void ProductSolver::MultTranspose(const Vector & x, Vector & y) const
3396 {
3397  y.SetSize(x.Size());
3398  y = 0.0;
3399  S1->MultTranspose(x, y);
3400 
3401  Vector z(x.Size());
3402  z = 0.0;
3403  A->MultTranspose(y, z);
3404  add(-1.0, z, 1.0, x, z); // z = (I - A^T * S1^T) x
3405 
3406  Vector S0Tz(x.Size());
3407  S0Tz = 0.0;
3408  S0->MultTranspose(z, S0Tz);
3409  y += S0Tz;
3410 }
3411 
3413  : Solver(0, false), global_size(-1)
3414 #ifdef MFEM_USE_MPI
3415  , parallel(false)
3416 #endif
3417 { }
3418 
3419 #ifdef MFEM_USE_MPI
3420 OrthoSolver::OrthoSolver(MPI_Comm mycomm_)
3421  : Solver(0, false), mycomm(mycomm_), global_size(-1), parallel(true) { }
3422 #endif
3423 
3425 {
3426  solver = &s;
3427  height = s.Height();
3428  width = s.Width();
3429  MFEM_VERIFY(height == width, "Solver must be a square Operator!");
3430  global_size = -1; // lazy evaluated
3431 }
3432 
3434 {
3435  MFEM_VERIFY(solver, "Solver hasn't been set, call SetSolver() first.");
3436  solver->SetOperator(op);
3437  height = solver->Height();
3438  width = solver->Width();
3439  MFEM_VERIFY(height == width, "Solver must be a square Operator!");
3440  global_size = -1; // lazy evaluated
3441 }
3442 
3443 void OrthoSolver::Mult(const Vector &b, Vector &x) const
3444 {
3445  MFEM_VERIFY(solver, "Solver hasn't been set, call SetSolver() first.");
3446  MFEM_VERIFY(height == solver->Height(),
3447  "solver was modified externally! call SetSolver() again!");
3448  MFEM_VERIFY(height == b.Size(), "incompatible input Vector size!");
3449  MFEM_VERIFY(height == x.Size(), "incompatible output Vector size!");
3450 
3451  // Orthogonalize input
3452  Orthogonalize(b, b_ortho);
3453 
3454  // Propagate iterative_mode to the solver:
3455  solver->iterative_mode = iterative_mode;
3456 
3457  // Apply the Solver
3458  solver->Mult(b_ortho, x);
3459 
3460  // Orthogonalize output
3461  Orthogonalize(x, x);
3462 }
3463 
3464 void OrthoSolver::Orthogonalize(const Vector &v, Vector &v_ortho) const
3465 {
3466  if (global_size == -1)
3467  {
3468  global_size = height;
3469 #ifdef MFEM_USE_MPI
3470  if (parallel)
3471  {
3472  MPI_Allreduce(MPI_IN_PLACE, &global_size, 1, HYPRE_MPI_BIG_INT,
3473  MPI_SUM, mycomm);
3474  }
3475 #endif
3476  }
3477 
3478  // TODO: GPU/device implementation
3479 
3480  double global_sum = v.Sum();
3481 
3482 #ifdef MFEM_USE_MPI
3483  if (parallel)
3484  {
3485  MPI_Allreduce(MPI_IN_PLACE, &global_sum, 1, MPI_DOUBLE, MPI_SUM, mycomm);
3486  }
3487 #endif
3488 
3489  double ratio = global_sum / static_cast<double>(global_size);
3490  v_ortho.SetSize(v.Size());
3491  v.HostRead();
3492  v_ortho.HostWrite();
3493  for (int i = 0; i < v_ortho.Size(); ++i)
3494  {
3495  v_ortho(i) = v(i) - ratio;
3496  }
3497 }
3498 
3499 #ifdef MFEM_USE_MPI
3501  HypreParMatrix *aux_map,
3502  bool op_is_symmetric,
3503  bool own_aux_map)
3504  : Solver(op.NumRows()), aux_map_(aux_map, own_aux_map)
3505 {
3506  aux_system_.Reset(RAP(&op, aux_map));
3507  aux_system_.As<HypreParMatrix>()->EliminateZeroRows();
3508  aux_smoother_.Reset(new HypreSmoother(*aux_system_.As<HypreParMatrix>()));
3509  aux_smoother_.As<HypreSmoother>()->SetOperatorSymmetry(op_is_symmetric);
3510 }
3511 
3512 void AuxSpaceSmoother::Mult(const Vector &x, Vector &y, bool transpose) const
3513 {
3514  Vector aux_rhs(aux_map_->NumCols());
3515  aux_map_->MultTranspose(x, aux_rhs);
3516 
3517  Vector aux_sol(aux_rhs.Size());
3518  if (transpose)
3519  {
3520  aux_smoother_->MultTranspose(aux_rhs, aux_sol);
3521  }
3522  else
3523  {
3524  aux_smoother_->Mult(aux_rhs, aux_sol);
3525  }
3526 
3527  y.SetSize(aux_map_->NumRows());
3528  aux_map_->Mult(aux_sol, y);
3529 }
3530 #endif // MFEM_USE_MPI
3531 
3532 }
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:959
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:344
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:3222
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:2991
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
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:3363
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
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:2635
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:391
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:99
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:73
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:812
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:227
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:547
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 double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:461
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:200
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
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:3324
void SetOperator(const Operator &op)
Definition: solvers.cpp:2760
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:659
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:209
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:457
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:118
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2882
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:313
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:94
const Vector * d_hi
Definition: solvers.hpp:790
virtual void SetOperator(const Operator &op)
Set the Operator that is the OrthoSolver is to invert (approximately).
Definition: solvers.cpp:3433
double GetFinalNorm() const
Definition: solvers.hpp:248
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
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:482
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:2282
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:1609
General product operator: x -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:747
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:2335
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:87
Data type sparse matrix.
Definition: sparsemat.hpp:50
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
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:1775
virtual ~UMFPackSolver()
Definition: solvers.cpp:3257
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:457
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:542
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:277
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
double delta
Definition: lissajous.cpp:43
const int * HostReadJ() const
Definition: sparsemat.hpp:261
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:2362
DirectSubBlockSolver(const SparseMatrix &A, const SparseMatrix &block_dof)
Definition: solvers.cpp:3348
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:2368
double * GetData(int k)
Definition: densemat.hpp:1058
void Setup(const Vector &diag)
Definition: solvers.cpp:277
Parallel smoothers in hypre.
Definition: hypre.hpp:962
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:590
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:1005
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1480
virtual void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3381
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:3279
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
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3213
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:639
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:453
BlockILU(int block_size_, Reordering reordering_=Reordering::MINIMUM_DISCARDED_FILL, int k_fill_=0)
Definition: solvers.cpp:2742
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1814
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:551
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:1982
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:2317
void SetEqualityConstraint(const Vector &c)
Definition: solvers.cpp:2309
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:3395
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:2002
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:729
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:70
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3310
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:460
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:1922
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:76
GMRES method.
Definition: solvers.hpp:497
void Mult(const Vector &b, Vector &x) const
Perform the action of the OrthoSolver: P * solver * P where P is the projection to the subspace of ve...
Definition: solvers.cpp:3443
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:157
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:2343
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:413
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:465
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:665
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:2299
void SetSolutionBounds(const Vector &xl, const Vector &xh)
Definition: solvers.cpp:2327
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition: operator.hpp:108
AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map, bool op_is_symmetric=true, bool own_aux_map=false)
Definition: solvers.cpp:3500
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:479
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:1935
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:904
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:2377
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:260
Vector data type.
Definition: vector.hpp:60
const int * HostReadI() const
Definition: sparsemat.hpp:245
void SetSolver(Solver &s)
Set the solver used by the OrthoSolver.
Definition: solvers.cpp:3424
void MonitorResidual(int it, double norm, const Vector &r, bool final) override
Monitor the residual vector r.
Definition: solvers.cpp:3042
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:2144
virtual ~KLUSolver()
Definition: solvers.cpp:3338
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:655
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:343
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition: solvers.cpp:2356
void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
Definition: densemat.cpp:2375
const Operator * D
Definition: solvers.hpp:789
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:920
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:3378
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:449
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:469
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:3189
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:1803
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3094