MFEM  v4.6.0 Finite element discretization library
solvers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "linalg.hpp"
13 #include "../general/annotation.hpp"
14 #include "../general/forall.hpp"
15 #include "../general/globals.hpp"
16 #include "../fem/bilinearform.hpp"
17 #include <iostream>
18 #include <iomanip>
19 #include <algorithm>
20 #include <cmath>
21 #include <set>
22
23 namespace mfem
24 {
25
26 using namespace std;
27
29  : Solver(0, true)
30 {
31  oper = NULL;
32  prec = NULL;
33  max_iter = 10;
34  rel_tol = abs_tol = 0.0;
35 #ifdef MFEM_USE_MPI
36  dot_prod_type = 0;
37 #endif
38 }
39
40 #ifdef MFEM_USE_MPI
41
43  : Solver(0, true)
44 {
45  oper = NULL;
46  prec = NULL;
47  max_iter = 10;
48  rel_tol = abs_tol = 0.0;
49  dot_prod_type = 1;
50  comm = comm_;
51 }
52
53 #endif // MFEM_USE_MPI
54
55 double IterativeSolver::Dot(const Vector &x, const Vector &y) const
56 {
57 #ifndef MFEM_USE_MPI
58  return (x * y);
59 #else
60  if (dot_prod_type == 0)
61  {
62  return (x * y);
63  }
64  else
65  {
66  return InnerProduct(comm, x, y);
67  }
68 #endif
69 }
70
71 void IterativeSolver::SetPrintLevel(int print_lvl)
72 {
74  int print_level_ = print_lvl;
75
76 #ifdef MFEM_USE_MPI
77  if (dot_prod_type != 0)
78  {
79  int rank;
80  MPI_Comm_rank(comm, &rank);
81  if (rank != 0) // Suppress output.
82  {
83  print_level_ = -1;
85  }
86  }
87 #endif
88
89  print_level = print_level_;
90 }
91
93 {
94  print_options = options;
95
96  int derived_print_level = GuessLegacyPrintLevel(options);
97
98 #ifdef MFEM_USE_MPI
99  if (dot_prod_type != 0)
100  {
101  int rank;
102  MPI_Comm_rank(comm, &rank);
103  if (rank != 0)
104  {
105  derived_print_level = -1;
107  }
108  }
109 #endif
110
111  print_level = derived_print_level;
112 }
113
115  int print_level_)
116 {
117 #ifdef MFEM_USE_MPI
118  int rank = 0;
119  if (comm != MPI_COMM_NULL)
120  {
121  MPI_Comm_rank(comm, &rank);
122  }
123 #endif
124
125  switch (print_level_)
126  {
127  case -1:
128  return PrintLevel();
129  case 0:
130  return PrintLevel().Errors().Warnings();
131  case 1:
132  return PrintLevel().Errors().Warnings().Iterations();
133  case 2:
134  return PrintLevel().Errors().Warnings().Summary();
135  case 3:
136  return PrintLevel().Errors().Warnings().FirstAndLast();
137  default:
138 #ifdef MFEM_USE_MPI
139  if (rank == 0)
140 #endif
141  {
142  MFEM_WARNING("Unknown print level " << print_level_ <<
143  ". Defaulting to level 0.");
144  }
145  return PrintLevel().Errors().Warnings();
146  }
147 }
148
150 {
151  if (print_options_.iterations)
152  {
153  return 1;
154  }
155  else if (print_options_.first_and_last)
156  {
157  return 3;
158  }
159  else if (print_options_.summary)
160  {
161  return 2;
162  }
163  else if (print_options_.errors && print_options_.warnings)
164  {
165  return 0;
166  }
167  else
168  {
169  return -1;
170  }
171 }
172
174 {
175  prec = &pr;
176  prec->iterative_mode = false;
177 }
178
180 {
181  oper = &op;
182  height = op.Height();
183  width = op.Width();
184  if (prec)
185  {
186  prec->SetOperator(*oper);
187  }
188 }
189
190 void IterativeSolver::Monitor(int it, double norm, const Vector& r,
191  const Vector& x, bool final) const
192 {
193  if (monitor != nullptr)
194  {
195  monitor->MonitorResidual(it, norm, r, final);
196  monitor->MonitorSolution(it, norm, x, final);
197  }
198 }
199
201  : damping(dmpng),
202  ess_tdof_list(nullptr),
203  oper(nullptr),
204  allow_updates(true)
205 { }
206
208  const Array<int> &ess_tdofs,
209  const double dmpng)
210  :
211  Solver(a.FESpace()->GetTrueVSize()),
212  dinv(height),
213  damping(dmpng),
214  ess_tdof_list(&ess_tdofs),
215  residual(height),
216  allow_updates(false)
217 {
218  Vector &diag(residual);
219  a.AssembleDiagonal(diag);
220  // 'a' cannot be used for iterative_mode == true because its size may be
221  // different.
222  oper = nullptr;
223  Setup(diag);
224 }
225
227  const Array<int> &ess_tdofs,
228  const double dmpng)
229  :
230  Solver(d.Size()),
231  dinv(height),
232  damping(dmpng),
233  ess_tdof_list(&ess_tdofs),
234  residual(height),
235  oper(NULL),
236  allow_updates(false)
237 {
238  Setup(d);
239 }
240
242 {
243  if (!allow_updates)
244  {
245  // original behavior of this method
246  oper = &op; return;
247  }
248
249  // Treat (Par)BilinearForm objects as a special case since their
250  // AssembleDiagonal method returns the true-dof diagonal whereas the form
251  // itself may act as an ldof operator. This is for compatibility with the
252  // constructor that takes a BilinearForm parameter.
253  const BilinearForm *blf = dynamic_cast<const BilinearForm *>(&op);
254  if (blf)
255  {
256  // 'a' cannot be used for iterative_mode == true because its size may be
257  // different.
258  oper = nullptr;
259  height = width = blf->FESpace()->GetTrueVSize();
260  }
261  else
262  {
263  oper = &op;
264  height = op.Height();
265  width = op.Width();
266  MFEM_ASSERT(height == width, "not a square matrix!");
267  // ess_tdof_list is only used with BilinearForm
268  ess_tdof_list = nullptr;
269  }
270  dinv.SetSize(height);
271  residual.SetSize(height);
272  Vector &diag(residual);
273  op.AssembleDiagonal(diag);
274  Setup(diag);
275 }
276
278 {
279  residual.UseDevice(true);
280  const double delta = damping;
281  auto D = diag.Read();
282  auto DI = dinv.Write();
283  const bool use_abs_diag_ = use_abs_diag;
284  mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
285  {
286  if (D[i] == 0.0)
287  {
288  MFEM_ABORT_KERNEL("Zero diagonal entry in OperatorJacobiSmoother");
289  }
290  if (!use_abs_diag_) { DI[i] = delta / D[i]; }
291  else { DI[i] = delta / std::abs(D[i]); }
292  });
293  if (ess_tdof_list && ess_tdof_list->Size() > 0)
294  {
295  auto I = ess_tdof_list->Read();
296  mfem::forall(ess_tdof_list->Size(), [=] MFEM_HOST_DEVICE (int i)
297  {
298  DI[I[i]] = delta;
299  });
300  }
301 }
302
304 {
305  // For empty MPI ranks, height may be 0:
306  // MFEM_VERIFY(Height() > 0, "The diagonal hasn't been computed.");
307  MFEM_ASSERT(x.Size() == Width(), "invalid input vector");
308  MFEM_ASSERT(y.Size() == Height(), "invalid output vector");
309
310  if (iterative_mode)
311  {
312  MFEM_VERIFY(oper, "iterative_mode == true requires the forward operator");
313  oper->Mult(y, residual); // r = A y
314  subtract(x, residual, residual); // r = x - A y
315  }
316  else
317  {
318  residual = x;
319  y.UseDevice(true);
320  y = 0.0;
321  }
322  auto DI = dinv.Read();
323  auto R = residual.Read();
324  auto Y = y.ReadWrite();
325  mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
326  {
327  Y[i] += DI[i] * R[i];
328  });
329 }
330
332  const Vector &d,
333  const Array<int>& ess_tdofs,
334  int order_, double max_eig_estimate_)
335  :
336  Solver(d.Size()),
337  order(order_),
338  max_eig_estimate(max_eig_estimate_),
339  N(d.Size()),
340  dinv(N),
341  diag(d),
342  coeffs(order),
343  ess_tdof_list(ess_tdofs),
344  residual(N),
345  oper(&oper_) { Setup(); }
346
347 #ifdef MFEM_USE_MPI
349  const Vector &d,
350  const Array<int>& ess_tdofs,
351  int order_, MPI_Comm comm, int power_iterations, double power_tolerance)
352 #else
354  const Vector &d,
355  const Array<int>& ess_tdofs,
356  int order_, int power_iterations, double power_tolerance)
357 #endif
358  : Solver(d.Size()),
359  order(order_),
360  N(d.Size()),
361  dinv(N),
362  diag(d),
363  coeffs(order),
364  ess_tdof_list(ess_tdofs),
365  residual(N),
366  oper(&oper_)
367 {
368  OperatorJacobiSmoother invDiagOperator(diag, ess_tdofs, 1.0);
369  ProductOperator diagPrecond(&invDiagOperator, oper, false, false);
370
371 #ifdef MFEM_USE_MPI
372  PowerMethod powerMethod(comm);
373 #else
374  PowerMethod powerMethod;
375 #endif
376  Vector ev(oper->Width());
377  max_eig_estimate = powerMethod.EstimateLargestEigenvalue(diagPrecond, ev,
378  power_iterations, power_tolerance);
379
380  Setup();
381 }
382
384  const Vector &d,
385  const Array<int>& ess_tdofs,
386  int order_, double max_eig_estimate_)
387  : OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, max_eig_estimate_) { }
388
389 #ifdef MFEM_USE_MPI
391  const Vector &d,
392  const Array<int>& ess_tdofs,
393  int order_, MPI_Comm comm, int power_iterations, double power_tolerance)
394  : OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, comm,
395  power_iterations, power_tolerance) { }
396 #else
398  const Vector &d,
399  const Array<int>& ess_tdofs,
400  int order_, int power_iterations, double power_tolerance)
401  : OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, power_iterations,
402  power_tolerance) { }
403 #endif
404
406 {
407  // Invert diagonal
408  residual.UseDevice(true);
409  auto D = diag.Read();
410  auto X = dinv.Write();
411  mfem::forall(N, [=] MFEM_HOST_DEVICE (int i) { X[i] = 1.0 / D[i]; });
412  auto I = ess_tdof_list.Read();
413  mfem::forall(ess_tdof_list.Size(), [=] MFEM_HOST_DEVICE (int i)
414  {
415  X[I[i]] = 1.0;
416  });
417
418  // Set up Chebyshev coefficients
419  // For reference, see e.g., Parallel multigrid smoothing: polynomial versus
420  // Gauss-Seidel by Adams et al.
421  double upper_bound = 1.2 * max_eig_estimate;
422  double lower_bound = 0.3 * max_eig_estimate;
423  double theta = 0.5 * (upper_bound + lower_bound);
424  double delta = 0.5 * (upper_bound - lower_bound);
425
426  switch (order-1)
427  {
428  case 0:
429  {
430  coeffs[0] = 1.0 / theta;
431  break;
432  }
433  case 1:
434  {
435  double tmp_0 = 1.0/(pow(delta, 2) - 2*pow(theta, 2));
436  coeffs[0] = -4*theta*tmp_0;
437  coeffs[1] = 2*tmp_0;
438  break;
439  }
440  case 2:
441  {
442  double tmp_0 = 3*pow(delta, 2);
443  double tmp_1 = pow(theta, 2);
444  double tmp_2 = 1.0/(-4*pow(theta, 3) + theta*tmp_0);
445  coeffs[0] = tmp_2*(tmp_0 - 12*tmp_1);
446  coeffs[1] = 12/(tmp_0 - 4*tmp_1);
447  coeffs[2] = -4*tmp_2;
448  break;
449  }
450  case 3:
451  {
452  double tmp_0 = pow(delta, 2);
453  double tmp_1 = pow(theta, 2);
454  double tmp_2 = 8*tmp_0;
455  double tmp_3 = 1.0/(pow(delta, 4) + 8*pow(theta, 4) - tmp_1*tmp_2);
456  coeffs[0] = tmp_3*(32*pow(theta, 3) - 16*theta*tmp_0);
457  coeffs[1] = tmp_3*(-48*tmp_1 + tmp_2);
458  coeffs[2] = 32*theta*tmp_3;
459  coeffs[3] = -8*tmp_3;
460  break;
461  }
462  case 4:
463  {
464  double tmp_0 = 5*pow(delta, 4);
465  double tmp_1 = pow(theta, 4);
466  double tmp_2 = pow(theta, 2);
467  double tmp_3 = pow(delta, 2);
468  double tmp_4 = 60*tmp_3;
469  double tmp_5 = 20*tmp_3;
470  double tmp_6 = 1.0/(16*pow(theta, 5) - pow(theta, 3)*tmp_5 + theta*tmp_0);
471  double tmp_7 = 160*tmp_2;
472  double tmp_8 = 1.0/(tmp_0 + 16*tmp_1 - tmp_2*tmp_5);
473  coeffs[0] = tmp_6*(tmp_0 + 80*tmp_1 - tmp_2*tmp_4);
474  coeffs[1] = tmp_8*(tmp_4 - tmp_7);
475  coeffs[2] = tmp_6*(-tmp_5 + tmp_7);
476  coeffs[3] = -80*tmp_8;
477  coeffs[4] = 16*tmp_6;
478  break;
479  }
480  default:
481  MFEM_ABORT("Chebyshev smoother not implemented for order = " << order);
482  }
483 }
484
486 {
487  if (iterative_mode)
488  {
489  MFEM_ABORT("Chebyshev smoother not implemented for iterative mode");
490  }
491
492  if (!oper)
493  {
494  MFEM_ABORT("Chebyshev smoother requires operator");
495  }
496
497  residual = x;
498  helperVector.SetSize(x.Size());
499
500  y.UseDevice(true);
501  y = 0.0;
502
503  for (int k = 0; k < order; ++k)
504  {
505  // Apply
506  if (k > 0)
507  {
508  oper->Mult(residual, helperVector);
509  residual = helperVector;
510  }
511
512  // Scale residual by inverse diagonal
513  const int n = N;
514  auto Dinv = dinv.Read();
515  auto R = residual.ReadWrite();
516  mfem::forall(n, [=] MFEM_HOST_DEVICE (int i) { R[i] *= Dinv[i]; });
517
518  // Add weighted contribution to y
519  auto Y = y.ReadWrite();
520  auto C = coeffs.Read();
521  mfem::forall(n, [=] MFEM_HOST_DEVICE (int i) { Y[i] += C[k] * R[i]; });
522  }
523 }
524
526 {
527  r.SetSize(width);
528  z.SetSize(width);
529 }
530
531 void SLISolver::Mult(const Vector &b, Vector &x) const
532 {
533  int i;
534
535  // Optimized preconditioned SLI with fixed number of iterations and given
536  // initial guess
537  if (rel_tol == 0.0 && iterative_mode && prec)
538  {
539  for (i = 0; i < max_iter; i++)
540  {
541  oper->Mult(x, r); // r = A x
542  subtract(b, r, r); // r = b - A x
543  prec->Mult(r, z); // z = B r
544  add(x, 1.0, z, x); // x = x + B (b - A x)
545  }
546  converged = true;
547  final_iter = i;
548  return;
549  }
550
551  // Optimized preconditioned SLI with fixed number of iterations and zero
552  // initial guess
553  if (rel_tol == 0.0 && !iterative_mode && prec)
554  {
555  prec->Mult(b, x); // x = B b (initial guess 0)
556  for (i = 1; i < max_iter; i++)
557  {
558  oper->Mult(x, r); // r = A x
559  subtract(b, r, r); // r = b - A x
560  prec->Mult(r, z); // z = B r
561  add(x, 1.0, z, x); // x = x + B (b - A x)
562  }
563  converged = true;
564  final_iter = i;
565  return;
566  }
567
568  // General version of SLI with a relative tolerance, optional preconditioner
569  // and optional initial guess
570  double r0, nom, nom0, nomold = 1, cf;
571
572  if (iterative_mode)
573  {
574  oper->Mult(x, r);
575  subtract(b, r, r); // r = b - A x
576  }
577  else
578  {
579  r = b;
580  x = 0.0;
581  }
582
583  if (prec)
584  {
585  prec->Mult(r, z); // z = B r
586  nom0 = nom = sqrt(Dot(z, z));
587  }
588  else
589  {
590  nom0 = nom = sqrt(Dot(r, r));
591  }
592  initial_norm = nom0;
593
595  {
596  mfem::out << " Iteration : " << setw(3) << right << 0 << " ||Br|| = "
597  << nom << (print_options.first_and_last ? " ..." : "") << '\n';
598  }
599
600  r0 = std::max(nom*rel_tol, abs_tol);
601  if (nom <= r0)
602  {
603  converged = true;
604  final_iter = 0;
605  final_norm = nom;
606  return;
607  }
608
609  // start iteration
610  converged = false;
612  for (i = 1; true; )
613  {
614  if (prec) // x = x + B (b - A x)
615  {
616  add(x, 1.0, z, x);
617  }
618  else
619  {
620  add(x, 1.0, r, x);
621  }
622
623  oper->Mult(x, r);
624  subtract(b, r, r); // r = b - A x
625
626  if (prec)
627  {
628  prec->Mult(r, z); // z = B r
629  nom = sqrt(Dot(z, z));
630  }
631  else
632  {
633  nom = sqrt(Dot(r, r));
634  }
635
636  cf = nom/nomold;
637  nomold = nom;
638
639  bool done = false;
640  if (nom < r0)
641  {
642  converged = true;
643  final_iter = i;
644  done = true;
645  }
646
647  if (++i > max_iter)
648  {
649  done = true;
650  }
651
653  {
654  mfem::out << " Iteration : " << setw(3) << right << (i-1)
655  << " ||Br|| = " << setw(11) << left << nom
656  << "\tConv. rate: " << cf << '\n';
657  }
658
659  if (done) { break; }
660  }
661
663  {
664  const auto rf = pow (nom/nom0, 1.0/final_iter);
665  mfem::out << "SLI: Number of iterations: " << final_iter << '\n'
666  << "Conv. rate: " << cf << '\n'
667  << "Average reduction factor: "<< rf << '\n';
668  }
670  {
671  mfem::out << "SLI: No convergence!" << '\n';
672  }
673
674  final_norm = nom;
675 }
676
677 void SLI(const Operator &A, const Vector &b, Vector &x,
678  int print_iter, int max_num_iter,
679  double RTOLERANCE, double ATOLERANCE)
680 {
681  MFEM_PERF_FUNCTION;
682
683  SLISolver sli;
684  sli.SetPrintLevel(print_iter);
685  sli.SetMaxIter(max_num_iter);
686  sli.SetRelTol(sqrt(RTOLERANCE));
687  sli.SetAbsTol(sqrt(ATOLERANCE));
688  sli.SetOperator(A);
689  sli.Mult(b, x);
690 }
691
692 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
693  int print_iter, int max_num_iter,
694  double RTOLERANCE, double ATOLERANCE)
695 {
696  MFEM_PERF_FUNCTION;
697
698  SLISolver sli;
699  sli.SetPrintLevel(print_iter);
700  sli.SetMaxIter(max_num_iter);
701  sli.SetRelTol(sqrt(RTOLERANCE));
702  sli.SetAbsTol(sqrt(ATOLERANCE));
703  sli.SetOperator(A);
704  sli.SetPreconditioner(B);
705  sli.Mult(b, x);
706 }
707
708
710 {
712
713  r.SetSize(width, mt); r.UseDevice(true);
714  d.SetSize(width, mt); d.UseDevice(true);
715  z.SetSize(width, mt); z.UseDevice(true);
716 }
717
718 void CGSolver::Mult(const Vector &b, Vector &x) const
719 {
720  int i;
721  double r0, den, nom, nom0, betanom, alpha, beta;
722
723  x.UseDevice(true);
724  if (iterative_mode)
725  {
726  oper->Mult(x, r);
727  subtract(b, r, r); // r = b - A x
728  }
729  else
730  {
731  r = b;
732  x = 0.0;
733  }
734
735  if (prec)
736  {
737  prec->Mult(r, z); // z = B r
738  d = z;
739  }
740  else
741  {
742  d = r;
743  }
744  nom0 = nom = Dot(d, r);
745  if (nom0 >= 0.0) { initial_norm = sqrt(nom0); }
746  MFEM_ASSERT(IsFinite(nom), "nom = " << nom);
748  {
749  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
750  << nom << (print_options.first_and_last ? " ...\n" : "\n");
751  }
752  Monitor(0, nom, r, x);
753
754  if (nom < 0.0)
755  {
757  {
758  mfem::out << "PCG: The preconditioner is not positive definite. (Br, r) = "
759  << nom << '\n';
760  }
761  converged = false;
762  final_iter = 0;
763  initial_norm = nom;
764  final_norm = nom;
765  return;
766  }
767  r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
768  if (nom <= r0)
769  {
770  converged = true;
771  final_iter = 0;
772  final_norm = sqrt(nom);
773  return;
774  }
775
776  oper->Mult(d, z); // z = A d
777  den = Dot(z, d);
778  MFEM_ASSERT(IsFinite(den), "den = " << den);
779  if (den <= 0.0)
780  {
781  if (Dot(d, d) > 0.0 && print_options.warnings)
782  {
783  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
784  << den << '\n';
785  }
786  if (den == 0.0)
787  {
788  converged = false;
789  final_iter = 0;
790  final_norm = sqrt(nom);
791  return;
792  }
793  }
794
795  // start iteration
796  converged = false;
798  for (i = 1; true; )
799  {
800  alpha = nom/den;
801  add(x, alpha, d, x); // x = x + alpha d
802  add(r, -alpha, z, r); // r = r - alpha A d
803
804  if (prec)
805  {
806  prec->Mult(r, z); // z = B r
807  betanom = Dot(r, z);
808  }
809  else
810  {
811  betanom = Dot(r, r);
812  }
813  MFEM_ASSERT(IsFinite(betanom), "betanom = " << betanom);
814  if (betanom < 0.0)
815  {
817  {
818  mfem::out << "PCG: The preconditioner is not positive definite. (Br, r) = "
819  << betanom << '\n';
820  }
821  converged = false;
822  final_iter = i;
823  break;
824  }
825
827  {
828  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
829  << betanom << std::endl;
830  }
831
832  Monitor(i, betanom, r, x);
833
834  if (betanom <= r0)
835  {
836  converged = true;
837  final_iter = i;
838  break;
839  }
840
841  if (++i > max_iter)
842  {
843  break;
844  }
845
846  beta = betanom/nom;
847  if (prec)
848  {
849  add(z, beta, d, d); // d = z + beta d
850  }
851  else
852  {
853  add(r, beta, d, d);
854  }
855  oper->Mult(d, z); // z = A d
856  den = Dot(d, z);
857  MFEM_ASSERT(IsFinite(den), "den = " << den);
858  if (den <= 0.0)
859  {
860  if (Dot(d, d) > 0.0 && print_options.warnings)
861  {
862  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
863  << den << '\n';
864  }
865  if (den == 0.0)
866  {
867  final_iter = i;
868  break;
869  }
870  }
871  nom = betanom;
872  }
874  {
875  mfem::out << " Iteration : " << setw(3) << final_iter << " (B r, r) = "
876  << betanom << '\n';
877  }
879  {
880  mfem::out << "PCG: Number of iterations: " << final_iter << '\n';
881  }
884  {
885  const auto arf = pow (betanom/nom0, 0.5/final_iter);
886  mfem::out << "Average reduction factor = " << arf << '\n';
887  }
889  {
890  mfem::out << "PCG: No convergence!" << '\n';
891  }
892
893  final_norm = sqrt(betanom);
894
895  Monitor(final_iter, final_norm, r, x, true);
896 }
897
898 void CG(const Operator &A, const Vector &b, Vector &x,
899  int print_iter, int max_num_iter,
900  double RTOLERANCE, double ATOLERANCE)
901 {
902  MFEM_PERF_FUNCTION;
903
904  CGSolver cg;
905  cg.SetPrintLevel(print_iter);
906  cg.SetMaxIter(max_num_iter);
907  cg.SetRelTol(sqrt(RTOLERANCE));
908  cg.SetAbsTol(sqrt(ATOLERANCE));
909  cg.SetOperator(A);
910  cg.Mult(b, x);
911 }
912
913 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
914  int print_iter, int max_num_iter,
915  double RTOLERANCE, double ATOLERANCE)
916 {
917  MFEM_PERF_FUNCTION;
918
919  CGSolver pcg;
920  pcg.SetPrintLevel(print_iter);
921  pcg.SetMaxIter(max_num_iter);
922  pcg.SetRelTol(sqrt(RTOLERANCE));
923  pcg.SetAbsTol(sqrt(ATOLERANCE));
924  pcg.SetOperator(A);
925  pcg.SetPreconditioner(B);
926  pcg.Mult(b, x);
927 }
928
929
930 inline void GeneratePlaneRotation(double &dx, double &dy,
931  double &cs, double &sn)
932 {
933  if (dy == 0.0)
934  {
935  cs = 1.0;
936  sn = 0.0;
937  }
938  else if (fabs(dy) > fabs(dx))
939  {
940  double temp = dx / dy;
941  sn = 1.0 / sqrt( 1.0 + temp*temp );
942  cs = temp * sn;
943  }
944  else
945  {
946  double temp = dy / dx;
947  cs = 1.0 / sqrt( 1.0 + temp*temp );
948  sn = temp * cs;
949  }
950 }
951
952 inline void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
953 {
954  double temp = cs * dx + sn * dy;
955  dy = -sn * dx + cs * dy;
956  dx = temp;
957 }
958
959 inline void Update(Vector &x, int k, DenseMatrix &h, Vector &s,
960  Array<Vector*> &v)
961 {
962  Vector y(s);
963
964  // Backsolve:
965  for (int i = k; i >= 0; i--)
966  {
967  y(i) /= h(i,i);
968  for (int j = i - 1; j >= 0; j--)
969  {
970  y(j) -= h(j,i) * y(i);
971  }
972  }
973
974  for (int j = 0; j <= k; j++)
975  {
976  x.Add(y(j), *v[j]);
977  }
978 }
979
980 void GMRESSolver::Mult(const Vector &b, Vector &x) const
981 {
982  // Generalized Minimum Residual method following the algorithm
983  // on p. 20 of the SIAM Templates book.
984
985  int n = width;
986
987  DenseMatrix H(m+1, m);
988  Vector s(m+1), cs(m+1), sn(m+1);
989  Vector r(n), w(n);
990  Array<Vector *> v;
991
992  int i, j, k;
993
994  if (iterative_mode)
995  {
996  oper->Mult(x, r);
997  }
998  else
999  {
1000  x = 0.0;
1001  }
1002
1003  if (prec)
1004  {
1005  if (iterative_mode)
1006  {
1007  subtract(b, r, w);
1008  prec->Mult(w, r); // r = M (b - A x)
1009  }
1010  else
1011  {
1012  prec->Mult(b, r);
1013  }
1014  }
1015  else
1016  {
1017  if (iterative_mode)
1018  {
1019  subtract(b, r, r);
1020  }
1021  else
1022  {
1023  r = b;
1024  }
1025  }
1026  double beta = initial_norm = Norm(r); // beta = ||r||
1027  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1028
1029  final_norm = std::max(rel_tol*beta, abs_tol);
1030
1031  if (beta <= final_norm)
1032  {
1033  final_norm = beta;
1034  final_iter = 0;
1035  converged = true;
1036  j = 0;
1037  goto finish;
1038  }
1039
1041  {
1042  mfem::out << " Pass : " << setw(2) << 1
1043  << " Iteration : " << setw(3) << 0
1044  << " ||B r|| = " << beta
1045  << (print_options.first_and_last ? " ...\n" : "\n");
1046  }
1047
1048  Monitor(0, beta, r, x);
1049
1050  v.SetSize(m+1, NULL);
1051
1052  for (j = 1; j <= max_iter; )
1053  {
1054  if (v[0] == NULL) { v[0] = new Vector(n); }
1055  v[0]->Set(1.0/beta, r);
1056  s = 0.0; s(0) = beta;
1057
1058  for (i = 0; i < m && j <= max_iter; i++, j++)
1059  {
1060  if (prec)
1061  {
1062  oper->Mult(*v[i], r);
1063  prec->Mult(r, w); // w = M A v[i]
1064  }
1065  else
1066  {
1067  oper->Mult(*v[i], w);
1068  }
1069
1070  for (k = 0; k <= i; k++)
1071  {
1072  H(k,i) = Dot(w, *v[k]); // H(k,i) = w * v[k]
1073  w.Add(-H(k,i), *v[k]); // w -= H(k,i) * v[k]
1074  }
1075
1076  H(i+1,i) = Norm(w); // H(i+1,i) = ||w||
1077  MFEM_ASSERT(IsFinite(H(i+1,i)), "Norm(w) = " << H(i+1,i));
1078  if (v[i+1] == NULL) { v[i+1] = new Vector(n); }
1079  v[i+1]->Set(1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
1080
1081  for (k = 0; k < i; k++)
1082  {
1083  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
1084  }
1085
1086  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1087  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1088  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
1089
1090  const double resid = fabs(s(i+1));
1091  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1092
1093  if (resid <= final_norm)
1094  {
1095  Update(x, i, H, s, v);
1096  final_norm = resid;
1097  final_iter = j;
1098  converged = true;
1099  goto finish;
1100  }
1101
1103  {
1104  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1105  << " Iteration : " << setw(3) << j
1106  << " ||B r|| = " << resid << '\n';
1107  }
1108
1109  Monitor(j, resid, r, x);
1110  }
1111
1112  if (print_options.iterations && j <= max_iter)
1113  {
1114  mfem::out << "Restarting..." << '\n';
1115  }
1116
1117  Update(x, i-1, H, s, v);
1118
1119  oper->Mult(x, r);
1120  if (prec)
1121  {
1122  subtract(b, r, w);
1123  prec->Mult(w, r); // r = M (b - A x)
1124  }
1125  else
1126  {
1127  subtract(b, r, r);
1128  }
1129  beta = Norm(r); // beta = ||r||
1130  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1131  if (beta <= final_norm)
1132  {
1133  final_norm = beta;
1134  final_iter = j;
1135  converged = true;
1136  goto finish;
1137  }
1138  }
1139
1140  final_norm = beta;
1141  final_iter = max_iter;
1142  converged = false;
1143
1144 finish:
1146  {
1147  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1148  << " Iteration : " << setw(3) << final_iter
1149  << " ||B r|| = " << final_norm << '\n';
1150  }
1152  {
1153  mfem::out << "GMRES: Number of iterations: " << final_iter << '\n';
1154  }
1156  {
1157  mfem::out << "GMRES: No convergence!\n";
1158  }
1159
1160  Monitor(final_iter, final_norm, r, x, true);
1161
1162  for (i = 0; i < v.Size(); i++)
1163  {
1164  delete v[i];
1165  }
1166 }
1167
1168 void FGMRESSolver::Mult(const Vector &b, Vector &x) const
1169 {
1170  DenseMatrix H(m+1,m);
1171  Vector s(m+1), cs(m+1), sn(m+1);
1172  Vector r(b.Size());
1173
1174  int i, j, k;
1175
1176  if (iterative_mode)
1177  {
1178  oper->Mult(x, r);
1179  subtract(b,r,r);
1180  }
1181  else
1182  {
1183  x = 0.;
1184  r = b;
1185  }
1186  double beta = initial_norm = Norm(r); // beta = ||r||
1187  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1188
1189  final_norm = std::max(rel_tol*beta, abs_tol);
1190
1191  if (beta <= final_norm)
1192  {
1193  final_norm = beta;
1194  final_iter = 0;
1195  converged = true;
1196  return;
1197  }
1198
1200  {
1201  mfem::out << " Pass : " << setw(2) << 1
1202  << " Iteration : " << setw(3) << 0
1203  << " || r || = " << beta
1204  << (print_options.first_and_last ? " ...\n" : "\n");
1205  }
1206
1207  Monitor(0, beta, r, x);
1208
1209  Array<Vector*> v(m+1);
1210  Array<Vector*> z(m+1);
1211  for (i= 0; i<=m; i++)
1212  {
1213  v[i] = NULL;
1214  z[i] = NULL;
1215  }
1216
1217  j = 1;
1218  while (j <= max_iter)
1219  {
1220  if (v[0] == NULL) { v[0] = new Vector(b.Size()); }
1221  (*v[0]) = 0.0;
1222  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
1223  s = 0.0; s(0) = beta;
1224
1225  for (i = 0; i < m && j <= max_iter; i++, j++)
1226  {
1227
1228  if (z[i] == NULL) { z[i] = new Vector(b.Size()); }
1229  (*z[i]) = 0.0;
1230
1231  if (prec)
1232  {
1233  prec->Mult(*v[i], *z[i]);
1234  }
1235  else
1236  {
1237  (*z[i]) = (*v[i]);
1238  }
1239  oper->Mult(*z[i], r);
1240
1241  for (k = 0; k <= i; k++)
1242  {
1243  H(k,i) = Dot( r, *v[k]); // H(k,i) = r * v[k]
1244  r.Add(-H(k,i), (*v[k])); // r -= H(k,i) * v[k]
1245  }
1246
1247  H(i+1,i) = Norm(r); // H(i+1,i) = ||r||
1248  if (v[i+1] == NULL) { v[i+1] = new Vector(b.Size()); }
1249  (*v[i+1]) = 0.0;
1250  v[i+1] -> Add (1.0/H(i+1,i), r); // v[i+1] = r / H(i+1,i)
1251
1252  for (k = 0; k < i; k++)
1253  {
1254  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
1255  }
1256
1257  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1258  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1259  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
1260
1261  const double resid = fabs(s(i+1));
1262  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1264  resid <= final_norm))
1265  {
1266  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1267  << " Iteration : " << setw(3) << j
1268  << " || r || = " << resid << endl;
1269  }
1270  Monitor(j, resid, r, x, resid <= final_norm);
1271
1272  if (resid <= final_norm)
1273  {
1274  Update(x, i, H, s, z);
1275  final_norm = resid;
1276  final_iter = j;
1277  converged = true;
1278
1279  if (print_options.summary)
1280  {
1281  mfem::out << "FGMRES: Number of iterations: " << final_iter << '\n';
1282  }
1283
1284  for (i= 0; i<=m; i++)
1285  {
1286  if (v[i]) { delete v[i]; }
1287  if (z[i]) { delete z[i]; }
1288  }
1289  return;
1290  }
1291  }
1292
1294  {
1295  mfem::out << "Restarting..." << endl;
1296  }
1297
1298  Update(x, i-1, H, s, z);
1299
1300  oper->Mult(x, r);
1301  subtract(b,r,r);
1302  beta = Norm(r);
1303  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1304  if (beta <= final_norm)
1305  {
1306  final_norm = beta;
1307  final_iter = j;
1308  converged = true;
1309
1310  break;
1311  }
1312  }
1313
1314  // Clean buffers up
1315  for (i = 0; i <= m; i++)
1316  {
1317  if (v[i]) { delete v[i]; }
1318  if (z[i]) { delete z[i]; }
1319  }
1320  converged = false;
1321
1322  // Note: j is off by one when we arrive here
1324  {
1325  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
1326  << " Iteration : " << setw(3) << j-1
1327  << " || r || = " << final_norm << endl;
1328  }
1330  {
1331  mfem::out << "FGMRES: Number of iterations: " << j-1 << '\n';
1332  }
1334  {
1335  mfem::out << "FGMRES: No convergence!\n";
1336  }
1337 }
1338
1339
1340 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
1341  int &max_iter, int m, double &tol, double atol, int printit)
1342 {
1343  MFEM_PERF_FUNCTION;
1344
1345  GMRESSolver gmres;
1346  gmres.SetPrintLevel(printit);
1347  gmres.SetMaxIter(max_iter);
1348  gmres.SetKDim(m);
1349  gmres.SetRelTol(sqrt(tol));
1350  gmres.SetAbsTol(sqrt(atol));
1351  gmres.SetOperator(A);
1352  gmres.SetPreconditioner(M);
1353  gmres.Mult(b, x);
1354  max_iter = gmres.GetNumIterations();
1355  tol = gmres.GetFinalNorm()*gmres.GetFinalNorm();
1356  return gmres.GetConverged();
1357 }
1358
1359 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
1360  int print_iter, int max_num_iter, int m, double rtol, double atol)
1361 {
1362  GMRES(A, x, b, B, max_num_iter, m, rtol, atol, print_iter);
1363 }
1364
1365
1367 {
1368  p.SetSize(width);
1369  phat.SetSize(width);
1370  s.SetSize(width);
1371  shat.SetSize(width);
1372  t.SetSize(width);
1373  v.SetSize(width);
1374  r.SetSize(width);
1375  rtilde.SetSize(width);
1376 }
1377
1378 void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
1379 {
1380  // BiConjugate Gradient Stabilized method following the algorithm
1381  // on p. 27 of the SIAM Templates book.
1382
1383  int i;
1384  double resid, tol_goal;
1385  double rho_1, rho_2=1.0, alpha=1.0, beta, omega=1.0;
1386
1387  if (iterative_mode)
1388  {
1389  oper->Mult(x, r);
1390  subtract(b, r, r); // r = b - A x
1391  }
1392  else
1393  {
1394  x = 0.0;
1395  r = b;
1396  }
1397  rtilde = r;
1398
1399  resid = initial_norm = Norm(r);
1400  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1402  {
1403  mfem::out << " Iteration : " << setw(3) << 0
1404  << " ||r|| = " << resid << (print_options.first_and_last ? " ...\n" : "\n");
1405  }
1406
1407  Monitor(0, resid, r, x);
1408
1409  tol_goal = std::max(resid*rel_tol, abs_tol);
1410
1411  if (resid <= tol_goal)
1412  {
1413  final_norm = resid;
1414  final_iter = 0;
1415  converged = true;
1416  return;
1417  }
1418
1419  for (i = 1; i <= max_iter; i++)
1420  {
1421  rho_1 = Dot(rtilde, r);
1422  if (rho_1 == 0)
1423  {
1425  {
1426  mfem::out << " Iteration : " << setw(3) << i
1427  << " ||r|| = " << resid << '\n';
1428  }
1429
1430  Monitor(i, resid, r, x);
1431
1432  final_norm = resid;
1433  final_iter = i;
1434  converged = false;
1436  {
1437  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1438  }
1439  if (print_options.warnings)
1440  {
1441  mfem::out << "BiCGStab: No convergence!\n";
1442  }
1443  return;
1444  }
1445  if (i == 1)
1446  {
1447  p = r;
1448  }
1449  else
1450  {
1451  beta = (rho_1/rho_2) * (alpha/omega);
1452  add(p, -omega, v, p); // p = p - omega * v
1453  add(r, beta, p, p); // p = r + beta * p
1454  }
1455  if (prec)
1456  {
1457  prec->Mult(p, phat); // phat = M^{-1} * p
1458  }
1459  else
1460  {
1461  phat = p;
1462  }
1463  oper->Mult(phat, v); // v = A * phat
1464  alpha = rho_1 / Dot(rtilde, v);
1465  add(r, -alpha, v, s); // s = r - alpha * v
1466  resid = Norm(s);
1467  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1468  if (resid < tol_goal)
1469  {
1470  x.Add(alpha, phat); // x = x + alpha * phat
1472  {
1473  mfem::out << " Iteration : " << setw(3) << i
1474  << " ||s|| = " << resid << '\n';
1475  }
1476  final_norm = resid;
1477  final_iter = i;
1478  converged = true;
1480  {
1481  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1482  }
1483  return;
1484  }
1486  {
1487  mfem::out << " Iteration : " << setw(3) << i
1488  << " ||s|| = " << resid;
1489  }
1490  Monitor(i, resid, r, x);
1491  if (prec)
1492  {
1493  prec->Mult(s, shat); // shat = M^{-1} * s
1494  }
1495  else
1496  {
1497  shat = s;
1498  }
1499  oper->Mult(shat, t); // t = A * shat
1500  omega = Dot(t, s) / Dot(t, t);
1501  x.Add(alpha, phat); // x += alpha * phat
1502  x.Add(omega, shat); // x += omega * shat
1503  add(s, -omega, t, r); // r = s - omega * t
1504
1505  rho_2 = rho_1;
1506  resid = Norm(r);
1507  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1509  {
1510  mfem::out << " ||r|| = " << resid << '\n';
1511  }
1512  Monitor(i, resid, r, x);
1513  if (resid < tol_goal)
1514  {
1515  final_norm = resid;
1516  final_iter = i;
1517  converged = true;
1519  {
1520  mfem::out << " Iteration : " << setw(3) << i
1521  << " ||r|| = " << resid << '\n';
1522  }
1524  {
1525  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1526  }
1527  return;
1528  }
1529  if (omega == 0)
1530  {
1531  final_norm = resid;
1532  final_iter = i;
1533  converged = false;
1535  {
1536  mfem::out << " Iteration : " << setw(3) << i
1537  << " ||r|| = " << resid << '\n';
1538  }
1540  {
1541  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1542  }
1543  if (print_options.warnings)
1544  {
1545  mfem::out << "BiCGStab: No convergence!\n";
1546  }
1547  return;
1548  }
1549  }
1550
1551  final_norm = resid;
1552  final_iter = max_iter;
1553  converged = false;
1554
1556  {
1557  mfem::out << " Iteration : " << setw(3) << final_iter
1558  << " ||r|| = " << resid << '\n';
1559  }
1561  {
1562  mfem::out << "BiCGStab: Number of iterations: " << final_iter << '\n';
1563  }
1564  if (print_options.warnings)
1565  {
1566  mfem::out << "BiCGStab: No convergence!\n";
1567  }
1568 }
1569
1570 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
1571  int &max_iter, double &tol, double atol, int printit)
1572 {
1573  BiCGSTABSolver bicgstab;
1574  bicgstab.SetPrintLevel(printit);
1575  bicgstab.SetMaxIter(max_iter);
1576  bicgstab.SetRelTol(sqrt(tol));
1577  bicgstab.SetAbsTol(sqrt(atol));
1578  bicgstab.SetOperator(A);
1579  bicgstab.SetPreconditioner(M);
1580  bicgstab.Mult(b, x);
1581  max_iter = bicgstab.GetNumIterations();
1582  tol = bicgstab.GetFinalNorm()*bicgstab.GetFinalNorm();
1583  return bicgstab.GetConverged();
1584 }
1585
1586 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
1587  int print_iter, int max_num_iter, double rtol, double atol)
1588 {
1589  BiCGSTAB(A, x, b, B, max_num_iter, rtol, atol, print_iter);
1590 }
1591
1592
1594 {
1596  v0.SetSize(width);
1597  v1.SetSize(width);
1598  w0.SetSize(width);
1599  w1.SetSize(width);
1600  q.SetSize(width);
1601  if (prec)
1602  {
1603  u1.SetSize(width);
1604  }
1605
1606  v0.UseDevice(true);
1607  v1.UseDevice(true);
1608  w0.UseDevice(true);
1609  w1.UseDevice(true);
1610  q.UseDevice(true);
1611  u1.UseDevice(true);
1612 }
1613
1614 void MINRESSolver::Mult(const Vector &b, Vector &x) const
1615 {
1616  // Based on the MINRES algorithm on p. 86, Fig. 6.9 in
1617  // "Iterative Krylov Methods for Large Linear Systems",
1618  // by Henk A. van der Vorst, 2003.
1619  // Extended to support an SPD preconditioner.
1620
1621  b.UseDevice(true);
1622  x.UseDevice(true);
1623
1624  int it;
1625  double beta, eta, gamma0, gamma1, sigma0, sigma1;
1626  double alpha, delta, rho1, rho2, rho3, norm_goal;
1627  Vector *z = (prec) ? &u1 : &v1;
1628
1629  converged = true;
1630
1631  if (!iterative_mode)
1632  {
1633  v1 = b;
1634  x = 0.;
1635  }
1636  else
1637  {
1638  oper->Mult(x, v1);
1639  subtract(b, v1, v1);
1640  }
1641
1642  if (prec)
1643  {
1644  prec->Mult(v1, u1);
1645  }
1646  eta = beta = initial_norm = sqrt(Dot(*z, v1));
1647  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1648  gamma0 = gamma1 = 1.;
1649  sigma0 = sigma1 = 0.;
1650
1651  norm_goal = std::max(rel_tol*eta, abs_tol);
1652
1653  if (eta <= norm_goal)
1654  {
1655  it = 0;
1656  goto loop_end;
1657  }
1658
1660  {
1661  mfem::out << "MINRES: iteration " << setw(3) << 0 << ": ||r||_B = "
1662  << eta << (print_options.first_and_last ? " ..." : "") << '\n';
1663  }
1664  Monitor(0, eta, *z, x);
1665
1666  for (it = 1; it <= max_iter; it++)
1667  {
1668  v1 /= beta;
1669  if (prec)
1670  {
1671  u1 /= beta;
1672  }
1673  oper->Mult(*z, q);
1674  alpha = Dot(*z, q);
1675  MFEM_ASSERT(IsFinite(alpha), "alpha = " << alpha);
1676  if (it > 1) // (v0 == 0) for (it == 1)
1677  {
1678  q.Add(-beta, v0);
1679  }
1680  add(q, -alpha, v1, v0);
1681
1682  delta = gamma1*alpha - gamma0*sigma1*beta;
1683  rho3 = sigma0*beta;
1684  rho2 = sigma1*alpha + gamma0*gamma1*beta;
1685  if (!prec)
1686  {
1687  beta = Norm(v0);
1688  }
1689  else
1690  {
1691  prec->Mult(v0, q);
1692  beta = sqrt(Dot(v0, q));
1693  }
1694  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1695  rho1 = hypot(delta, beta);
1696
1697  if (it == 1)
1698  {
1699  w0.Set(1./rho1, *z); // (w0 == 0) and (w1 == 0)
1700  }
1701  else if (it == 2)
1702  {
1703  add(1./rho1, *z, -rho2/rho1, w1, w0); // (w0 == 0)
1704  }
1705  else
1706  {
1707  add(-rho3/rho1, w0, -rho2/rho1, w1, w0);
1708  w0.Add(1./rho1, *z);
1709  }
1710
1711  gamma0 = gamma1;
1712  gamma1 = delta/rho1;
1713
1714  x.Add(gamma1*eta, w0);
1715
1716  sigma0 = sigma1;
1717  sigma1 = beta/rho1;
1718
1719  eta = -sigma1*eta;
1720  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1721
1722  if (fabs(eta) <= norm_goal)
1723  {
1724  goto loop_end;
1725  }
1726
1728  {
1729  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1730  << fabs(eta) << '\n';
1731  }
1732  Monitor(it, fabs(eta), *z, x);
1733
1734  if (prec)
1735  {
1736  Swap(u1, q);
1737  }
1738  Swap(v0, v1);
1739  Swap(w0, w1);
1740  }
1741  converged = false;
1742  it--;
1743
1744 loop_end:
1745  final_iter = it;
1746  final_norm = fabs(eta);
1747
1749  {
1750  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1751  << fabs(eta) << '\n';
1752  }
1753
1755  {
1756  mfem::out << "MINRES: Number of iterations: " << setw(3) << final_iter << '\n';
1757  }
1758
1759  Monitor(final_iter, final_norm, *z, x, true);
1760
1761  // if (print_options.iteration_details || (!converged && print_options.errors))
1762  // {
1763  // oper->Mult(x, v1);
1764  // subtract(b, v1, v1);
1765  // if (prec)
1766  // {
1767  // prec->Mult(v1, u1);
1768  // }
1769  // eta = sqrt(Dot(*z, v1));
1770  // mfem::out << "MINRES: iteration " << setw(3) << it << '\n'
1771  // << " ||r||_B = " << eta << " (re-computed)" << '\n';
1772  // }
1773
1774  if (!converged && (print_options.warnings))
1775  {
1776  mfem::out << "MINRES: No convergence!\n";
1777  }
1778 }
1779
1780 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it,
1781  int max_it, double rtol, double atol)
1782 {
1783  MFEM_PERF_FUNCTION;
1784
1785  MINRESSolver minres;
1786  minres.SetPrintLevel(print_it);
1787  minres.SetMaxIter(max_it);
1788  minres.SetRelTol(sqrt(rtol));
1789  minres.SetAbsTol(sqrt(atol));
1790  minres.SetOperator(A);
1791  minres.Mult(b, x);
1792 }
1793
1794 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
1795  int print_it, int max_it, double rtol, double atol)
1796 {
1797  MINRESSolver minres;
1798  minres.SetPrintLevel(print_it);
1799  minres.SetMaxIter(max_it);
1800  minres.SetRelTol(sqrt(rtol));
1801  minres.SetAbsTol(sqrt(atol));
1802  minres.SetOperator(A);
1803  minres.SetPreconditioner(B);
1804  minres.Mult(b, x);
1805 }
1806
1807
1809 {
1810  oper = &op;
1811  height = op.Height();
1812  width = op.Width();
1813  MFEM_ASSERT(height == width, "square Operator is required.");
1814
1815  r.SetSize(width);
1816  c.SetSize(width);
1817 }
1818
1819 void NewtonSolver::Mult(const Vector &b, Vector &x) const
1820 {
1821  MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
1822  MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver).");
1823
1824  int it;
1825  double norm0, norm, norm_goal;
1826  const bool have_b = (b.Size() == Height());
1827
1828  if (!iterative_mode)
1829  {
1830  x = 0.0;
1831  }
1832
1833  ProcessNewState(x);
1834
1835  oper->Mult(x, r);
1836  if (have_b)
1837  {
1838  r -= b;
1839  }
1840
1841  norm0 = norm = initial_norm = Norm(r);
1843  {
1844  mfem::out << "Newton iteration " << setw(2) << 0
1845  << " : ||r|| = " << norm << "...\n";
1846  }
1847  norm_goal = std::max(rel_tol*norm, abs_tol);
1848
1849  prec->iterative_mode = false;
1850
1851  // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b]
1852  for (it = 0; true; it++)
1853  {
1854  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
1856  {
1857  mfem::out << "Newton iteration " << setw(2) << it
1858  << " : ||r|| = " << norm;
1859  if (it > 0)
1860  {
1861  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
1862  }
1863  mfem::out << '\n';
1864  }
1865  Monitor(it, norm, r, x);
1866
1867  if (norm <= norm_goal)
1868  {
1869  converged = true;
1870  break;
1871  }
1872
1873  if (it >= max_iter)
1874  {
1875  converged = false;
1876  break;
1877  }
1878
1879  grad = &oper->GetGradient(x);
1880  prec->SetOperator(*grad);
1881
1882  if (lin_rtol_type)
1883  {
1884  AdaptiveLinRtolPreSolve(x, it, norm);
1885  }
1886
1887  prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b]
1888
1889  if (lin_rtol_type)
1890  {
1891  AdaptiveLinRtolPostSolve(c, r, it, norm);
1892  }
1893
1894  const double c_scale = ComputeScalingFactor(x, b);
1895  if (c_scale == 0.0)
1896  {
1897  converged = false;
1898  break;
1899  }
1900  add(x, -c_scale, c, x);
1901
1902  ProcessNewState(x);
1903
1904  oper->Mult(x, r);
1905  if (have_b)
1906  {
1907  r -= b;
1908  }
1909  norm = Norm(r);
1910  }
1911
1912  final_iter = it;
1913  final_norm = norm;
1914
1917  {
1918  mfem::out << "Newton: Number of iterations: " << final_iter << '\n'
1919  << " ||r|| = " << final_norm << '\n';
1920  }
1922  {
1923  mfem::out << "Newton: No convergence!\n";
1924  }
1925 }
1926
1928  const double rtol0,
1929  const double rtol_max,
1930  const double alpha_,
1931  const double gamma_)
1932 {
1933  lin_rtol_type = type;
1934  lin_rtol0 = rtol0;
1935  lin_rtol_max = rtol_max;
1936  this->alpha = alpha_;
1937  this->gamma = gamma_;
1938 }
1939
1941  const int it,
1942  const double fnorm) const
1943 {
1944  // Assume that when adaptive linear solver relative tolerance is activated,
1945  // we are working with an iterative solver.
1946  auto iterative_solver = static_cast<IterativeSolver *>(prec);
1947  // Adaptive linear solver relative tolerance
1948  double eta;
1949  // Safeguard threshold
1950  double sg_threshold = 0.1;
1951
1952  if (it == 0)
1953  {
1954  eta = lin_rtol0;
1955  }
1956  else
1957  {
1958  if (lin_rtol_type == 1)
1959  {
1960  // eta = gamma * abs(||F(x1)|| - ||F(x0) + DF(x0) s0||) / ||F(x0)||
1961  eta = gamma * abs(fnorm - lnorm_last) / fnorm_last;
1962  }
1963  else if (lin_rtol_type == 2)
1964  {
1965  // eta = gamma * (||F(x1)|| / ||F(x0)||)^alpha
1966  eta = gamma * pow(fnorm / fnorm_last, alpha);
1967  }
1968  else
1969  {
1970  MFEM_ABORT("Unknown adaptive linear solver rtol version");
1971  }
1972
1973  // Safeguard rtol from "oversolving" ?!
1974  const double sg_eta = gamma * pow(eta_last, alpha);
1975  if (sg_eta > sg_threshold) { eta = std::max(eta, sg_eta); }
1976  }
1977
1978  eta = std::min(eta, lin_rtol_max);
1979  iterative_solver->SetRelTol(eta);
1980  eta_last = eta;
1982  {
1983  mfem::out << "Eisenstat-Walker rtol = " << eta << "\n";
1984  }
1985 }
1986
1988  const Vector &b,
1989  const int it,
1990  const double fnorm) const
1991 {
1992  fnorm_last = fnorm;
1993
1994  // If version 1 is chosen, the true linear residual norm has to be computed
1995  // and in most cases we can only retrieve the preconditioned linear residual
1996  // norm.
1997  if (lin_rtol_type == 1)
1998  {
1999  // lnorm_last = ||F(x0) + DF(x0) s0||
2000  Vector linres(x.Size());
2001  grad->Mult(x, linres);
2002  linres -= b;
2003  lnorm_last = Norm(linres);
2004  }
2005 }
2006
2007 void LBFGSSolver::Mult(const Vector &b, Vector &x) const
2008 {
2009  MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");
2010
2011  // Quadrature points that are checked for negative Jacobians etc.
2012  Vector sk, rk, yk, rho, alpha;
2013
2014  // r - r_{k+1}, c - descent direction
2015  sk.SetSize(width); // x_{k+1}-x_k
2016  rk.SetSize(width); // nabla(f(x_{k}))
2017  yk.SetSize(width); // r_{k+1}-r_{k}
2018  rho.SetSize(m); // 1/(dot(yk,sk)
2019  alpha.SetSize(m); // rhok*sk'*c
2020  int last_saved_id = -1;
2021
2022  int it;
2023  double norm0, norm, norm_goal;
2024  const bool have_b = (b.Size() == Height());
2025
2026  if (!iterative_mode)
2027  {
2028  x = 0.0;
2029  }
2030
2031  ProcessNewState(x);
2032
2033  // r = F(x)-b
2034  oper->Mult(x, r);
2035  if (have_b) { r -= b; }
2036
2037  c = r; // initial descent direction
2038
2039  norm0 = norm = initial_norm = Norm(r);
2041  {
2042  mfem::out << "LBFGS iteration " << setw(2) << 0
2043  << " : ||r|| = " << norm << "...\n";
2044  }
2045  norm_goal = std::max(rel_tol*norm, abs_tol);
2046  for (it = 0; true; it++)
2047  {
2048  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
2050  {
2051  mfem::out << "LBFGS iteration " << it
2052  << " : ||r|| = " << norm;
2053  if (it > 0)
2054  {
2055  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
2056  }
2057  mfem::out << '\n';
2058  }
2059
2060  if (norm <= norm_goal)
2061  {
2062  converged = true;
2063  break;
2064  }
2065
2066  if (it >= max_iter)
2067  {
2068  converged = false;
2069  break;
2070  }
2071
2072  rk = r;
2073  const double c_scale = ComputeScalingFactor(x, b);
2074  if (c_scale == 0.0)
2075  {
2076  converged = false;
2077  break;
2078  }
2079  add(x, -c_scale, c, x); // x_{k+1} = x_k - c_scale*c
2080
2081  ProcessNewState(x);
2082
2083  oper->Mult(x, r);
2084  if (have_b)
2085  {
2086  r -= b;
2087  }
2088
2089  // LBFGS - construct descent direction
2090  subtract(r, rk, yk); // yk = r_{k+1} - r_{k}
2091  sk = c; sk *= -c_scale; //sk = x_{k+1} - x_{k} = -c_scale*c
2092  const double gamma = Dot(sk, yk)/Dot(yk, yk);
2093
2094  // Save last m vectors
2095  last_saved_id = (last_saved_id == m-1) ? 0 : last_saved_id+1;
2096  *skArray[last_saved_id] = sk;
2097  *ykArray[last_saved_id] = yk;
2098
2099  c = r;
2100  for (int i = last_saved_id; i > -1; i--)
2101  {
2102  rho(i) = 1.0/Dot((*skArray[i]),(*ykArray[i]));
2103  alpha(i) = rho(i)*Dot((*skArray[i]),c);
2104  add(c, -alpha(i), (*ykArray[i]), c);
2105  }
2106  if (it > m-1)
2107  {
2108  for (int i = m-1; i > last_saved_id; i--)
2109  {
2110  rho(i) = 1./Dot((*skArray[i]), (*ykArray[i]));
2111  alpha(i) = rho(i)*Dot((*skArray[i]),c);
2112  add(c, -alpha(i), (*ykArray[i]), c);
2113  }
2114  }
2115
2116  c *= gamma; // scale search direction
2117  if (it > m-1)
2118  {
2119  for (int i = last_saved_id+1; i < m ; i++)
2120  {
2121  double betai = rho(i)*Dot((*ykArray[i]), c);
2122  add(c, alpha(i)-betai, (*skArray[i]), c);
2123  }
2124  }
2125  for (int i = 0; i < last_saved_id+1 ; i++)
2126  {
2127  double betai = rho(i)*Dot((*ykArray[i]), c);
2128  add(c, alpha(i)-betai, (*skArray[i]), c);
2129  }
2130
2131  norm = Norm(r);
2132  }
2133
2134  final_iter = it;
2135  final_norm = norm;
2136
2139  {
2140  mfem::out << "LBFGS: Number of iterations: " << final_iter << '\n'
2141  << " ||r|| = " << final_norm << '\n';
2142  }
2144  {
2145  mfem::out << "LBFGS: No convergence!\n";
2146  }
2147 }
2148
2149 int aGMRES(const Operator &A, Vector &x, const Vector &b,
2150  const Operator &M, int &max_iter,
2151  int m_max, int m_min, int m_step, double cf,
2152  double &tol, double &atol, int printit)
2153 {
2154  int n = A.Width();
2155
2156  int m = m_max;
2157
2158  DenseMatrix H(m+1,m);
2159  Vector s(m+1), cs(m+1), sn(m+1);
2160  Vector w(n), av(n);
2161
2162  double r1, resid;
2163  int i, j, k;
2164
2165  M.Mult(b,w);
2166  double normb = w.Norml2(); // normb = ||M b||
2167  if (normb == 0.0)
2168  {
2169  normb = 1;
2170  }
2171
2172  Vector r(n);
2173  A.Mult(x, r);
2174  subtract(b,r,w);
2175  M.Mult(w, r); // r = M (b - A x)
2176  double beta = r.Norml2(); // beta = ||r||
2177
2178  resid = beta / normb;
2179
2180  if (resid * resid <= tol)
2181  {
2182  tol = resid * resid;
2183  max_iter = 0;
2184  return 0;
2185  }
2186
2187  if (printit)
2188  {
2189  mfem::out << " Pass : " << setw(2) << 1
2190  << " Iteration : " << setw(3) << 0
2191  << " (r, r) = " << beta*beta << '\n';
2192  }
2193
2194  tol *= (normb*normb);
2195  tol = (atol > tol) ? atol : tol;
2196
2197  m = m_max;
2198  Array<Vector *> v(m+1);
2199  for (i= 0; i<=m; i++)
2200  {
2201  v[i] = new Vector(n);
2202  (*v[i]) = 0.0;
2203  }
2204
2205  j = 1;
2206  while (j <= max_iter)
2207  {
2208  (*v[0]) = 0.0;
2209  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
2210  s = 0.0; s(0) = beta;
2211
2212  r1 = beta;
2213
2214  for (i = 0; i < m && j <= max_iter; i++)
2215  {
2216  A.Mult((*v[i]),av);
2217  M.Mult(av,w); // w = M A v[i]
2218
2219  for (k = 0; k <= i; k++)
2220  {
2221  H(k,i) = w * (*v[k]); // H(k,i) = w * v[k]
2222  w.Add(-H(k,i), (*v[k])); // w -= H(k,i) * v[k]
2223  }
2224
2225  H(i+1,i) = w.Norml2(); // H(i+1,i) = ||w||
2226  (*v[i+1]) = 0.0;
2227  v[i+1] -> Add (1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
2228
2229  for (k = 0; k < i; k++)
2230  {
2231  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
2232  }
2233
2234  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
2235  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
2236  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
2237
2238  resid = fabs(s(i+1));
2239  if (printit)
2240  {
2241  mfem::out << " Pass : " << setw(2) << j
2242  << " Iteration : " << setw(3) << i+1
2243  << " (r, r) = " << resid*resid << '\n';
2244  }
2245
2246  if ( resid*resid < tol)
2247  {
2248  Update(x, i, H, s, v);
2249  tol = resid * resid;
2250  max_iter = j;
2251  for (i= 0; i<=m; i++)
2252  {
2253  delete v[i];
2254  }
2255  return 0;
2256  }
2257  }
2258
2259  if (printit)
2260  {
2261  mfem::out << "Restarting..." << '\n';
2262  }
2263
2264  Update(x, i-1, H, s, v);
2265
2266  A.Mult(x, r);
2267  subtract(b,r,w);
2268  M.Mult(w, r); // r = M (b - A x)
2269  beta = r.Norml2(); // beta = ||r||
2270  if ( resid*resid < tol)
2271  {
2272  tol = resid * resid;
2273  max_iter = j;
2274  for (i= 0; i<=m; i++)
2275  {
2276  delete v[i];
2277  }
2278  return 0;
2279  }
2280
2281  if (beta/r1 > cf)
2282  {
2283  if (m - m_step >= m_min)
2284  {
2285  m -= m_step;
2286  }
2287  else
2288  {
2289  m = m_max;
2290  }
2291  }
2292
2293  j++;
2294  }
2295
2296  tol = resid * resid;
2297  for (i= 0; i<=m; i++)
2298  {
2299  delete v[i];
2300  }
2301  return 1;
2302 }
2303
2305  const Operator *C_,
2306  const Operator *D_)
2307  : C(C_), D(D_), c_e(NULL), d_lo(NULL), d_hi(NULL), x_lo(NULL), x_hi(NULL),
2308  input_size(insize)
2309 {
2310  if (C) { MFEM_ASSERT(C->Width() == input_size, "Wrong width of C."); }
2311  if (D) { MFEM_ASSERT(D->Width() == input_size, "Wrong width of D."); }
2312 }
2313
2315 {
2316  MFEM_ASSERT(C, "The C operator is unspecified -- can't set constraints.");
2317  MFEM_ASSERT(c.Size() == C->Height(), "Wrong size of the constraint.");
2318
2319  c_e = &c;
2320 }
2321
2323  const Vector &dh)
2324 {
2325  MFEM_ASSERT(D, "The D operator is unspecified -- can't set constraints.");
2326  MFEM_ASSERT(dl.Size() == D->Height() && dh.Size() == D->Height(),
2327  "Wrong size of the constraint.");
2328
2329  d_lo = &dl; d_hi = &dh;
2330 }
2331
2333 {
2334  MFEM_ASSERT(xl.Size() == input_size && xh.Size() == input_size,
2335  "Wrong size of the constraint.");
2336
2337  x_lo = &xl; x_hi = &xh;
2338 }
2339
2341 {
2342  int m = 0;
2343  if (C) { m += C->Height(); }
2344  if (D) { m += D->Height(); }
2345  return m;
2346 }
2347
2349 {
2350  if (print_options.warnings)
2351  {
2352  MFEM_WARNING("Objective functional is ignored as SLBQP always minimizes"
2353  "the l2 norm of (x - x_target).");
2354  }
2355  MFEM_ASSERT(prob.GetC(), "Linear constraint is not set.");
2356  MFEM_ASSERT(prob.GetC()->Height() == 1, "Solver expects scalar constraint.");
2357
2358  problem = &prob;
2359 }
2360
2361 void SLBQPOptimizer::SetBounds(const Vector &lo_, const Vector &hi_)
2362 {
2363  lo.SetDataAndSize(lo_.GetData(), lo_.Size());
2364  hi.SetDataAndSize(hi_.GetData(), hi_.Size());
2365 }
2366
2368 {
2369  w.SetDataAndSize(w_.GetData(), w_.Size());
2370  a = a_;
2371 }
2372
2373 inline void SLBQPOptimizer::print_iteration(int it, double r, double l) const
2374 {
2376  {
2377  mfem::out << "SLBQP iteration " << it << ": residual = " << r
2378  << ", lambda = " << l << '\n';
2379  }
2380 }
2381
2382 void SLBQPOptimizer::Mult(const Vector& xt, Vector& x) const
2383 {
2384  // Based on code provided by Denis Ridzal, dridzal@sandia.gov.
2385  // Algorithm adapted from Dai and Fletcher, "New Algorithms for
2386  // Singly Linearly Constrained Quadratic Programs Subject to Lower
2387  // and Upper Bounds", Numerical Analysis Report NA/216, 2003.
2388
2389  // Set some algorithm-specific constants and temporaries.
2390  int nclip = 0;
2391  double l = 0;
2392  double llow = 0;
2393  double lupp = 0;
2394  double lnew = 0;
2395  double dl = 2;
2396  double r = 0;
2397  double rlow = 0;
2398  double rupp = 0;
2399  double s = 0;
2400
2401  const double smin = 0.1;
2402
2403  const double tol = max(abs_tol, rel_tol*a);
2404
2405  // *** Start bracketing phase of SLBQP ***
2407  {
2408  mfem::out << "SLBQP bracketing phase" << '\n';
2409  }
2410
2411  // Solve QP with fixed Lagrange multiplier
2412  r = initial_norm = solve(l,xt,x,nclip);
2413  print_iteration(nclip, r, l);
2414
2415
2416  // If x=xt was already within bounds and satisfies the linear
2417  // constraint, then we already have the solution.
2418  if (fabs(r) <= tol)
2419  {
2420  converged = true;
2421  goto slbqp_done;
2422  }
2423
2424  if (r < 0)
2425  {
2426  llow = l; rlow = r; l = l + dl;
2427
2428  // Solve QP with fixed Lagrange multiplier
2429  r = solve(l,xt,x,nclip);
2430  print_iteration(nclip, r, l);
2431
2432  while ((r < 0) && (nclip < max_iter))
2433  {
2434  llow = l;
2435  s = rlow/r - 1.0;
2436  if (s < smin) { s = smin; }
2437  dl = dl + dl/s;
2438  l = l + dl;
2439
2440  // Solve QP with fixed Lagrange multiplier
2441  r = solve(l,xt,x,nclip);
2442  print_iteration(nclip, r, l);
2443  }
2444
2445  lupp = l; rupp = r;
2446  }
2447  else
2448  {
2449  lupp = l; rupp = r; l = l - dl;
2450
2451  // Solve QP with fixed Lagrange multiplier
2452  r = solve(l,xt,x,nclip);
2453  print_iteration(nclip, r, l);
2454
2455  while ((r > 0) && (nclip < max_iter))
2456  {
2457  lupp = l;
2458  s = rupp/r - 1.0;
2459  if (s < smin) { s = smin; }
2460  dl = dl + dl/s;
2461  l = l - dl;
2462
2463  // Solve QP with fixed Lagrange multiplier
2464  r = solve(l,xt,x,nclip);
2465  print_iteration(nclip, r, l);
2466  }
2467
2468  llow = l; rlow = r;
2469  }
2470
2471  // *** Stop bracketing phase of SLBQP ***
2472
2473
2474  // *** Start secant phase of SLBQP ***
2476  {
2477  mfem::out << "SLBQP secant phase" << '\n';
2478  }
2479
2480  s = 1.0 - rlow/rupp; dl = dl/s; l = lupp - dl;
2481
2482  // Solve QP with fixed Lagrange multiplier
2483  r = solve(l,xt,x,nclip);
2484  print_iteration(nclip, r, l);
2485
2486  while ( (fabs(r) > tol) && (nclip < max_iter) )
2487  {
2488  if (r > 0)
2489  {
2490  if (s <= 2.0)
2491  {
2492  lupp = l; rupp = r; s = 1.0 - rlow/rupp;
2493  dl = (lupp - llow)/s; l = lupp - dl;
2494  }
2495  else
2496  {
2497  s = rupp/r - 1.0;
2498  if (s < smin) { s = smin; }
2499  dl = (lupp - l)/s;
2500  lnew = 0.75*llow + 0.25*l;
2501  if (lnew < l-dl) { lnew = l-dl; }
2502  lupp = l; rupp = r; l = lnew;
2503  s = (lupp - llow)/(lupp - l);
2504  }
2505
2506  }
2507  else
2508  {
2509  if (s >= 2.0)
2510  {
2511  llow = l; rlow = r; s = 1.0 - rlow/rupp;
2512  dl = (lupp - llow)/s; l = lupp - dl;
2513  }
2514  else
2515  {
2516  s = rlow/r - 1.0;
2517  if (s < smin) { s = smin; }
2518  dl = (l - llow)/s;
2519  lnew = 0.75*lupp + 0.25*l;
2520  if (lnew < l+dl) { lnew = l+dl; }
2521  llow = l; rlow = r; l = lnew;
2522  s = (lupp - llow)/(lupp - l);
2523  }
2524  }
2525
2526  // Solve QP with fixed Lagrange multiplier
2527  r = solve(l,xt,x,nclip);
2528  print_iteration(nclip, r, l);
2529  }
2530
2531  // *** Stop secant phase of SLBQP ***
2532  converged = (fabs(r) <= tol);
2533
2534 slbqp_done:
2535
2536  final_iter = nclip;
2537  final_norm = r;
2538
2541  {
2542  mfem::out << "SLBQP: Number of iterations: " << final_iter << '\n'
2543  << " lambda = " << l << '\n'
2544  << " ||r|| = " << final_norm << '\n';
2545  }
2547  {
2548  mfem::out << "SLBQP: No convergence!" << '\n';
2549  }
2550 }
2551
2552 struct WeightMinHeap
2553 {
2554  const std::vector<double> &w;
2555  std::vector<size_t> c;
2556  std::vector<int> loc;
2557
2558  WeightMinHeap(const std::vector<double> &w_) : w(w_)
2559  {
2560  c.reserve(w.size());
2561  loc.resize(w.size());
2562  for (size_t i=0; i<w.size(); ++i) { push(i); }
2563  }
2564
2565  size_t percolate_up(size_t pos, double val)
2566  {
2567  for (; pos > 0 && w[c[(pos-1)/2]] > val; pos = (pos-1)/2)
2568  {
2569  c[pos] = c[(pos-1)/2];
2570  loc[c[(pos-1)/2]] = pos;
2571  }
2572  return pos;
2573  }
2574
2575  size_t percolate_down(size_t pos, double val)
2576  {
2577  while (2*pos+1 < c.size())
2578  {
2579  size_t left = 2*pos+1;
2580  size_t right = left+1;
2581  size_t tgt;
2582  if (right < c.size() && w[c[right]] < w[c[left]]) { tgt = right; }
2583  else { tgt = left; }
2584  if (w[c[tgt]] < val)
2585  {
2586  c[pos] = c[tgt];
2587  loc[c[tgt]] = pos;
2588  pos = tgt;
2589  }
2590  else
2591  {
2592  break;
2593  }
2594  }
2595  return pos;
2596  }
2597
2598  void push(size_t i)
2599  {
2600  double val = w[i];
2601  c.push_back(0);
2602  size_t pos = c.size()-1;
2603  pos = percolate_up(pos, val);
2604  c[pos] = i;
2605  loc[i] = pos;
2606  }
2607
2608  int pop()
2609  {
2610  size_t i = c[0];
2611  size_t j = c.back();
2612  c.pop_back();
2613  // Mark as removed
2614  loc[i] = -1;
2615  if (c.empty()) { return i; }
2616  double val = w[j];
2617  size_t pos = 0;
2618  pos = percolate_down(pos, val);
2619  c[pos] = j;
2620  loc[j] = pos;
2621  return i;
2622  }
2623
2624  void update(size_t i)
2625  {
2626  size_t pos = loc[i];
2627  double val = w[i];
2628  pos = percolate_up(pos, val);
2629  pos = percolate_down(pos, val);
2630  c[pos] = i;
2631  loc[i] = pos;
2632  }
2633
2634  bool picked(size_t i)
2635  {
2636  return loc[i] < 0;
2637  }
2638 };
2639
2641 {
2642  int n = C.Width();
2643  // Scale rows by reciprocal of diagonal and take absolute value
2644  Vector D;
2645  C.GetDiag(D);
2646  int *I = C.GetI();
2647  int *J = C.GetJ();
2648  double *V = C.GetData();
2649  for (int i=0; i<n; ++i)
2650  {
2651  for (int j=I[i]; j<I[i+1]; ++j)
2652  {
2653  V[j] = abs(V[j]/D[i]);
2654  }
2655  }
2656
2657  std::vector<double> w(n, 0.0);
2658  for (int k=0; k<n; ++k)
2659  {
2660  // Find all neighbors i of k
2661  for (int ii=I[k]; ii<I[k+1]; ++ii)
2662  {
2663  int i = J[ii];
2664  // Find value of (i,k)
2665  double C_ik = 0.0;
2666  for (int kk=I[i]; kk<I[i+1]; ++kk)
2667  {
2668  if (J[kk] == k)
2669  {
2670  C_ik = V[kk];
2671  break;
2672  }
2673  }
2674  for (int jj=I[k]; jj<I[k+1]; ++jj)
2675  {
2676  int j = J[jj];
2677  if (j == k) { continue; }
2678  double C_kj = V[jj];
2679  bool ij_exists = false;
2680  for (int jj2=I[i]; jj2<I[i+1]; ++jj2)
2681  {
2682  if (J[jj2] == j)
2683  {
2684  ij_exists = true;
2685  break;
2686  }
2687  }
2688  if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2689  }
2690  }
2691  w[k] = sqrt(w[k]);
2692  }
2693
2694  WeightMinHeap w_heap(w);
2695
2696  // Compute ordering
2697  p.SetSize(n);
2698  for (int ii=0; ii<n; ++ii)
2699  {
2700  int pi = w_heap.pop();
2701  p[ii] = pi;
2702  w[pi] = -1;
2703  for (int kk=I[pi]; kk<I[pi+1]; ++kk)
2704  {
2705  int k = J[kk];
2706  if (w_heap.picked(k)) { continue; }
2707  // Recompute weight
2708  w[k] = 0.0;
2709  // Find all neighbors i of k
2710  for (int ii2=I[k]; ii2<I[k+1]; ++ii2)
2711  {
2712  int i = J[ii2];
2713  if (w_heap.picked(i)) { continue; }
2714  // Find value of (i,k)
2715  double C_ik = 0.0;
2716  for (int kk2=I[i]; kk2<I[i+1]; ++kk2)
2717  {
2718  if (J[kk2] == k)
2719  {
2720  C_ik = V[kk2];
2721  break;
2722  }
2723  }
2724  for (int jj=I[k]; jj<I[k+1]; ++jj)
2725  {
2726  int j = J[jj];
2727  if (j == k || w_heap.picked(j)) { continue; }
2728  double C_kj = V[jj];
2729  bool ij_exists = false;
2730  for (int jj2=I[i]; jj2<I[i+1]; ++jj2)
2731  {
2732  if (J[jj2] == j)
2733  {
2734  ij_exists = true;
2735  break;
2736  }
2737  }
2738  if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2739  }
2740  }
2741  w[k] = sqrt(w[k]);
2742  w_heap.update(k);
2743  }
2744  }
2745 }
2746
2747 BlockILU::BlockILU(int block_size_,
2748  Reordering reordering_,
2749  int k_fill_)
2750  : Solver(0),
2751  block_size(block_size_),
2752  k_fill(k_fill_),
2753  reordering(reordering_)
2754 { }
2755
2757  int block_size_,
2758  Reordering reordering_,
2759  int k_fill_)
2760  : BlockILU(block_size_, reordering_, k_fill_)
2761 {
2762  SetOperator(op);
2763 }
2764
2766 {
2767  const SparseMatrix *A = NULL;
2768 #ifdef MFEM_USE_MPI
2769  const HypreParMatrix *A_par = dynamic_cast<const HypreParMatrix *>(&op);
2770  SparseMatrix A_par_diag;
2771  if (A_par != NULL)
2772  {
2773  A_par->GetDiag(A_par_diag);
2774  A = &A_par_diag;
2775  }
2776 #endif
2777  if (A == NULL)
2778  {
2779  A = dynamic_cast<const SparseMatrix *>(&op);
2780  if (A == NULL)
2781  {
2782  MFEM_ABORT("BlockILU must be created with a SparseMatrix or HypreParMatrix");
2783  }
2784  }
2785  height = op.Height();
2786  width = op.Width();
2787  MFEM_ASSERT(A->Finalized(), "Matrix must be finalized.");
2788  CreateBlockPattern(*A);
2789  Factorize();
2790 }
2791
2792 void BlockILU::CreateBlockPattern(const SparseMatrix &A)
2793 {
2794  MFEM_VERIFY(k_fill == 0, "Only block ILU(0) is currently supported.");
2795  if (A.Height() % block_size != 0)
2796  {
2797  MFEM_ABORT("BlockILU: block size must evenly divide the matrix size");
2798  }
2799
2800  int nrows = A.Height();
2801  const int *I = A.GetI();
2802  const int *J = A.GetJ();
2803  const double *V = A.GetData();
2804  int nnz = 0;
2805  int nblockrows = nrows / block_size;
2806
2807  std::vector<std::set<int>> unique_block_cols(nblockrows);
2808
2809  for (int iblock = 0; iblock < nblockrows; ++iblock)
2810  {
2811  for (int bi = 0; bi < block_size; ++bi)
2812  {
2813  int i = iblock * block_size + bi;
2814  for (int k = I[i]; k < I[i + 1]; ++k)
2815  {
2816  unique_block_cols[iblock].insert(J[k] / block_size);
2817  }
2818  }
2819  nnz += unique_block_cols[iblock].size();
2820  }
2821
2822  if (reordering != Reordering::NONE)
2823  {
2824  SparseMatrix C(nblockrows, nblockrows);
2825  for (int iblock = 0; iblock < nblockrows; ++iblock)
2826  {
2827  for (int jblock : unique_block_cols[iblock])
2828  {
2829  for (int bi = 0; bi < block_size; ++bi)
2830  {
2831  int i = iblock * block_size + bi;
2832  for (int k = I[i]; k < I[i + 1]; ++k)
2833  {
2834  int j = J[k];
2835  if (j >= jblock * block_size && j < (jblock + 1) * block_size)
2836  {
2837  C.Add(iblock, jblock, V[k]*V[k]);
2838  }
2839  }
2840  }
2841  }
2842  }
2843  C.Finalize(false);
2844  double *CV = C.GetData();
2845  for (int i=0; i<C.NumNonZeroElems(); ++i)
2846  {
2847  CV[i] = sqrt(CV[i]);
2848  }
2849
2850  switch (reordering)
2851  {
2854  break;
2855  default:
2856  MFEM_ABORT("BlockILU: unknown reordering")
2857  }
2858  }
2859  else
2860  {
2861  // No reordering: permutation is identity
2862  P.SetSize(nblockrows);
2863  for (int i=0; i<nblockrows; ++i)
2864  {
2865  P[i] = i;
2866  }
2867  }
2868
2869  // Compute inverse permutation
2870  Pinv.SetSize(nblockrows);
2871  for (int i=0; i<nblockrows; ++i)
2872  {
2873  Pinv[P[i]] = i;
2874  }
2875
2876  // Permute columns
2877  std::vector<std::vector<int>> unique_block_cols_perminv(nblockrows);
2878  for (int i=0; i<nblockrows; ++i)
2879  {
2880  std::vector<int> &cols = unique_block_cols_perminv[i];
2881  for (int j : unique_block_cols[P[i]])
2882  {
2883  cols.push_back(Pinv[j]);
2884  }
2885  std::sort(cols.begin(), cols.end());
2886  }
2887
2888  ID.SetSize(nblockrows);
2889  IB.SetSize(nblockrows + 1);
2890  IB[0] = 0;
2891  JB.SetSize(nnz);
2892  AB.SetSize(block_size, block_size, nnz);
2893  DB.SetSize(block_size, block_size, nblockrows);
2894  AB = 0.0;
2895  DB = 0.0;
2896  ipiv.SetSize(block_size*nblockrows);
2897  int counter = 0;
2898
2899  for (int iblock = 0; iblock < nblockrows; ++iblock)
2900  {
2901  int iblock_perm = P[iblock];
2902  for (int jblock : unique_block_cols_perminv[iblock])
2903  {
2904  int jblock_perm = P[jblock];
2905  if (iblock == jblock)
2906  {
2907  ID[iblock] = counter;
2908  }
2909  JB[counter] = jblock;
2910  for (int bi = 0; bi < block_size; ++bi)
2911  {
2912  int i = iblock_perm*block_size + bi;
2913  for (int k = I[i]; k < I[i + 1]; ++k)
2914  {
2915  int j = J[k];
2916  if (j >= jblock_perm*block_size && j < (jblock_perm + 1)*block_size)
2917  {
2918  int bj = j - jblock_perm*block_size;
2919  double val = V[k];
2920  AB(bi, bj, counter) = val;
2921  // Extract the diagonal
2922  if (iblock == jblock)
2923  {
2924  DB(bi, bj, iblock) = val;
2925  }
2926  }
2927  }
2928  }
2929  ++counter;
2930  }
2931  IB[iblock + 1] = counter;
2932  }
2933 }
2934
2935 void BlockILU::Factorize()
2936 {
2937  int nblockrows = Height()/block_size;
2938
2939  // Precompute LU factorization of diagonal blocks
2940  for (int i=0; i<nblockrows; ++i)
2941  {
2942  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2943  factorization.Factor(block_size);
2944  }
2945
2946  // Note: we use UseExternalData to extract submatrices from the tensor AB
2947  // instead of the DenseTensor call operator, because the call operator does
2948  // not allow for two simultaneous submatrix views into the same tensor
2949  DenseMatrix A_ik, A_ij, A_kj;
2950  // Loop over block rows (starting with second block row)
2951  for (int i=1; i<nblockrows; ++i)
2952  {
2953  // Find all nonzeros to the left of the diagonal in row i
2954  for (int kk=IB[i]; kk<IB[i+1]; ++kk)
2955  {
2956  int k = JB[kk];
2957  // Make sure we're still to the left of the diagonal
2958  if (k == i) { break; }
2959  if (k > i)
2960  {
2961  MFEM_ABORT("Matrix must be sorted with nonzero diagonal");
2962  }
2963  LUFactors A_kk_inv(DB.GetData(k), &ipiv[k*block_size]);
2964  A_ik.UseExternalData(&AB(0,0,kk), block_size, block_size);
2965  // A_ik = A_ik * A_kk^{-1}
2966  A_kk_inv.RightSolve(block_size, block_size, A_ik.GetData());
2967  // Modify everything to the right of k in row i
2968  for (int jj=kk+1; jj<IB[i+1]; ++jj)
2969  {
2970  int j = JB[jj];
2971  if (j <= k) { continue; } // Superfluous because JB is sorted?
2972  A_ij.UseExternalData(&AB(0,0,jj), block_size, block_size);
2973  for (int ll=IB[k]; ll<IB[k+1]; ++ll)
2974  {
2975  int l = JB[ll];
2976  if (l == j)
2977  {
2978  A_kj.UseExternalData(&AB(0,0,ll), block_size, block_size);
2979  // A_ij = A_ij - A_ik*A_kj;
2980  AddMult_a(-1.0, A_ik, A_kj, A_ij);
2981  // If we need to, update diagonal factorization
2982  if (j == i)
2983  {
2984  DB(i) = A_ij;
2985  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2986  factorization.Factor(block_size);
2987  }
2988  break;
2989  }
2990  }
2991  }
2992  }
2993  }
2994 }
2995
2996 void BlockILU::Mult(const Vector &b, Vector &x) const
2997 {
2998  MFEM_ASSERT(height > 0, "BlockILU(0) preconditioner is not constructed");
2999  int nblockrows = Height()/block_size;
3000  y.SetSize(Height());
3001
3002  DenseMatrix B;
3003  Vector yi, yj, xi, xj;
3004  Vector tmp(block_size);
3005  // Forward substitute to solve Ly = b
3006  // Implicitly, L has identity on the diagonal
3007  y = 0.0;
3008  for (int i=0; i<nblockrows; ++i)
3009  {
3010  yi.SetDataAndSize(&y[i*block_size], block_size);
3011  for (int ib=0; ib<block_size; ++ib)
3012  {
3013  yi[ib] = b[ib + P[i]*block_size];
3014  }
3015  for (int k=IB[i]; k<ID[i]; ++k)
3016  {
3017  int j = JB[k];
3018  const DenseMatrix &L_ij = AB(k);
3019  yj.SetDataAndSize(&y[j*block_size], block_size);
3020  // y_i = y_i - L_ij*y_j
3021  L_ij.AddMult_a(-1.0, yj, yi);
3022  }
3023  }
3024  // Backward substitution to solve Ux = y
3025  for (int i=nblockrows-1; i >= 0; --i)
3026  {
3027  xi.SetDataAndSize(&x[P[i]*block_size], block_size);
3028  for (int ib=0; ib<block_size; ++ib)
3029  {
3030  xi[ib] = y[ib + i*block_size];
3031  }
3032  for (int k=ID[i]+1; k<IB[i+1]; ++k)
3033  {
3034  int j = JB[k];
3035  const DenseMatrix &U_ij = AB(k);
3036  xj.SetDataAndSize(&x[P[j]*block_size], block_size);
3037  // x_i = x_i - U_ij*x_j
3038  U_ij.AddMult_a(-1.0, xj, xi);
3039  }
3040  LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
3041  // x_i = D_ii^{-1} x_i
3042  A_ii_inv.Solve(block_size, 1, xi.GetData());
3043  }
3044 }
3045
3046
3048  int it, double norm, const Vector &r, bool final)
3049 {
3050  if (!ess_dofs_list) { return; }
3051
3052  double bc_norm_squared = 0.0;
3053  r.HostRead();
3055  for (int i = 0; i < ess_dofs_list->Size(); i++)
3056  {
3057  const double r_entry = r((*ess_dofs_list)[i]);
3058  bc_norm_squared += r_entry*r_entry;
3059  }
3060  bool print = true;
3061 #ifdef MFEM_USE_MPI
3062  MPI_Comm comm = iter_solver->GetComm();
3063  if (comm != MPI_COMM_NULL)
3064  {
3065  double glob_bc_norm_squared = 0.0;
3066  MPI_Reduce(&bc_norm_squared, &glob_bc_norm_squared, 1, MPI_DOUBLE,
3067  MPI_SUM, 0, comm);
3068  bc_norm_squared = glob_bc_norm_squared;
3069  int rank;
3070  MPI_Comm_rank(comm, &rank);
3071  print = (rank == 0);
3072  }
3073 #endif
3074  if ((it == 0 || final || bc_norm_squared > 0.0) && print)
3075  {
3076  mfem::out << " ResidualBCMonitor : b.c. residual norm = "
3077  << sqrt(bc_norm_squared) << endl;
3078  }
3079 }
3080
3081
3082 #ifdef MFEM_USE_SUITESPARSE
3083
3085 {
3086  mat = NULL;
3087  Numeric = NULL;
3088  AI = AJ = NULL;
3089  if (!use_long_ints)
3090  {
3091  umfpack_di_defaults(Control);
3092  }
3093  else
3094  {
3095  umfpack_dl_defaults(Control);
3096  }
3097 }
3098
3100 {
3101  void *Symbolic;
3102
3103  if (Numeric)
3104  {
3105  if (!use_long_ints)
3106  {
3107  umfpack_di_free_numeric(&Numeric);
3108  }
3109  else
3110  {
3111  umfpack_dl_free_numeric(&Numeric);
3112  }
3113  }
3114
3115  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
3116  MFEM_VERIFY(mat, "not a SparseMatrix");
3117
3118  // UMFPack requires that the column-indices in mat corresponding to each
3119  // row be sorted.
3120  // Generally, this will modify the ordering of the entries of mat.
3122
3123  height = mat->Height();
3124  width = mat->Width();
3125  MFEM_VERIFY(width == height, "not a square matrix");
3126
3127  const int * Ap = mat->HostReadI();
3128  const int * Ai = mat->HostReadJ();
3129  const double * Ax = mat->HostReadData();
3130
3131  if (!use_long_ints)
3132  {
3133  int status = umfpack_di_symbolic(width, width, Ap, Ai, Ax, &Symbolic,
3134  Control, Info);
3135  if (status < 0)
3136  {
3137  umfpack_di_report_info(Control, Info);
3138  umfpack_di_report_status(Control, status);
3139  mfem_error("UMFPackSolver::SetOperator :"
3140  " umfpack_di_symbolic() failed!");
3141  }
3142
3143  status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric,
3144  Control, Info);
3145  if (status < 0)
3146  {
3147  umfpack_di_report_info(Control, Info);
3148  umfpack_di_report_status(Control, status);
3149  mfem_error("UMFPackSolver::SetOperator :"
3150  " umfpack_di_numeric() failed!");
3151  }
3152  umfpack_di_free_symbolic(&Symbolic);
3153  }
3154  else
3155  {
3156  SuiteSparse_long status;
3157
3158  delete [] AJ;
3159  delete [] AI;
3160  AI = new SuiteSparse_long[width + 1];
3161  AJ = new SuiteSparse_long[Ap[width]];
3162  for (int i = 0; i <= width; i++)
3163  {
3164  AI[i] = (SuiteSparse_long)(Ap[i]);
3165  }
3166  for (int i = 0; i < Ap[width]; i++)
3167  {
3168  AJ[i] = (SuiteSparse_long)(Ai[i]);
3169  }
3170
3171  status = umfpack_dl_symbolic(width, width, AI, AJ, Ax, &Symbolic,
3172  Control, Info);
3173  if (status < 0)
3174  {
3175  umfpack_dl_report_info(Control, Info);
3176  umfpack_dl_report_status(Control, status);
3177  mfem_error("UMFPackSolver::SetOperator :"
3178  " umfpack_dl_symbolic() failed!");
3179  }
3180
3181  status = umfpack_dl_numeric(AI, AJ, Ax, Symbolic, &Numeric,
3182  Control, Info);
3183  if (status < 0)
3184  {
3185  umfpack_dl_report_info(Control, Info);
3186  umfpack_dl_report_status(Control, status);
3187  mfem_error("UMFPackSolver::SetOperator :"
3188  " umfpack_dl_numeric() failed!");
3189  }
3190  umfpack_dl_free_symbolic(&Symbolic);
3191  }
3192 }
3193
3194 void UMFPackSolver::Mult(const Vector &b, Vector &x) const
3195 {
3196  if (mat == NULL)
3197  mfem_error("UMFPackSolver::Mult : matrix is not set!"
3198  " Call SetOperator first!");
3199  b.HostRead();
3200  x.HostReadWrite();
3201  if (!use_long_ints)
3202  {
3203  int status =
3204  umfpack_di_solve(UMFPACK_At, mat->HostReadI(), mat->HostReadJ(),
3205  mat->HostReadData(), x.HostWrite(), b.HostRead(),
3206  Numeric, Control, Info);
3207  umfpack_di_report_info(Control, Info);
3208  if (status < 0)
3209  {
3210  umfpack_di_report_status(Control, status);
3211  mfem_error("UMFPackSolver::Mult : umfpack_di_solve() failed!");
3212  }
3213  }
3214  else
3215  {
3216  SuiteSparse_long status =
3217  umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->HostReadData(),
3218  x.HostWrite(), b.HostRead(), Numeric, Control,
3219  Info);
3220  umfpack_dl_report_info(Control, Info);
3221  if (status < 0)
3222  {
3223  umfpack_dl_report_status(Control, status);
3224  mfem_error("UMFPackSolver::Mult : umfpack_dl_solve() failed!");
3225  }
3226  }
3227 }
3228
3230 {
3231  if (mat == NULL)
3232  mfem_error("UMFPackSolver::MultTranspose : matrix is not set!"
3233  " Call SetOperator first!");
3234  b.HostRead();
3235  x.HostReadWrite();
3236  if (!use_long_ints)
3237  {
3238  int status =
3239  umfpack_di_solve(UMFPACK_A, mat->HostReadI(), mat->HostReadJ(),
3240  mat->HostReadData(), x.HostWrite(), b.HostRead(),
3241  Numeric, Control, Info);
3242  umfpack_di_report_info(Control, Info);
3243  if (status < 0)
3244  {
3245  umfpack_di_report_status(Control, status);
3246  mfem_error("UMFPackSolver::MultTranspose :"
3247  " umfpack_di_solve() failed!");
3248  }
3249  }
3250  else
3251  {
3252  SuiteSparse_long status =
3253  umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->HostReadData(),
3254  x.HostWrite(), b.HostRead(), Numeric, Control,
3255  Info);
3256  umfpack_dl_report_info(Control, Info);
3257  if (status < 0)
3258  {
3259  umfpack_dl_report_status(Control, status);
3260  mfem_error("UMFPackSolver::MultTranspose :"
3261  " umfpack_dl_solve() failed!");
3262  }
3263  }
3264 }
3265
3267 {
3268  delete [] AJ;
3269  delete [] AI;
3270  if (Numeric)
3271  {
3272  if (!use_long_ints)
3273  {
3274  umfpack_di_free_numeric(&Numeric);
3275  }
3276  else
3277  {
3278  umfpack_dl_free_numeric(&Numeric);
3279  }
3280  }
3281 }
3282
3284 {
3285  klu_defaults(&Common);
3286 }
3287
3289 {
3290  if (Numeric)
3291  {
3292  MFEM_ASSERT(Symbolic != 0,
3293  "Had Numeric pointer in KLU, but not Symbolic");
3294  klu_free_symbolic(&Symbolic, &Common);
3295  Symbolic = 0;
3296  klu_free_numeric(&Numeric, &Common);
3297  Numeric = 0;
3298  }
3299
3300  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
3301  MFEM_VERIFY(mat != NULL, "not a SparseMatrix");
3302
3303  // KLU requires that the column-indices in mat corresponding to each row be
3304  // sorted. Generally, this will modify the ordering of the entries of mat.
3306
3307  height = mat->Height();
3308  width = mat->Width();
3309  MFEM_VERIFY(width == height, "not a square matrix");
3310
3311  int * Ap = mat->GetI();
3312  int * Ai = mat->GetJ();
3313  double * Ax = mat->GetData();
3314
3315  Symbolic = klu_analyze( height, Ap, Ai, &Common);
3316  Numeric = klu_factor(Ap, Ai, Ax, Symbolic, &Common);
3317 }
3318
3319 void KLUSolver::Mult(const Vector &b, Vector &x) const
3320 {
3321  MFEM_VERIFY(mat != NULL,
3322  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3323
3324  int n = mat->Height();
3325  int numRhs = 1;
3326  // Copy B into X, so we can pass it in and overwrite it.
3327  x = b;
3328  // Solve the transpose, since KLU thinks the matrix is compressed column
3329  // format.
3330  klu_tsolve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
3331 }
3332
3333 void KLUSolver::MultTranspose(const Vector &b, Vector &x) const
3334 {
3335  MFEM_VERIFY(mat != NULL,
3336  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3337
3338  int n = mat->Height();
3339  int numRhs = 1;
3340  // Copy B into X, so we can pass it in and overwrite it.
3341  x = b;
3342  // Solve the regular matrix, not the transpose, since KLU thinks the matrix
3343  // is compressed column format.
3344  klu_solve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
3345 }
3346
3348 {
3349  klu_free_symbolic (&Symbolic, &Common) ;
3350  klu_free_numeric (&Numeric, &Common) ;
3351  Symbolic = 0;
3352  Numeric = 0;
3353 }
3354
3355 #endif // MFEM_USE_SUITESPARSE
3356
3358  const SparseMatrix &block_dof_)
3359  : Solver(A.NumRows()), block_dof(const_cast<SparseMatrix&>(block_dof_)),
3360  block_solvers(new DenseMatrixInverse[block_dof.NumRows()])
3361 {
3362  DenseMatrix sub_A;
3363  for (int i = 0; i < block_dof.NumRows(); ++i)
3364  {
3365  local_dofs.MakeRef(block_dof.GetRowColumns(i), block_dof.RowSize(i));
3366  sub_A.SetSize(local_dofs.Size());
3367  A.GetSubMatrix(local_dofs, local_dofs, sub_A);
3368  block_solvers[i].SetOperator(sub_A);
3369  }
3370 }
3371
3372 void DirectSubBlockSolver::Mult(const Vector &x, Vector &y) const
3373 {
3374  y.SetSize(x.Size());
3375  y = 0.0;
3376
3377  for (int i = 0; i < block_dof.NumRows(); ++i)
3378  {
3379  local_dofs.MakeRef(block_dof.GetRowColumns(i), block_dof.RowSize(i));
3380  x.GetSubVector(local_dofs, sub_rhs);
3381  sub_sol.SetSize(local_dofs.Size());
3382  block_solvers[i].Mult(sub_rhs, sub_sol);
3383  y.AddElementVector(local_dofs, sub_sol);
3384  }
3385 }
3386
3387 void ProductSolver::Mult(const Vector & x, Vector & y) const
3388 {
3389  y.SetSize(x.Size());
3390  y = 0.0;
3391  S0->Mult(x, y);
3392
3393  Vector z(x.Size());
3394  z = 0.0;
3395  A->Mult(y, z);
3396  add(-1.0, z, 1.0, x, z); // z = (I - A * S0) x
3397
3398  Vector S1z(x.Size());
3399  S1z = 0.0;
3400  S1->Mult(z, S1z);
3401  y += S1z;
3402 }
3403
3404 void ProductSolver::MultTranspose(const Vector & x, Vector & y) const
3405 {
3406  y.SetSize(x.Size());
3407  y = 0.0;
3408  S1->MultTranspose(x, y);
3409
3410  Vector z(x.Size());
3411  z = 0.0;
3412  A->MultTranspose(y, z);
3413  add(-1.0, z, 1.0, x, z); // z = (I - A^T * S1^T) x
3414
3415  Vector S0Tz(x.Size());
3416  S0Tz = 0.0;
3417  S0->MultTranspose(z, S0Tz);
3418  y += S0Tz;
3419 }
3420
3422  : Solver(0, false), global_size(-1)
3423 #ifdef MFEM_USE_MPI
3424  , parallel(false)
3425 #endif
3426 { }
3427
3428 #ifdef MFEM_USE_MPI
3429 OrthoSolver::OrthoSolver(MPI_Comm mycomm_)
3430  : Solver(0, false), mycomm(mycomm_), global_size(-1), parallel(true) { }
3431 #endif
3432
3434 {
3435  solver = &s;
3436  height = s.Height();
3437  width = s.Width();
3438  MFEM_VERIFY(height == width, "Solver must be a square Operator!");
3439  global_size = -1; // lazy evaluated
3440 }
3441
3443 {
3444  MFEM_VERIFY(solver, "Solver hasn't been set, call SetSolver() first.");
3445  solver->SetOperator(op);
3446  height = solver->Height();
3447  width = solver->Width();
3448  MFEM_VERIFY(height == width, "Solver must be a square Operator!");
3449  global_size = -1; // lazy evaluated
3450 }
3451
3452 void OrthoSolver::Mult(const Vector &b, Vector &x) const
3453 {
3454  MFEM_VERIFY(solver, "Solver hasn't been set, call SetSolver() first.");
3455  MFEM_VERIFY(height == solver->Height(),
3456  "solver was modified externally! call SetSolver() again!");
3457  MFEM_VERIFY(height == b.Size(), "incompatible input Vector size!");
3458  MFEM_VERIFY(height == x.Size(), "incompatible output Vector size!");
3459
3460  // Orthogonalize input
3461  Orthogonalize(b, b_ortho);
3462
3463  // Propagate iterative_mode to the solver:
3464  solver->iterative_mode = iterative_mode;
3465
3466  // Apply the Solver
3467  solver->Mult(b_ortho, x);
3468
3469  // Orthogonalize output
3470  Orthogonalize(x, x);
3471 }
3472
3473 void OrthoSolver::Orthogonalize(const Vector &v, Vector &v_ortho) const
3474 {
3475  if (global_size == -1)
3476  {
3477  global_size = height;
3478 #ifdef MFEM_USE_MPI
3479  if (parallel)
3480  {
3481  MPI_Allreduce(MPI_IN_PLACE, &global_size, 1, HYPRE_MPI_BIG_INT,
3482  MPI_SUM, mycomm);
3483  }
3484 #endif
3485  }
3486
3487  // TODO: GPU/device implementation
3488
3489  double global_sum = v.Sum();
3490
3491 #ifdef MFEM_USE_MPI
3492  if (parallel)
3493  {
3494  MPI_Allreduce(MPI_IN_PLACE, &global_sum, 1, MPI_DOUBLE, MPI_SUM, mycomm);
3495  }
3496 #endif
3497
3498  double ratio = global_sum / static_cast<double>(global_size);
3499  v_ortho.SetSize(v.Size());
3500  v.HostRead();
3501  v_ortho.HostWrite();
3502  for (int i = 0; i < v_ortho.Size(); ++i)
3503  {
3504  v_ortho(i) = v(i) - ratio;
3505  }
3506 }
3507
3508 #ifdef MFEM_USE_MPI
3510  HypreParMatrix *aux_map,
3511  bool op_is_symmetric,
3512  bool own_aux_map)
3513  : Solver(op.NumRows()), aux_map_(aux_map, own_aux_map)
3514 {
3515  aux_system_.Reset(RAP(&op, aux_map));
3516  aux_system_.As<HypreParMatrix>()->EliminateZeroRows();
3517  aux_smoother_.Reset(new HypreSmoother(*aux_system_.As<HypreParMatrix>()));
3518  aux_smoother_.As<HypreSmoother>()->SetOperatorSymmetry(op_is_symmetric);
3519 }
3520
3521 void AuxSpaceSmoother::Mult(const Vector &x, Vector &y, bool transpose) const
3522 {
3523  Vector aux_rhs(aux_map_->NumCols());
3524  aux_map_->MultTranspose(x, aux_rhs);
3525
3526  Vector aux_sol(aux_rhs.Size());
3527  if (transpose)
3528  {
3529  aux_smoother_->MultTranspose(aux_rhs, aux_sol);
3530  }
3531  else
3532  {
3533  aux_smoother_->Mult(aux_rhs, aux_sol);
3534  }
3535
3536  y.SetSize(aux_map_->NumRows());
3537  aux_map_->Mult(aux_sol, y);
3538 }
3539 #endif // MFEM_USE_MPI
3540
3541 #ifdef MFEM_USE_LAPACK
3542 // LAPACK routines for NNLSSolver
3543 extern "C" void
3544 dormqr_(char *, char *, int *, int *, int *, double *, int*, double *,
3545  double *, int *, double *, int*, int*);
3546
3547 extern "C" void
3548 dgeqrf_(int *, int *, double *, int *, double *, double *, int *, int *);
3549
3550 extern "C" void
3551 dgemv_(char *, int *, int *, double *, double *, int *, double *, int *,
3552  double *, double *, int *);
3553
3554 extern "C" void
3555 dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m, int *n,
3556  double *alpha, double *a, int *lda, double *b, int *ldb);
3557
3559  : Solver(0), mat(nullptr), const_tol_(1.0e-14), min_nnz_(0),
3560  max_nnz_(0), verbosity_(0), res_change_termination_tol_(1.0e-4),
3561  zero_tol_(1.0e-14), rhs_delta_(1.0e-11), n_outer_(100000),
3562  n_inner_(100000), nStallCheck_(100), normalize_(true),
3563  NNLS_qrres_on_(false), qr_residual_mode_(QRresidualMode::hybrid)
3564 {}
3565
3567 {
3568  mat = dynamic_cast<const DenseMatrix*>(&op);
3569  MFEM_VERIFY(mat, "NNLSSolver operator must be of type DenseMatrix");
3570
3571  // The size of this operator is that of the transpose of op.
3572  height = op.Width();
3573  width = op.Height();
3574
3575  row_scaling_.SetSize(mat->NumRows());
3576  row_scaling_ = 1.0;
3577 }
3578
3579 void NNLSSolver::SetQRResidualMode(const QRresidualMode qr_residual_mode)
3580 {
3581  qr_residual_mode_ = qr_residual_mode;
3582  if (qr_residual_mode_ == QRresidualMode::on)
3583  {
3584  NNLS_qrres_on_ = true;
3585  }
3586 }
3587
3588 void NNLSSolver::NormalizeConstraints(Vector& rhs_lb, Vector& rhs_ub) const
3589 {
3590  // Scale everything so that rescaled half gap is the same for all constraints
3591  const int m = mat->NumRows();
3592
3593  MFEM_VERIFY(rhs_lb.Size() == m && rhs_ub.Size() == m, "");
3594
3595  Vector rhs_avg = rhs_ub;
3596  rhs_avg += rhs_lb;
3597  rhs_avg *= 0.5;
3598
3599  Vector rhs_halfgap = rhs_ub;
3600  rhs_halfgap -= rhs_lb;
3601  rhs_halfgap *= 0.5;
3602
3603  Vector rhs_avg_glob = rhs_avg;
3604  Vector rhs_halfgap_glob = rhs_halfgap;
3605  Vector halfgap_target(m);
3606  halfgap_target = 1.0e3 * const_tol_;
3607
3608  row_scaling_.SetSize(m);
3609
3610  for (int i=0; i<m; ++i)
3611  {
3612  const double s = halfgap_target(i) / rhs_halfgap_glob(i);
3613  row_scaling_[i] = s;
3614
3615  rhs_lb(i) = (rhs_avg(i) * s) - halfgap_target(i);
3616  rhs_ub(i) = (rhs_avg(i) * s) + halfgap_target(i);
3617  }
3618 }
3619
3620 void NNLSSolver::Mult(const Vector &w, Vector &sol) const
3621 {
3622  MFEM_VERIFY(mat, "NNLSSolver operator must be of type DenseMatrix");
3623  Vector rhs_ub(mat->NumRows());
3624  mat->Mult(w, rhs_ub);
3625  rhs_ub *= row_scaling_;
3626
3627  Vector rhs_lb(rhs_ub);
3628  Vector rhs_Gw(rhs_ub);
3629
3630  for (int i=0; i<rhs_ub.Size(); ++i)
3631  {
3632  rhs_lb(i) -= rhs_delta_;
3633  rhs_ub(i) += rhs_delta_;
3634  }
3635
3636  if (normalize_) { NormalizeConstraints(rhs_lb, rhs_ub); }
3637  Solve(rhs_lb, rhs_ub, sol);
3638
3639  if (verbosity_ > 1)
3640  {
3641  int nnz = 0;
3642  for (int i=0; i<sol.Size(); ++i)
3643  {
3644  if (sol(i) != 0.0)
3645  {
3646  nnz++;
3647  }
3648  }
3649
3650  mfem::out << "Number of nonzeros in NNLSSolver solution: " << nnz
3651  << ", out of " << sol.Size() << endl;
3652
3653  // Check residual of NNLS solution
3654  Vector res(mat->NumRows());
3655  mat->Mult(sol, res);
3656  res *= row_scaling_;
3657
3658  const double normGsol = res.Norml2();
3659  const double normRHS = rhs_Gw.Norml2();
3660
3661  res -= rhs_Gw;
3662  const double relNorm = res.Norml2() / std::max(normGsol, normRHS);
3663  mfem::out << "Relative residual norm for NNLSSolver solution of Gs = Gw: "
3664  << relNorm << endl;
3665  }
3666 }
3667
3668 void NNLSSolver::Solve(const Vector& rhs_lb, const Vector& rhs_ub,
3669  Vector& soln) const
3670 {
3671  int m = mat->NumRows();
3672  int n = mat->NumCols();
3673
3674  MFEM_VERIFY(rhs_lb.Size() == m && rhs_lb.Size() == m && soln.Size() == n, "");
3675  MFEM_VERIFY(n >= m, "NNLSSolver system cannot be over-determined.");
3676
3677  if (max_nnz_ == 0)
3678  {
3679  max_nnz_ = mat->NumCols();
3680  }
3681
3682  // Prepare right hand side
3683  Vector rhs_avg(rhs_ub);
3684  rhs_avg += rhs_lb;
3685  rhs_avg *= 0.5;
3686
3687  Vector rhs_halfgap(rhs_ub);
3688  rhs_halfgap -= rhs_lb;
3689  rhs_halfgap *= 0.5;
3690
3691  Vector rhs_avg_glob(rhs_avg);
3692  Vector rhs_halfgap_glob(rhs_halfgap);
3693
3694  int ione = 1;
3695  double fone = 1.0;
3696
3697  char lside = 'L';
3698  char trans = 'T';
3699  char notrans = 'N';
3700
3701  std::vector<unsigned int> nz_ind(m);
3702  Vector res_glob(m);
3703  Vector mu(n);
3704  Vector mu2(n);
3705  int n_nz_ind = 0;
3706  int n_glob = 0;
3707  int m_update;
3708  int min_nnz_cap = std::min(static_cast<int>(min_nnz_), std::min(m,n));
3709  int info;
3710  std::vector<double> l2_res_hist;
3711  std::vector<unsigned int> stalled_indices;
3712  int stalledFlag = 0;
3713  int num_stalled = 0;
3714  int nz_ind_zero = 0;
3715
3716  Vector soln_nz_glob(m);
3717  Vector soln_nz_glob_up(m);
3718
3719  // The following matrices are stored in column-major format as Vectors
3720  Vector mat_0_data(m * n);
3721  Vector mat_qr_data(m * n);
3722  Vector submat_data(m * n);
3723
3724  Vector tau(n);
3725  Vector sub_tau = tau;
3726  Vector vec1(m);
3727
3728  // Temporary work arrays
3729  int lwork;
3730  std::vector<double> work;
3731  int n_outer_iter = 0;
3732  int n_total_inner_iter = 0;
3733  int i_qr_start;
3734  int n_update;
3735  // 0 = converged; 1 = maximum iterations reached;
3736  // 2 = NNLS stalled (no change in residual for many iterations)
3737  int exit_flag = 1;
3738
3739  res_glob = rhs_avg_glob;
3740  Vector qt_rhs_glob = rhs_avg_glob;
3741  Vector qqt_rhs_glob = qt_rhs_glob;
3742  Vector sub_qt = rhs_avg_glob;
3743
3744  // Compute threshold tolerance for the Lagrange multiplier mu
3745  double mu_tol = 0.0;
3746
3747  {
3748  Vector rhs_scaled(rhs_halfgap_glob);
3749  Vector tmp(n);
3750  rhs_scaled *= row_scaling_;
3751  mat->MultTranspose(rhs_scaled, tmp);
3752
3753  mu_tol = 1.0e-15 * tmp.Max();
3754  }
3755
3756  double rmax = 0.0;
3757  double mumax = 0.0;
3758
3759  for (int oiter = 0; oiter < n_outer_; ++oiter)
3760  {
3761  stalledFlag = 0;
3762
3763  rmax = fabs(res_glob(0)) - rhs_halfgap_glob(0);
3764  for (int i=1; i<m; ++i)
3765  {
3766  rmax = std::max(rmax, fabs(res_glob(i)) - rhs_halfgap_glob(i));
3767  }
3768
3769  l2_res_hist.push_back(res_glob.Norml2());
3770
3771  if (verbosity_ > 1)
3772  {
3773  mfem::out << "NNLS " << oiter << " " << n_total_inner_iter << " " << m
3774  << " " << n << " " << n_glob << " " << rmax << " "
3775  << l2_res_hist[oiter] << endl;
3776  }
3777  if (rmax <= const_tol_ && n_glob >= min_nnz_cap)
3778  {
3779  if (verbosity_ > 1)
3780  {
3781  mfem::out << "NNLS target tolerance met" << endl;
3782  }
3783  exit_flag = 0;
3784  break;
3785  }
3786
3787  if (n_glob >= max_nnz_)
3788  {
3789  if (verbosity_ > 1)
3790  {
3791  mfem::out << "NNLS target nnz met" << endl;
3792  }
3793  exit_flag = 0;
3794  break;
3795  }
3796
3797  if (n_glob >= m)
3798  {
3799  if (verbosity_ > 1)
3800  {
3801  mfem::out << "NNLS system is square... exiting" << endl;
3802  }
3803  exit_flag = 3;
3804  break;
3805  }
3806
3807  // Check for stall after the first nStallCheck iterations
3808  if (oiter > nStallCheck_)
3809  {
3810  double mean0 = 0.0;
3811  double mean1 = 0.0;
3812  for (int i=0; i<nStallCheck_/2; ++i)
3813  {
3814  mean0 += l2_res_hist[oiter - i];
3815  mean1 += l2_res_hist[oiter - (nStallCheck_) - i];
3816  }
3817
3818  double mean_res_change = (mean1 / mean0) - 1.0;
3819  if (std::abs(mean_res_change) < res_change_termination_tol_)
3820  {
3821  if (verbosity_ > 1)
3822  {
3823  mfem::out << "NNLSSolver stall detected... exiting" << endl;
3824  }
3825  exit_flag = 2;
3826  break;
3827  }
3828  }
3829
3830  // Find the next index
3831  res_glob *= row_scaling_;
3832  mat->MultTranspose(res_glob, mu);
3833
3834  for (int i = 0; i < n_nz_ind; ++i)
3835  {
3836  mu(nz_ind[i]) = 0.0;
3837  }
3838  for (unsigned int i = 0; i < stalled_indices.size(); ++i)
3839  {
3840  mu(stalled_indices[i]) = 0.0;
3841  }
3842
3843  mumax = mu.Max();
3844
3845  if (mumax < mu_tol)
3846  {
3847  num_stalled = stalled_indices.size();
3848  if (num_stalled > 0)
3849  {
3850  if (verbosity_ > 0)
3851  {
3852  mfem::out << "NNLS Lagrange multiplier is below the minimum "
3853  << "threshold: mumax = " << mumax << ", mutol = "
3854  << mu_tol << "\n" << " Resetting stalled indices "
3855  << "vector of size " << num_stalled << "\n";
3856  }
3857  stalled_indices.resize(0);
3858
3859  mat->MultTranspose(res_glob, mu);
3860
3861  for (int i = 0; i < n_nz_ind; ++i)
3862  {
3863  mu(nz_ind[i]) = 0.0;
3864  }
3865
3866  mumax = mu.Max();
3867  }
3868  }
3869
3870  int imax = 0;
3871  {
3872  double tmax = mu(0);
3873  for (int i=1; i<n; ++i)
3874  {
3875  if (mu(i) > tmax)
3876  {
3877  tmax = mu(i);
3878  imax = i;
3879  }
3880  }
3881  }
3882
3883  // Record the local value of the next index
3884  nz_ind[n_nz_ind] = imax;
3885  ++n_nz_ind;
3886
3887  if (verbosity_ > 2)
3888  {
3889  mfem::out << "Found next index: " << imax << " " << mumax << endl;
3890  }
3891
3892  for (int i=0; i<m; ++i)
3893  {
3894  mat_0_data(i + (n_glob*m)) = (*mat)(i,imax) * row_scaling_[i];
3895  mat_qr_data(i + (n_glob*m)) = mat_0_data(i + (n_glob*m));
3896  }
3897
3898  i_qr_start = n_glob;
3899  ++n_glob; // Increment the size of the global matrix
3900
3901  if (verbosity_ > 2)
3902  {
3903  mfem::out << "Updated matrix with new index" << endl;
3904  }
3905
3906  for (int iiter = 0; iiter < n_inner_; ++iiter)
3907  {
3908  ++n_total_inner_iter;
3909
3910  // Initialize
3911  const bool incremental_update = true;
3912  n_update = n_glob - i_qr_start;
3913  m_update = m - i_qr_start;
3914  if (incremental_update)
3915  {
3916  // Apply Householder reflectors to compute Q^T new_cols
3917  lwork = -1;
3918  work.resize(10);
3919
3920  dormqr_(&lside, &trans, &m, &n_update, &i_qr_start,
3921  mat_qr_data.GetData(), &m, tau.GetData(),
3922  mat_qr_data.GetData() + (i_qr_start * m), &m,
3923  work.data(), &lwork, &info);
3924  MFEM_VERIFY(info == 0, ""); // Q^T A update work calculation failed
3925  lwork = static_cast<int>(work[0]);
3926  work.resize(lwork);
3927  dormqr_(&lside, &trans, &m, &n_update, &i_qr_start,
3928  mat_qr_data.GetData(), &m, tau.GetData(),
3929  mat_qr_data.GetData() + (i_qr_start * m), &m,
3930  work.data(), &lwork, &info);
3931  MFEM_VERIFY(info == 0, ""); // Q^T A update failed
3932  // Compute QR factorization of the submatrix
3933  lwork = -1;
3934  work.resize(10);
3935
3936  // Copy m_update-by-n_update submatrix of mat_qr_data,
3937  // starting at (i_qr_start, i_qr_start)
3938  for (int i=0; i<m_update; ++i)
3939  for (int j=0; j<n_update; ++j)
3940  {
3941  submat_data[i + (j * m_update)] =
3942  mat_qr_data[i + i_qr_start + ((j + i_qr_start) * m)];
3943  }
3944
3945  // Copy tau subvector of length n_update, starting at i_qr_start
3946  for (int j=0; j<n_update; ++j)
3947  {
3948  sub_tau[j] = tau[i_qr_start + j];
3949  }
3950
3951  dgeqrf_(&m_update, &n_update,
3952  submat_data.GetData(), &m_update, sub_tau.GetData(),
3953  work.data(), &lwork, &info);
3954  MFEM_VERIFY(info == 0, ""); // QR update factorization work calc
3955  lwork = static_cast<int>(work[0]);
3956  if (lwork == 0) { lwork = 1; }
3957  work.resize(lwork);
3958  dgeqrf_(&m_update, &n_update,
3959  submat_data.GetData(), &m_update, sub_tau.GetData(),
3960  work.data(), &lwork, &info);
3961  MFEM_VERIFY(info == 0, ""); // QR update factorization failed
3962
3963  // Copy result back
3964  for (int i=0; i<m_update; ++i)
3965  for (int j=0; j<n_update; ++j)
3966  {
3967  mat_qr_data[i + i_qr_start + ((j + i_qr_start)* m)] =
3968  submat_data[i + (j * m_update)];
3969  }
3970
3971  for (int j=0; j<n_update; ++j)
3972  {
3973  tau[i_qr_start + j] = sub_tau[j];
3974  }
3975  }
3976  else
3977  {
3978  // Copy everything to mat_qr then do full QR
3979  for (int i=0; i<m; ++i)
3980  for (int j=0; j<n_glob; ++j)
3981  {
3982  mat_qr_data(i + (j*m)) = mat_0_data(i + (j*m));
3983  }
3984
3985  // Compute qr factorization (first find the size of work and then
3986  // perform qr)
3987  lwork = -1;
3988  work.resize(10);
3989  dgeqrf_(&m, &n_glob,
3990  mat_qr_data.GetData(), &m, tau.GetData(),
3991  work.data(), &lwork, &info);
3992  MFEM_VERIFY(info == 0, ""); // QR factorization work calculation
3993  lwork = static_cast<int>(work[0]);
3994  work.resize(lwork);
3995  dgeqrf_(&m, &n_glob,
3996  mat_qr_data.GetData(), &m, tau.GetData(),
3997  work.data(), &lwork, &info);
3998  MFEM_VERIFY(info == 0, ""); // QR factorization failed
3999  }
4000
4001  if (verbosity_ > 2)
4002  {
4003  mfem::out << "Updated QR " << iiter << endl;
4004  }
4005
4006  // Apply Householder reflectors to compute Q^T b
4007  if (incremental_update && iiter == 0)
4008  {
4009  lwork = -1;
4010  work.resize(10);
4011
4012  // Copy submatrix of mat_qr_data starting at
4013  // (i_qr_start, i_qr_start), of size m_update-by-1
4014  // Copy submatrix of qt_rhs_glob starting at (i_qr_start, 0),
4015  // of size m_update-by-1
4016
4017  for (int i=0; i<m_update; ++i)
4018  {
4019  submat_data[i] = mat_qr_data[i + i_qr_start + (i_qr_start * m)];
4020  sub_qt[i] = qt_rhs_glob[i + i_qr_start];
4021  }
4022
4023  sub_tau[0] = tau[i_qr_start];
4024
4025  dormqr_(&lside, &trans, &m_update, &ione, &ione,
4026  submat_data.GetData(), &m_update, sub_tau.GetData(),
4027  sub_qt.GetData(), &m_update,
4028  work.data(), &lwork, &info);
4029  MFEM_VERIFY(info == 0, ""); // H_last y work calculation failed
4030  lwork = static_cast<int>(work[0]);
4031  work.resize(lwork);
4032  dormqr_(&lside, &trans, &m_update, &ione, &ione,
4033  submat_data.GetData(), &m_update, sub_tau.GetData(),
4034  sub_qt.GetData(), &m_update,
4035  work.data(), &lwork, &info);
4036  MFEM_VERIFY(info == 0, ""); // H_last y failed
4037  // Copy result back
4038  for (int i=0; i<m_update; ++i)
4039  {
4040  qt_rhs_glob[i + i_qr_start] = sub_qt[i];
4041  }
4042  }
4043  else
4044  {
4045  // Compute Q^T b from scratch
4046  qt_rhs_glob = rhs_avg_glob;
4047  lwork = -1;
4048  work.resize(10);
4049  dormqr_(&lside, &trans, &m, &ione, &n_glob,
4050  mat_qr_data.GetData(), &m, tau.GetData(),
4051  qt_rhs_glob.GetData(), &m,
4052  work.data(), &lwork, &info);
4053  MFEM_VERIFY(info == 0, ""); // Q^T b work calculation failed
4054  lwork = static_cast<int>(work[0]);
4055  work.resize(lwork);
4056  dormqr_(&lside, &trans, &m, &ione, &n_glob,
4057  mat_qr_data.GetData(), &m, tau.GetData(),
4058  qt_rhs_glob.GetData(), &m,
4059  work.data(), &lwork, &info);
4060  MFEM_VERIFY(info == 0, ""); // Q^T b failed
4061  }
4062
4063  if (verbosity_ > 2)
4064  {
4065  mfem::out << "Updated rhs " << iiter << endl;
4066  }
4067
4068  // Apply R^{-1}; first n_glob entries of vec1 are overwritten
4069  char upper = 'U';
4070  char nounit = 'N';
4071  vec1 = qt_rhs_glob;
4072  dtrsm_(&lside, &upper, &notrans, &nounit,
4073  &n_glob, &ione, &fone,
4074  mat_qr_data.GetData(), &m,
4075  vec1.GetData(), &n_glob);
4076
4077  if (verbosity_ > 2)
4078  {
4079  mfem::out << "Solved triangular system " << iiter << endl;
4080  }
4081
4082  // Check if all entries are positive
4083  int pos_ibool = 0;
4084  double smin = n_glob > 0 ? vec1(0) : 0.0;
4085  for (int i=0; i<n_glob; ++i)
4086  {
4087  soln_nz_glob_up(i) = vec1(i);
4088  smin = std::min(smin, soln_nz_glob_up(i));
4089  }
4090
4091  if (smin > zero_tol_)
4092  {
4093  pos_ibool = 1;
4094  for (int i=0; i<n_glob; ++i)
4095  {
4096  soln_nz_glob(i) = soln_nz_glob_up(i);
4097  }
4098  }
4099
4100  if (pos_ibool == 1)
4101  {
4102  break;
4103  }
4104
4105  if (verbosity_ > 2)
4106  {
4107  mfem::out << "Start pruning " << iiter << endl;
4108  for (int i = 0; i < n_glob; ++i)
4109  {
4110  if (soln_nz_glob_up(i) <= zero_tol_)
4111  {
4112  mfem::out << i << " " << n_glob << " " << soln_nz_glob_up(i) << endl;
4113  }
4114  }
4115  }
4116
4117  if (soln_nz_glob_up(n_glob - 1) <= zero_tol_)
4118  {
4119  stalledFlag = 1;
4120  if (verbosity_ > 2)
4121  {
4122  if (qr_residual_mode_ == QRresidualMode::hybrid)
4123  {
4124  mfem::out << "Detected stall due to adding and removing same "
4125  << "column. Switching to QR residual calculation "
4126  << "method." << endl;
4127  }
4128  else
4129  {
4130  mfem::out << "Detected stall due to adding and removing same"
4131  << " column. Exiting now." << endl;
4132  }
4133  }
4134  }
4135
4136  if (stalledFlag == 1 && qr_residual_mode_ == QRresidualMode::hybrid)
4137  {
4138  NNLS_qrres_on_ = true;
4139  break;
4140  }
4141
4142  double alpha = 1.0e300;
4143  // Find maximum permissible step
4144  for (int i = 0; i < n_glob; ++i)
4145  {
4146  if (soln_nz_glob_up(i) <= zero_tol_)
4147  {
4148  alpha = std::min(alpha, soln_nz_glob(i)/(soln_nz_glob(i) - soln_nz_glob_up(i)));
4149  }
4150  }
4151  // Update solution
4152  smin = 0.0;
4153  for (int i = 0; i < n_glob; ++i)
4154  {
4155  soln_nz_glob(i) += alpha*(soln_nz_glob_up(i) - soln_nz_glob(i));
4156  if (i == 0 || soln_nz_glob(i) < smin)
4157  {
4158  smin = soln_nz_glob(i);
4159  }
4160  }
4161
4162  while (smin > zero_tol_)
4163  {
4164  // This means there was a rounding error, as we should have
4165  // a zero element by definition. Recalculate alpha based on
4166  // the index that corresponds to the element that should be
4167  // zero.
4168
4169  int index_min = 0;
4170  smin = soln_nz_glob(0);
4171  for (int i = 1; i < n_glob; ++i)
4172  {
4173  if (soln_nz_glob(i) < smin)
4174  {
4175  smin = soln_nz_glob(i);
4176  index_min = i;
4177  }
4178  }
4179
4180  alpha = soln_nz_glob(index_min)/(soln_nz_glob(index_min)
4181  - soln_nz_glob_up(index_min));
4182
4183  // Reupdate solution
4184  for (int i = 0; i < n_glob; ++i)
4185  {
4186  soln_nz_glob(i) += alpha*(soln_nz_glob_up(i) - soln_nz_glob(i));
4187  }
4188  }
4189
4190  // Clean up zeroed entry
4191  i_qr_start = n_glob+1;
4192  while (true)
4193  {
4194  // Check if there is a zero entry
4195  int zero_ibool;
4196
4197  smin = n_glob > 0 ? soln_nz_glob(0) : 0.0;
4198  for (int i=1; i<n_glob; ++i)
4199  {
4200  smin = std::min(smin, soln_nz_glob(i));
4201  }
4202
4203  if (smin < zero_tol_)
4204  {
4205  zero_ibool = 1;
4206  }
4207  else
4208  {
4209  zero_ibool = 0;
4210  }
4211
4212  if (zero_ibool == 0) // Break if there is no more zero entry
4213  {
4214  break;
4215  }
4216
4217  int ind_zero = -1; // Index where the first zero is encountered
4218  nz_ind_zero = 0;
4219
4220  // Identify global index of the zeroed element
4221  for (int i = 0; i < n_glob; ++i)
4222  {
4223  if (soln_nz_glob(i) < zero_tol_)
4224  {
4225  ind_zero = i;
4226  break;
4227  }
4228  }
4229  MFEM_VERIFY(ind_zero != -1, "");
4230  // Identify the local index for nz_ind to which the zeroed entry
4231  // belongs
4232  for (int i = 0; i < ind_zero; ++i)
4233  {
4234  ++nz_ind_zero;
4235  }
4236
4237  {
4238  // Copy mat_0.cols[ind_zero+1,n_glob) to mat_qr.cols[ind_zero,n_glob-1)
4239  for (int i=0; i<m; ++i)
4240  for (int j=ind_zero; j<n_glob-1; ++j)
4241  {
4242  mat_qr_data(i + (j*m)) = mat_0_data(i + ((j+1)*m));
4243  }
4244
4245  // Copy mat_qr.cols[ind_zero,n_glob-1) to
4246  // mat_0.cols[ind_zero,n_glob-1)
4247  for (int i=0; i<m; ++i)
4248  for (int j=ind_zero; j<n_glob-1; ++j)
4249  {
4250  mat_0_data(i + (j*m)) = mat_qr_data(i + (j*m));
4251  }
4252  }
4253
4254  // Remove the zeroed entry from the local matrix index
4255  for (int i = nz_ind_zero; i < n_nz_ind-1; ++i)
4256  {
4257  nz_ind[i] = nz_ind[i+1];
4258  }
4259  --n_nz_ind;
4260
4261  // Shift soln_nz_glob and proc_index
4262  for (int i = ind_zero; i < n_glob-1; ++i)
4263  {
4264  soln_nz_glob(i) = soln_nz_glob(i+1);
4265  }
4266
4267  i_qr_start = std::min(i_qr_start, ind_zero);
4268  --n_glob;
4269  } // End of pruning loop
4270
4271  if (verbosity_ > 2)
4272  {
4273  mfem::out << "Finished pruning " << iiter << endl;
4274  }
4275  } // End of inner loop
4276
4277  // Check if we have stalled
4278  if (stalledFlag == 1)
4279  {
4280  --n_glob;
4281  --n_nz_ind;
4282  num_stalled = stalled_indices.size();
4283  stalled_indices.resize(num_stalled + 1);
4284  stalled_indices[num_stalled] = imax;
4285  if (verbosity_ > 2)
4286  {
4287  mfem::out << "Adding index " << imax << " to stalled index list "
4288  << "of size " << num_stalled << endl;
4289  }
4290  }
4291
4292  // Compute residual
4293  if (!NNLS_qrres_on_)
4294  {
4295  res_glob = rhs_avg_glob;
4296  double fmone = -1.0;
4297  dgemv_(&notrans, &m, &n_glob, &fmone,
4298  mat_0_data.GetData(), &m,
4299  soln_nz_glob.GetData(), &ione, &fone,
4300  res_glob.GetData(), &ione);
4301  }
4302  else
4303  {
4304  // Compute residual using res = b - Q*Q^T*b, where Q is from an
4305  // economical QR decomposition
4306  lwork = -1;
4307  work.resize(10);
4308  qqt_rhs_glob = 0.0;
4309  for (int i=0; i<n_glob; ++i)
4310  {
4311  qqt_rhs_glob(i) = qt_rhs_glob(i);
4312  }
4313
4314  dormqr_(&lside, &notrans, &m, &ione, &n_glob, mat_qr_data.GetData(), &m,
4315  tau.GetData(), qqt_rhs_glob.GetData(), &m,
4316  work.data(), &lwork, &info);
4317
4318  MFEM_VERIFY(info == 0, ""); // Q Q^T b work calculation failed.
4319  lwork = static_cast<int>(work[0]);
4320  work.resize(lwork);
4321  dormqr_(&lside, &notrans, &m, &ione, &n_glob, mat_qr_data.GetData(), &m,
4322  tau.GetData(), qqt_rhs_glob.GetData(), &m,
4323  work.data(), &lwork, &info);
4324  MFEM_VERIFY(info == 0, ""); // Q Q^T b calculation failed.
4325  res_glob = rhs_avg_glob;
4326  res_glob -= qqt_rhs_glob;
4327  }
4328
4329  if (verbosity_ > 2)
4330  {
4331  mfem::out << "Computed residual" << endl;
4332  }
4333
4334  ++n_outer_iter;
4335  } // End of outer loop
4336
4337  // Insert the solutions
4338  MFEM_VERIFY(n_glob == n_nz_ind, "");
4339  soln = 0.0;
4340  for (int i = 0; i < n_glob; ++i)
4341  {
4342  soln(nz_ind[i]) = soln_nz_glob(i);
4343  }
4344
4345  if (verbosity_ > 0)
4346  {
4347  mfem::out << "NNLS solver: m = " << m << ", n = " << n
4348  << ", outer_iter = " << n_outer_iter << ", inner_iter = "
4349  << n_total_inner_iter;
4350
4351  if (exit_flag == 0)
4352  {
4353  mfem::out << ": converged" << endl;
4354  }
4355  else
4356  {
4357  mfem::out << endl << "Warning, NNLS convergence stalled: "
4358  << (exit_flag == 2) << endl;
4359  mfem::out << "resErr = " << rmax << " vs tol = " << const_tol_
4360  << "; mumax = " << mumax << " vs tol = " << mu_tol << endl;
4361  }
4362  }
4363 }
4364 #endif // MFEM_USE_LAPACK
4365
4366 }
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
const Array< int > * ess_dofs_list
Not owned.
Definition: solvers.hpp:1056
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:952
PowerMethod helper class to estimate the largest eigenvalue of an operator using the iterative power ...
Definition: operator.hpp:986
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 double * HostReadData() const
Definition: sparsemat.hpp:277
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
void dgemv_(char *, int *, int *, double *, double *, int *, double *, int *, double *, double *, int *)
Conjugate gradient method.
Definition: solvers.hpp:493
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
double Info[UMFPACK_INFO]
Definition: solvers.hpp:1082
Chebyshev accelerated smoothing with given vector, no matrix necessary.
Definition: solvers.hpp:379
SparseMatrix * mat
Definition: solvers.hpp:1074
const Operator * oper
Definition: solvers.hpp:120
void dgeqrf_(int *, int *, double *, int *, double *, double *, int *, int *)
Array< Vector * > ykArray
Definition: solvers.hpp:730
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:122
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:669
SuiteSparse_long * AI
Definition: solvers.hpp:1076
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void Mult(const Vector &w, Vector &sol) const override
Operator application: y=A(x).
Definition: solvers.cpp:3620
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:457
void MinimumDiscardedFillOrdering(SparseMatrix &C, Array< int > &p)
Definition: solvers.cpp:2640
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:391
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:344
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Definition: solvers.hpp:91
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:303
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:704
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:2007
bool first_and_last
Information about the first and last iteration will be printed to mfem::out.
Definition: solvers.hpp:94
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1168
OperatorJacobiSmoother(const double damping=1.0)
Default constructor: the diagonal will be computed by subsequent calls to SetOperator() using the Ope...
Definition: solvers.cpp:200
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
void dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m, int *n, double *alpha, double *a, int *lda, double *b, int *ldb)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1614
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
virtual void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3400
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:980
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:465
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:86
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:531
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
bool errors
If a fatal problem has been detected the failure will be reported to mfem::err.
Definition: solvers.hpp:82
void SetQRResidualMode(const QRresidualMode qr_residual_mode)
Definition: solvers.cpp:3579
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 &x, Vector &y) const =0
Operator application: y=A(x).
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3319
Vector beta
void SetOperator(const Operator &op)
Definition: solvers.cpp:2765
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
STL namespace.
MINRES method.
Definition: solvers.hpp:603
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
Reordering
The reordering method used by the BlockILU factorization.
Definition: solvers.hpp:968
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2886
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:2996
int print_level
(DEPRECATED) Legacy print level definition, which is left for compatibility with custom iterative sol...
Definition: solvers.hpp:130
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:457
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
void Monitor(int it, double norm, const Vector &r, const Vector &x, bool final=false) const
Definition: solvers.cpp:190
const Vector * d_hi
Definition: solvers.hpp:829
double GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult()...
Definition: solvers.hpp:265
virtual void SetOperator(const Operator &op)
Set the Operator that is the OrthoSolver is to invert (approximately).
Definition: solvers.cpp:3442
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3372
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:1570
bool IsFinite(const double &val)
Definition: vector.hpp:486
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
const Vector * x_lo
Definition: solvers.hpp:829
double Norm(const Vector &x) const
Definition: solvers.hpp:172
Operator * grad
Definition: solvers.hpp:645
const Vector * x_hi
Definition: solvers.hpp:829
const class IterativeSolver * iter_solver
The last IterativeSolver to which this monitor was attached.
Definition: solvers.hpp:40
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2320
General product operator: x -> (A*B)(x) = A(B(x)).
Definition: operator.hpp:774
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:616
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
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:3333
void Solve(const Vector &rhs_lb, const Vector &rhs_ub, Vector &soln) const
Solve the NNLS problem. Specifically, we find a vector soln, such that rhs_lb < mat*soln < rhs_ub is ...
Definition: solvers.cpp:3668
int GetNumConstraints() const
Definition: solvers.cpp:2340
Data type sparse matrix.
Definition: sparsemat.hpp:50
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:904
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:1780
virtual ~UMFPackSolver()
Definition: solvers.cpp:3266
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:580
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:214
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition: solvers.hpp:250
double b
Definition: lissajous.cpp:42
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector *> &v)
Definition: solvers.cpp:959
const int * HostReadJ() const
Definition: sparsemat.hpp:261
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
double delta
Definition: lissajous.cpp:43
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d, const Array< int > &ess_tdof_list, int order, double max_eig_estimate)
Definition: solvers.cpp:331
void SetLinearConstraint(const Vector &w_, double a_)
Definition: solvers.cpp:2367
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:1286
DirectSubBlockSolver(const SparseMatrix &A, const SparseMatrix &block_dof)
Definition: solvers.cpp:3357
SuiteSparse_long * AJ
Definition: solvers.hpp:1076
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:538
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:2373
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:298
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:898
double * GetData(int k)
Definition: densemat.hpp:1201
void Setup(const Vector &diag)
Definition: solvers.cpp:277
Parallel smoothers in hypre.
Definition: hypre.hpp:971
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:913
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:1148
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:3288
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1593
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1481
Array< Vector * > skArray
Definition: solvers.hpp:730
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3232
prob_type prob
Definition: ex25.cpp:154
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:172
klu_symbolic * Symbolic
Definition: solvers.hpp:1115
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition: array.hpp:311
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:671
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:552
BlockILU(int block_size_, Reordering reordering_=Reordering::MINIMUM_DISCARDED_FILL, int k_fill_=0)
Definition: solvers.cpp:2747
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
Stationary linear iteration: x <- x + B (b - A x)
Definition: solvers.hpp:461
double omega
Definition: ex25.cpp:142
void SetInequalityConstraint(const Vector &dl, const Vector &dh)
Definition: solvers.cpp:2322
void SetEqualityConstraint(const Vector &c)
Definition: solvers.cpp:2314
void dormqr_(char *, char *, int *, int *, int *, double *, int *, double *, double *, int *, double *, int *, int *)
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:908
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1819
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:638
void SetAbsTol(double atol)
Definition: solvers.hpp:200
PrintLevel print_options
Output behavior for the iterative solver.
Definition: solvers.hpp:137
void SetRelTol(double rtol)
Definition: solvers.hpp:199
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1081
Abstract base class for iterative solver.
Definition: solvers.hpp:66
void forall(int N, lambda &&body)
Definition: forall.hpp:742
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:3229
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
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:473
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:1927
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
Definition: solvers.hpp:252
int max_iter
Limit for the number of iterations the solver is allowed to do.
Definition: solvers.hpp:151
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3194
GMRES method.
Definition: solvers.hpp:525
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:1987
static int GuessLegacyPrintLevel(PrintLevel)
Use some heuristics to guess a legacy print level corresponding to the given PrintLevel.
Definition: solvers.cpp:149
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:154
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:925
const int * HostReadI() const
Definition: sparsemat.hpp:245
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:263
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:2348
void NormalizeConstraints(Vector &rhs_lb, Vector &rhs_ub) const
Definition: solvers.cpp:3588
double a
Definition: lissajous.cpp:41
void UpdateVectors()
Definition: solvers.cpp:525
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:414
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:469
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
klu_common Common
Definition: solvers.hpp:1136
double rel_tol
Relative tolerance.
Definition: solvers.hpp:154
OptimizationProblem(int insize, const Operator *C_, const Operator *D_)
In parallel, insize is the number of the local true dofs.
Definition: solvers.cpp:2304
void SetSolutionBounds(const Vector &xl, const Vector &xh)
Definition: solvers.cpp:2332
virtual void Mult(const Vector &xt, Vector &x) const
Definition: solvers.cpp:2382
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map, bool op_is_symmetric=true, bool own_aux_map=false)
Definition: solvers.cpp:3509
double mu
Definition: ex25.cpp:140
void SetOperator(const Operator &op)
Recompute the diagonal using the method AssembleDiagonal of the given new Operator, op.
Definition: solvers.cpp:241
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:930
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:60
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition: operator.hpp:131
PrintLevel FromLegacyPrintLevel(int)
Convert a legacy print level integer to a PrintLevel object.
Definition: solvers.cpp:114
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:586
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:835
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:55
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:292
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:58
void SetSolver(Solver &s)
Set the solver used by the OrthoSolver.
Definition: solvers.cpp:3433
void Mult(const Vector &b, Vector &x) const
Perform the action of the OrthoSolver: P * solver * P where P is the projection to the subspace of ve...
Definition: solvers.cpp:3452
void MonitorResidual(int it, double norm, const Vector &r, bool final) override
Monitor the residual vector r.
Definition: solvers.cpp:3047
const Vector * d_lo
Definition: solvers.hpp:829
virtual void MonitorSolution(int it, double norm, const Vector &x, bool final)
Monitor the solution vector x.
Definition: solvers.hpp:54
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
int aGMRES(const Operator &A, Vector &x, const Vector &b, const Operator &M, int &max_iter, int m_max, int m_min, int m_step, double cf, double &tol, double &atol, int printit)
Definition: solvers.cpp:2149
virtual ~KLUSolver()
Definition: solvers.cpp:3347
RefCoord s[3]
void UpdateVectors()
Definition: solvers.cpp:709
BiCGSTAB method.
Definition: solvers.hpp:572
bool warnings
If a non-fatal problem has been detected some context-specific information will be reported to mfem::...
Definition: solvers.hpp:85
bool iterations
Detailed information about each iteration will be reported to mfem::out.
Definition: solvers.hpp:88
void SetOperator(const Operator &op) override
The operator must be a DenseMatrix.
Definition: solvers.cpp:3566
SparseMatrix * mat
Definition: solvers.hpp:1114
Base class for solvers.
Definition: operator.hpp:682
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
const OptimizationProblem * problem
Definition: solvers.hpp:868
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1378
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:710
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:699
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition: solvers.cpp:2361
void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
Definition: densemat.cpp:2413
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3387
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:3404
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:1940
const Operator * D
Definition: solvers.hpp:828
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:485
const Vector * c_e
Definition: solvers.hpp:829
IterativeSolverMonitor * monitor
Definition: solvers.hpp:122
klu_numeric * Numeric
Definition: solvers.hpp:1116
virtual void MonitorResidual(int it, double norm, const Vector &r, bool final)
Monitor the residual vector r.
Definition: solvers.hpp:48
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:473
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:475
Settings for the output behavior of the IterativeSolver.
Definition: solvers.hpp:78
const Operator * C
Not owned, some can remain unused (NULL).
Definition: solvers.hpp:828
double abs_tol
Absolute tolerance.
Definition: solvers.hpp:157
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
void SLI(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Stationary linear iteration. (tolerances are squared)
Definition: solvers.cpp:677
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:1340
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1808
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3099