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