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