MFEM  v4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
solvers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "linalg.hpp"
13 #include "../general/globals.hpp"
14 #include <iostream>
15 #include <iomanip>
16 #include <algorithm>
17 #include <cmath>
18 
19 namespace mfem
20 {
21 
22 using namespace std;
23 
25  : Solver(0, true)
26 {
27  oper = NULL;
28  prec = NULL;
29  max_iter = 10;
30  print_level = -1;
31  rel_tol = abs_tol = 0.0;
32 #ifdef MFEM_USE_MPI
33  dot_prod_type = 0;
34 #endif
35 }
36 
37 #ifdef MFEM_USE_MPI
39  : Solver(0, true)
40 {
41  oper = NULL;
42  prec = NULL;
43  max_iter = 10;
44  print_level = -1;
45  rel_tol = abs_tol = 0.0;
46  dot_prod_type = 1;
47  comm = _comm;
48 }
49 #endif
50 
51 double IterativeSolver::Dot(const Vector &x, const Vector &y) const
52 {
53 #ifndef MFEM_USE_MPI
54  return (x * y);
55 #else
56  if (dot_prod_type == 0)
57  {
58  return (x * y);
59  }
60  else
61  {
62  return InnerProduct(comm, x, y);
63  }
64 #endif
65 }
66 
67 void IterativeSolver::SetPrintLevel(int print_lvl)
68 {
69 #ifndef MFEM_USE_MPI
70  print_level = print_lvl;
71 #else
72  if (dot_prod_type == 0)
73  {
74  print_level = print_lvl;
75  }
76  else
77  {
78  int rank;
79  MPI_Comm_rank(comm, &rank);
80  if (rank == 0)
81  {
82  print_level = print_lvl;
83  }
84  }
85 #endif
86 }
87 
89 {
90  prec = &pr;
91  prec->iterative_mode = false;
92 }
93 
95 {
96  oper = &op;
97  height = op.Height();
98  width = op.Width();
99  if (prec)
100  {
101  prec->SetOperator(*oper);
102  }
103 }
104 
105 
107 {
108  r.SetSize(width);
109  z.SetSize(width);
110 }
111 
112 void SLISolver::Mult(const Vector &b, Vector &x) const
113 {
114  int i;
115 
116  // Optimized preconditioned SLI with fixed number of iterations and given
117  // initial guess
118  if (!rel_tol && iterative_mode && prec)
119  {
120  for (i = 0; i < max_iter; i++)
121  {
122  oper->Mult(x, r); // r = A x
123  subtract(b, r, r); // r = b - A x
124  prec->Mult(r, z); // z = B r
125  add(x, 1.0, z, x); // x = x + B (b - A x)
126  }
127  converged = 1;
128  final_iter = i;
129  return;
130  }
131 
132  // Optimized preconditioned SLI with fixed number of iterations and zero
133  // initial guess
134  if (!rel_tol && !iterative_mode && prec)
135  {
136  prec->Mult(b, x); // x = B b (initial guess 0)
137  for (i = 1; i < max_iter; i++)
138  {
139  oper->Mult(x, r); // r = A x
140  subtract(b, r, r); // r = b - A x
141  prec->Mult(r, z); // z = B r
142  add(x, 1.0, z, x); // x = x + B (b - A x)
143  }
144  converged = 1;
145  final_iter = i;
146  return;
147  }
148 
149  // General version of SLI with a relative tolerance, optional preconditioner
150  // and optional initial guess
151  double r0, nom, nom0, nomold = 1, cf;
152 
153  if (iterative_mode)
154  {
155  oper->Mult(x, r);
156  subtract(b, r, r); // r = b - A x
157  }
158  else
159  {
160  r = b;
161  x = 0.0;
162  }
163 
164  if (prec)
165  {
166  prec->Mult(r, z); // z = B r
167  nom0 = nom = Dot(z, r);
168  }
169  else
170  {
171  nom0 = nom = Dot(r, r);
172  }
173 
174  if (print_level == 1)
175  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
176  << nom << '\n';
177 
178  r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
179  if (nom <= r0)
180  {
181  converged = 1;
182  final_iter = 0;
183  final_norm = sqrt(nom);
184  return;
185  }
186 
187  // start iteration
188  converged = 0;
190  for (i = 1; true; )
191  {
192  if (prec) // x = x + B (b - A x)
193  {
194  add(x, 1.0, z, x);
195  }
196  else
197  {
198  add(x, 1.0, r, x);
199  }
200 
201  oper->Mult(x, r);
202  subtract(b, r, r); // r = b - A x
203 
204  if (prec)
205  {
206  prec->Mult(r, z); // z = B r
207  nom = Dot(z, r);
208  }
209  else
210  {
211  nom = Dot(r, r);
212  }
213 
214  cf = sqrt(nom/nomold);
215  if (print_level == 1)
216  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
217  << nom << "\tConv. rate: " << cf << '\n';
218  nomold = nom;
219 
220  if (nom < r0)
221  {
222  if (print_level == 2)
223  mfem::out << "Number of SLI iterations: " << i << '\n'
224  << "Conv. rate: " << cf << '\n';
225  else if (print_level == 3)
226  mfem::out << "(B r_0, r_0) = " << nom0 << '\n'
227  << "(B r_N, r_N) = " << nom << '\n'
228  << "Number of SLI iterations: " << i << '\n';
229  converged = 1;
230  final_iter = i;
231  break;
232  }
233 
234  if (++i > max_iter)
235  {
236  break;
237  }
238  }
239 
240  if (print_level >= 0 && !converged)
241  {
242  mfem::err << "SLI: No convergence!" << '\n';
243  mfem::out << "(B r_0, r_0) = " << nom0 << '\n'
244  << "(B r_N, r_N) = " << nom << '\n'
245  << "Number of SLI iterations: " << final_iter << '\n';
246  }
247  if (print_level >= 1 || (print_level >= 0 && !converged))
248  {
249  mfem::out << "Average reduction factor = "
250  << pow (nom/nom0, 0.5/final_iter) << '\n';
251  }
252  final_norm = sqrt(nom);
253 }
254 
255 void SLI(const Operator &A, const Vector &b, Vector &x,
256  int print_iter, int max_num_iter,
257  double RTOLERANCE, double ATOLERANCE)
258 {
259  SLISolver sli;
260  sli.SetPrintLevel(print_iter);
261  sli.SetMaxIter(max_num_iter);
262  sli.SetRelTol(sqrt(RTOLERANCE));
263  sli.SetAbsTol(sqrt(ATOLERANCE));
264  sli.SetOperator(A);
265  sli.Mult(b, x);
266 }
267 
268 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
269  int print_iter, int max_num_iter,
270  double RTOLERANCE, double ATOLERANCE)
271 {
272  SLISolver sli;
273  sli.SetPrintLevel(print_iter);
274  sli.SetMaxIter(max_num_iter);
275  sli.SetRelTol(sqrt(RTOLERANCE));
276  sli.SetAbsTol(sqrt(ATOLERANCE));
277  sli.SetOperator(A);
278  sli.SetPreconditioner(B);
279  sli.Mult(b, x);
280 }
281 
282 
284 {
285  r.SetSize(width);
286  d.SetSize(width);
287  z.SetSize(width);
288 }
289 
290 void CGSolver::Mult(const Vector &b, Vector &x) const
291 {
292  int i;
293  double r0, den, nom, nom0, betanom, alpha, beta;
294 
295  if (iterative_mode)
296  {
297  oper->Mult(x, r);
298  subtract(b, r, r); // r = b - A x
299  }
300  else
301  {
302  r = b;
303  x = 0.0;
304  }
305 
306  if (prec)
307  {
308  prec->Mult(r, z); // z = B r
309  d = z;
310  }
311  else
312  {
313  d = r;
314  }
315  nom0 = nom = Dot(d, r);
316  MFEM_ASSERT(IsFinite(nom), "nom = " << nom);
317 
318  if (print_level == 1 || print_level == 3)
319  {
320  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
321  << nom << (print_level == 3 ? " ...\n" : "\n");
322  }
323 
324  r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
325  if (nom <= r0)
326  {
327  converged = 1;
328  final_iter = 0;
329  final_norm = sqrt(nom);
330  return;
331  }
332 
333  oper->Mult(d, z); // z = A d
334  den = Dot(z, d);
335  MFEM_ASSERT(IsFinite(den), "den = " << den);
336  if (den <= 0.0)
337  {
338  if (Dot(d, d) > 0.0 && print_level >= 0)
339  {
340  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
341  << den << '\n';
342  }
343  if (den == 0.0)
344  {
345  converged = 0;
346  final_iter = 0;
347  final_norm = sqrt(nom);
348  return;
349  }
350  }
351 
352  // start iteration
353  converged = 0;
355  for (i = 1; true; )
356  {
357  alpha = nom/den;
358  add(x, alpha, d, x); // x = x + alpha d
359  add(r, -alpha, z, r); // r = r - alpha A d
360 
361  if (prec)
362  {
363  prec->Mult(r, z); // z = B r
364  betanom = Dot(r, z);
365  }
366  else
367  {
368  betanom = Dot(r, r);
369  }
370  MFEM_ASSERT(IsFinite(betanom), "betanom = " << betanom);
371 
372  if (print_level == 1)
373  {
374  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
375  << betanom << '\n';
376  }
377 
378  if (betanom < r0)
379  {
380  if (print_level == 2)
381  {
382  mfem::out << "Number of PCG iterations: " << i << '\n';
383  }
384  else if (print_level == 3)
385  {
386  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
387  << betanom << '\n';
388  }
389  converged = 1;
390  final_iter = i;
391  break;
392  }
393 
394  if (++i > max_iter)
395  {
396  break;
397  }
398 
399  beta = betanom/nom;
400  if (prec)
401  {
402  add(z, beta, d, d); // d = z + beta d
403  }
404  else
405  {
406  add(r, beta, d, d);
407  }
408  oper->Mult(d, z); // z = A d
409  den = Dot(d, z);
410  MFEM_ASSERT(IsFinite(den), "den = " << den);
411  if (den <= 0.0)
412  {
413  if (Dot(d, d) > 0.0 && print_level >= 0)
414  {
415  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
416  << den << '\n';
417  }
418  if (den == 0.0)
419  {
420  final_iter = i;
421  break;
422  }
423  }
424  nom = betanom;
425  }
426  if (print_level >= 0 && !converged)
427  {
428  if (print_level != 1)
429  {
430  if (print_level != 3)
431  {
432  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
433  << nom0 << " ...\n";
434  }
435  mfem::out << " Iteration : " << setw(3) << final_iter << " (B r, r) = "
436  << betanom << '\n';
437  }
438  mfem::out << "PCG: No convergence!" << '\n';
439  }
440  if (print_level >= 1 || (print_level >= 0 && !converged))
441  {
442  mfem::out << "Average reduction factor = "
443  << pow (betanom/nom0, 0.5/final_iter) << '\n';
444  }
445  final_norm = sqrt(betanom);
446 }
447 
448 void CG(const Operator &A, const Vector &b, Vector &x,
449  int print_iter, int max_num_iter,
450  double RTOLERANCE, double ATOLERANCE)
451 {
452  CGSolver cg;
453  cg.SetPrintLevel(print_iter);
454  cg.SetMaxIter(max_num_iter);
455  cg.SetRelTol(sqrt(RTOLERANCE));
456  cg.SetAbsTol(sqrt(ATOLERANCE));
457  cg.SetOperator(A);
458  cg.Mult(b, x);
459 }
460 
461 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
462  int print_iter, int max_num_iter,
463  double RTOLERANCE, double ATOLERANCE)
464 {
465  CGSolver pcg;
466  pcg.SetPrintLevel(print_iter);
467  pcg.SetMaxIter(max_num_iter);
468  pcg.SetRelTol(sqrt(RTOLERANCE));
469  pcg.SetAbsTol(sqrt(ATOLERANCE));
470  pcg.SetOperator(A);
471  pcg.SetPreconditioner(B);
472  pcg.Mult(b, x);
473 }
474 
475 
476 inline void GeneratePlaneRotation(double &dx, double &dy,
477  double &cs, double &sn)
478 {
479  if (dy == 0.0)
480  {
481  cs = 1.0;
482  sn = 0.0;
483  }
484  else if (fabs(dy) > fabs(dx))
485  {
486  double temp = dx / dy;
487  sn = 1.0 / sqrt( 1.0 + temp*temp );
488  cs = temp * sn;
489  }
490  else
491  {
492  double temp = dy / dx;
493  cs = 1.0 / sqrt( 1.0 + temp*temp );
494  sn = temp * cs;
495  }
496 }
497 
498 inline void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
499 {
500  double temp = cs * dx + sn * dy;
501  dy = -sn * dx + cs * dy;
502  dx = temp;
503 }
504 
505 inline void Update(Vector &x, int k, DenseMatrix &h, Vector &s,
506  Array<Vector*> &v)
507 {
508  Vector y(s);
509 
510  // Backsolve:
511  for (int i = k; i >= 0; i--)
512  {
513  y(i) /= h(i,i);
514  for (int j = i - 1; j >= 0; j--)
515  {
516  y(j) -= h(j,i) * y(i);
517  }
518  }
519 
520  for (int j = 0; j <= k; j++)
521  {
522  x.Add(y(j), *v[j]);
523  }
524 }
525 
526 void GMRESSolver::Mult(const Vector &b, Vector &x) const
527 {
528  // Generalized Minimum Residual method following the algorithm
529  // on p. 20 of the SIAM Templates book.
530 
531  int n = width;
532 
533  DenseMatrix H(m+1, m);
534  Vector s(m+1), cs(m+1), sn(m+1);
535  Vector r(n), w(n);
536  Array<Vector *> v;
537 
538  double resid;
539  int i, j, k;
540 
541  if (iterative_mode)
542  {
543  oper->Mult(x, r);
544  }
545  else
546  {
547  x = 0.0;
548  }
549 
550  if (prec)
551  {
552  if (iterative_mode)
553  {
554  subtract(b, r, w);
555  prec->Mult(w, r); // r = M (b - A x)
556  }
557  else
558  {
559  prec->Mult(b, r);
560  }
561  }
562  else
563  {
564  if (iterative_mode)
565  {
566  subtract(b, r, r);
567  }
568  else
569  {
570  r = b;
571  }
572  }
573  double beta = Norm(r); // beta = ||r||
574  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
575 
576  final_norm = std::max(rel_tol*beta, abs_tol);
577 
578  if (beta <= final_norm)
579  {
580  final_norm = beta;
581  final_iter = 0;
582  converged = 1;
583  goto finish;
584  }
585 
586  if (print_level == 1 || print_level == 3)
587  {
588  mfem::out << " Pass : " << setw(2) << 1
589  << " Iteration : " << setw(3) << 0
590  << " ||B r|| = " << beta << (print_level == 3 ? " ...\n" : "\n");
591  }
592 
593  v.SetSize(m+1, NULL);
594 
595  for (j = 1; j <= max_iter; )
596  {
597  if (v[0] == NULL) { v[0] = new Vector(n); }
598  v[0]->Set(1.0/beta, r);
599  s = 0.0; s(0) = beta;
600 
601  for (i = 0; i < m && j <= max_iter; i++, j++)
602  {
603  if (prec)
604  {
605  oper->Mult(*v[i], r);
606  prec->Mult(r, w); // w = M A v[i]
607  }
608  else
609  {
610  oper->Mult(*v[i], w);
611  }
612 
613  for (k = 0; k <= i; k++)
614  {
615  H(k,i) = Dot(w, *v[k]); // H(k,i) = w * v[k]
616  w.Add(-H(k,i), *v[k]); // w -= H(k,i) * v[k]
617  }
618 
619  H(i+1,i) = Norm(w); // H(i+1,i) = ||w||
620  MFEM_ASSERT(IsFinite(H(i+1,i)), "Norm(w) = " << H(i+1,i));
621  if (v[i+1] == NULL) { v[i+1] = new Vector(n); }
622  v[i+1]->Set(1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
623 
624  for (k = 0; k < i; k++)
625  {
626  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
627  }
628 
629  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
630  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
631  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
632 
633  resid = fabs(s(i+1));
634  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
635 
636  if (resid <= final_norm)
637  {
638  Update(x, i, H, s, v);
639  final_norm = resid;
640  final_iter = j;
641  converged = 1;
642  goto finish;
643  }
644 
645  if (print_level == 1)
646  {
647  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
648  << " Iteration : " << setw(3) << j
649  << " ||B r|| = " << resid << '\n';
650  }
651  }
652 
653  if (print_level == 1 && j <= max_iter)
654  {
655  mfem::out << "Restarting..." << '\n';
656  }
657 
658  Update(x, i-1, H, s, v);
659 
660  oper->Mult(x, r);
661  if (prec)
662  {
663  subtract(b, r, w);
664  prec->Mult(w, r); // r = M (b - A x)
665  }
666  else
667  {
668  subtract(b, r, r);
669  }
670  beta = Norm(r); // beta = ||r||
671  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
672  if (beta <= final_norm)
673  {
674  final_norm = beta;
675  final_iter = j;
676  converged = 1;
677  goto finish;
678  }
679  }
680 
681  final_norm = beta;
683  converged = 0;
684 
685 finish:
686  if (print_level == 1 || print_level == 3)
687  {
688  mfem::out << " Pass : " << setw(2) << (final_iter-1)/m+1
689  << " Iteration : " << setw(3) << final_iter
690  << " ||B r|| = " << final_norm << '\n';
691  }
692  else if (print_level == 2)
693  {
694  mfem::out << "GMRES: Number of iterations: " << final_iter << '\n';
695  }
696  if (print_level >= 0 && !converged)
697  {
698  mfem::out << "GMRES: No convergence!\n";
699  }
700  for (i = 0; i < v.Size(); i++)
701  {
702  delete v[i];
703  }
704 }
705 
706 void FGMRESSolver::Mult(const Vector &b, Vector &x) const
707 {
708  DenseMatrix H(m+1,m);
709  Vector s(m+1), cs(m+1), sn(m+1);
710  Vector r(b.Size());
711 
712  int i, j, k;
713 
714 
715  if (iterative_mode)
716  {
717  oper->Mult(x, r);
718  subtract(b,r,r);
719  }
720  else
721  {
722  x = 0.;
723  r = b;
724  }
725  double beta = Norm(r); // beta = ||r||
726  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
727 
728  final_norm = std::max(rel_tol*beta, abs_tol);
729 
730  if (beta <= final_norm)
731  {
732  final_norm = beta;
733  final_iter = 0;
734  converged = 1;
735  return;
736  }
737 
738  if (print_level>=0)
739  mfem::out << " Pass : " << setw(2) << 1
740  << " Iteration : " << setw(3) << 0
741  << " || r || = " << beta << endl;
742 
743  Array<Vector*> v(m+1);
744  Array<Vector*> z(m+1);
745  for (i= 0; i<=m; i++)
746  {
747  v[i] = NULL;
748  z[i] = NULL;
749  }
750 
751  j = 1;
752  while (j <= max_iter)
753  {
754  if (v[0] == NULL) { v[0] = new Vector(b.Size()); }
755  (*v[0]) = 0.0;
756  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
757  s = 0.0; s(0) = beta;
758 
759  for (i = 0; i < m && j <= max_iter; i++, j++)
760  {
761 
762  if (z[i] == NULL) { z[i] = new Vector(b.Size()); }
763  (*z[i]) = 0.0;
764 
765  if (prec)
766  {
767  prec->Mult(*v[i], *z[i]);
768  }
769  else
770  {
771  (*z[i]) = (*v[i]);
772  }
773  oper->Mult(*z[i], r);
774 
775  for (k = 0; k <= i; k++)
776  {
777  H(k,i) = Dot( r, *v[k]); // H(k,i) = r * v[k]
778  r.Add(-H(k,i), (*v[k])); // r -= H(k,i) * v[k]
779  }
780 
781  H(i+1,i) = Norm(r); // H(i+1,i) = ||r||
782  if (v[i+1] == NULL) { v[i+1] = new Vector(b.Size()); }
783  (*v[i+1]) = 0.0;
784  v[i+1] -> Add (1.0/H(i+1,i), r); // v[i+1] = r / H(i+1,i)
785 
786  for (k = 0; k < i; k++)
787  {
788  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
789  }
790 
791  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
792  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
793  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
794 
795  double resid = fabs(s(i+1));
796  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
797  if (print_level >= 0)
798  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
799  << " Iteration : " << setw(3) << j
800  << " || r || = " << resid << endl;
801 
802  if ( resid <= final_norm)
803  {
804  Update(x, i, H, s, z);
805  final_norm = resid;
806  final_iter = j;
807  converged = 1;
808  for (i= 0; i<=m; i++)
809  {
810  if (v[i]) { delete v[i]; }
811  if (z[i]) { delete z[i]; }
812  }
813  return;
814  }
815  }
816 
817  if (print_level>=0)
818  {
819  mfem::out << "Restarting..." << endl;
820  }
821 
822  Update(x, i-1, H, s, z);
823 
824  oper->Mult(x, r);
825  subtract(b,r,r);
826  beta = Norm(r);
827  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
828  if ( beta <= final_norm)
829  {
830  final_norm = beta;
831  final_iter = j;
832  converged = 1;
833  for (i= 0; i<=m; i++)
834  {
835  if (v[i]) { delete v[i]; }
836  if (z[i]) { delete z[i]; }
837  }
838  return;
839  }
840  }
841 
842  for (i = 0; i <= m; i++)
843  {
844  if (v[i]) { delete v[i]; }
845  if (z[i]) { delete z[i]; }
846  }
847  converged = 0;
848  return;
849 
850 }
851 
852 
853 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
854  int &max_iter, int m, double &tol, double atol, int printit)
855 {
856  GMRESSolver gmres;
857  gmres.SetPrintLevel(printit);
858  gmres.SetMaxIter(max_iter);
859  gmres.SetKDim(m);
860  gmres.SetRelTol(sqrt(tol));
861  gmres.SetAbsTol(sqrt(atol));
862  gmres.SetOperator(A);
863  gmres.SetPreconditioner(M);
864  gmres.Mult(b, x);
865  max_iter = gmres.GetNumIterations();
866  tol = gmres.GetFinalNorm()*gmres.GetFinalNorm();
867  return gmres.GetConverged();
868 }
869 
870 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
871  int print_iter, int max_num_iter, int m, double rtol, double atol)
872 {
873  GMRES(A, x, b, B, max_num_iter, m, rtol, atol, print_iter);
874 }
875 
876 
878 {
879  p.SetSize(width);
880  phat.SetSize(width);
881  s.SetSize(width);
882  shat.SetSize(width);
883  t.SetSize(width);
884  v.SetSize(width);
885  r.SetSize(width);
887 }
888 
889 void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
890 {
891  // BiConjugate Gradient Stabilized method following the algorithm
892  // on p. 27 of the SIAM Templates book.
893 
894  int i;
895  double resid, tol_goal;
896  double rho_1, rho_2=1.0, alpha=1.0, beta, omega=1.0;
897 
898  if (iterative_mode)
899  {
900  oper->Mult(x, r);
901  subtract(b, r, r); // r = b - A x
902  }
903  else
904  {
905  x = 0.0;
906  r = b;
907  }
908  rtilde = r;
909 
910  resid = Norm(r);
911  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
912  if (print_level >= 0)
913  mfem::out << " Iteration : " << setw(3) << 0
914  << " ||r|| = " << resid << '\n';
915 
916  tol_goal = std::max(resid*rel_tol, abs_tol);
917 
918  if (resid <= tol_goal)
919  {
920  final_norm = resid;
921  final_iter = 0;
922  converged = 1;
923  return;
924  }
925 
926  for (i = 1; i <= max_iter; i++)
927  {
928  rho_1 = Dot(rtilde, r);
929  if (rho_1 == 0)
930  {
931  if (print_level >= 0)
932  mfem::out << " Iteration : " << setw(3) << i
933  << " ||r|| = " << resid << '\n';
934  final_norm = resid;
935  final_iter = i;
936  converged = 0;
937  return;
938  }
939  if (i == 1)
940  {
941  p = r;
942  }
943  else
944  {
945  beta = (rho_1/rho_2) * (alpha/omega);
946  add(p, -omega, v, p); // p = p - omega * v
947  add(r, beta, p, p); // p = r + beta * p
948  }
949  if (prec)
950  {
951  prec->Mult(p, phat); // phat = M^{-1} * p
952  }
953  else
954  {
955  phat = p;
956  }
957  oper->Mult(phat, v); // v = A * phat
958  alpha = rho_1 / Dot(rtilde, v);
959  add(r, -alpha, v, s); // s = r - alpha * v
960  resid = Norm(s);
961  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
962  if (resid < tol_goal)
963  {
964  x.Add(alpha, phat); // x = x + alpha * phat
965  if (print_level >= 0)
966  mfem::out << " Iteration : " << setw(3) << i
967  << " ||s|| = " << resid << '\n';
968  final_norm = resid;
969  final_iter = i;
970  converged = 1;
971  return;
972  }
973  if (print_level >= 0)
974  mfem::out << " Iteration : " << setw(3) << i
975  << " ||s|| = " << resid;
976  if (prec)
977  {
978  prec->Mult(s, shat); // shat = M^{-1} * s
979  }
980  else
981  {
982  shat = s;
983  }
984  oper->Mult(shat, t); // t = A * shat
985  omega = Dot(t, s) / Dot(t, t);
986  x.Add(alpha, phat); // x += alpha * phat
987  x.Add(omega, shat); // x += omega * shat
988  add(s, -omega, t, r); // r = s - omega * t
989 
990  rho_2 = rho_1;
991  resid = Norm(r);
992  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
993  if (print_level >= 0)
994  {
995  mfem::out << " ||r|| = " << resid << '\n';
996  }
997  if (resid < tol_goal)
998  {
999  final_norm = resid;
1000  final_iter = i;
1001  converged = 1;
1002  return;
1003  }
1004  if (omega == 0)
1005  {
1006  final_norm = resid;
1007  final_iter = i;
1008  converged = 0;
1009  return;
1010  }
1011  }
1012 
1013  final_norm = resid;
1014  final_iter = max_iter;
1015  converged = 0;
1016 }
1017 
1018 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
1019  int &max_iter, double &tol, double atol, int printit)
1020 {
1021  BiCGSTABSolver bicgstab;
1022  bicgstab.SetPrintLevel(printit);
1023  bicgstab.SetMaxIter(max_iter);
1024  bicgstab.SetRelTol(sqrt(tol));
1025  bicgstab.SetAbsTol(sqrt(atol));
1026  bicgstab.SetOperator(A);
1027  bicgstab.SetPreconditioner(M);
1028  bicgstab.Mult(b, x);
1029  max_iter = bicgstab.GetNumIterations();
1030  tol = bicgstab.GetFinalNorm()*bicgstab.GetFinalNorm();
1031  return bicgstab.GetConverged();
1032 }
1033 
1034 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
1035  int print_iter, int max_num_iter, double rtol, double atol)
1036 {
1037  BiCGSTAB(A, x, b, B, max_num_iter, rtol, atol, print_iter);
1038 }
1039 
1040 
1042 {
1044  v0.SetSize(width);
1045  v1.SetSize(width);
1046  w0.SetSize(width);
1047  w1.SetSize(width);
1048  q.SetSize(width);
1049  if (prec)
1050  {
1051  u1.SetSize(width);
1052  }
1053 }
1054 
1055 void MINRESSolver::Mult(const Vector &b, Vector &x) const
1056 {
1057  // Based on the MINRES algorithm on p. 86, Fig. 6.9 in
1058  // "Iterative Krylov Methods for Large Linear Systems",
1059  // by Henk A. van der Vorst, 2003.
1060  // Extended to support an SPD preconditioner.
1061 
1062  int it;
1063  double beta, eta, gamma0, gamma1, sigma0, sigma1;
1064  double alpha, delta, rho1, rho2, rho3, norm_goal;
1065  Vector *z = (prec) ? &u1 : &v1;
1066 
1067  converged = 1;
1068 
1069  if (!iterative_mode)
1070  {
1071  v1 = b;
1072  x = 0.;
1073  }
1074  else
1075  {
1076  oper->Mult(x, v1);
1077  subtract(b, v1, v1);
1078  }
1079 
1080  if (prec)
1081  {
1082  prec->Mult(v1, u1);
1083  }
1084  eta = beta = sqrt(Dot(*z, v1));
1085  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1086  gamma0 = gamma1 = 1.;
1087  sigma0 = sigma1 = 0.;
1088 
1089  norm_goal = std::max(rel_tol*eta, abs_tol);
1090 
1091  if (eta <= norm_goal)
1092  {
1093  it = 0;
1094  goto loop_end;
1095  }
1096 
1097  if (print_level == 1 || print_level == 3)
1098  {
1099  mfem::out << "MINRES: iteration " << setw(3) << 0 << ": ||r||_B = "
1100  << eta << (print_level == 3 ? " ...\n" : "\n");
1101  }
1102 
1103  for (it = 1; it <= max_iter; it++)
1104  {
1105  v1 /= beta;
1106  if (prec)
1107  {
1108  u1 /= beta;
1109  }
1110  oper->Mult(*z, q);
1111  alpha = Dot(*z, q);
1112  MFEM_ASSERT(IsFinite(alpha), "alpha = " << alpha);
1113  if (it > 1) // (v0 == 0) for (it == 1)
1114  {
1115  q.Add(-beta, v0);
1116  }
1117  add(q, -alpha, v1, v0);
1118 
1119  delta = gamma1*alpha - gamma0*sigma1*beta;
1120  rho3 = sigma0*beta;
1121  rho2 = sigma1*alpha + gamma0*gamma1*beta;
1122  if (!prec)
1123  {
1124  beta = Norm(v0);
1125  }
1126  else
1127  {
1128  prec->Mult(v0, q);
1129  beta = sqrt(Dot(v0, q));
1130  }
1131  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1132  rho1 = hypot(delta, beta);
1133 
1134  if (it == 1)
1135  {
1136  w0.Set(1./rho1, *z); // (w0 == 0) and (w1 == 0)
1137  }
1138  else if (it == 2)
1139  {
1140  add(1./rho1, *z, -rho2/rho1, w1, w0); // (w0 == 0)
1141  }
1142  else
1143  {
1144  add(-rho3/rho1, w0, -rho2/rho1, w1, w0);
1145  w0.Add(1./rho1, *z);
1146  }
1147 
1148  gamma0 = gamma1;
1149  gamma1 = delta/rho1;
1150 
1151  x.Add(gamma1*eta, w0);
1152 
1153  sigma0 = sigma1;
1154  sigma1 = beta/rho1;
1155 
1156  eta = -sigma1*eta;
1157  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1158 
1159  if (fabs(eta) <= norm_goal)
1160  {
1161  goto loop_end;
1162  }
1163 
1164  if (print_level == 1)
1165  {
1166  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1167  << fabs(eta) << '\n';
1168  }
1169 
1170  if (prec)
1171  {
1172  Swap(u1, q);
1173  }
1174  Swap(v0, v1);
1175  Swap(w0, w1);
1176  }
1177  converged = 0;
1178  it--;
1179 
1180 loop_end:
1181  final_iter = it;
1182  final_norm = fabs(eta);
1183 
1184  if (print_level == 1 || print_level == 3)
1185  {
1186  mfem::out << "MINRES: iteration " << setw(3) << final_iter << ": ||r||_B = "
1187  << final_norm << '\n';
1188  }
1189  else if (print_level == 2)
1190  {
1191  mfem::out << "MINRES: number of iterations: " << final_iter << '\n';
1192  }
1193 #if 0
1194  if (print_level >= 1)
1195  {
1196  oper->Mult(x, v1);
1197  subtract(b, v1, v1);
1198  if (prec)
1199  {
1200  prec->Mult(v1, u1);
1201  }
1202  eta = sqrt(Dot(*z, v1));
1203  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1204  << eta << " (re-computed)" << '\n';
1205  }
1206 #endif
1207  if (!converged && print_level >= 0)
1208  {
1209  mfem::out << "MINRES: No convergence!\n";
1210  }
1211 }
1212 
1213 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it,
1214  int max_it, double rtol, double atol)
1215 {
1216  MINRESSolver minres;
1217  minres.SetPrintLevel(print_it);
1218  minres.SetMaxIter(max_it);
1219  minres.SetRelTol(sqrt(rtol));
1220  minres.SetAbsTol(sqrt(atol));
1221  minres.SetOperator(A);
1222  minres.Mult(b, x);
1223 }
1224 
1225 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
1226  int print_it, int max_it, double rtol, double atol)
1227 {
1228  MINRESSolver minres;
1229  minres.SetPrintLevel(print_it);
1230  minres.SetMaxIter(max_it);
1231  minres.SetRelTol(sqrt(rtol));
1232  minres.SetAbsTol(sqrt(atol));
1233  minres.SetOperator(A);
1234  minres.SetPreconditioner(B);
1235  minres.Mult(b, x);
1236 }
1237 
1238 
1240 {
1241  oper = &op;
1242  height = op.Height();
1243  width = op.Width();
1244  MFEM_ASSERT(height == width, "square Operator is required.");
1245 
1246  r.SetSize(width);
1247  c.SetSize(width);
1248 }
1249 
1250 void NewtonSolver::Mult(const Vector &b, Vector &x) const
1251 {
1252  MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
1253  MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver).");
1254 
1255  int it;
1256  double norm0, norm, norm_goal;
1257  const bool have_b = (b.Size() == Height());
1258 
1259  if (!iterative_mode)
1260  {
1261  x = 0.0;
1262  }
1263 
1264  oper->Mult(x, r);
1265  if (have_b)
1266  {
1267  r -= b;
1268  }
1269 
1270  norm0 = norm = Norm(r);
1271  norm_goal = std::max(rel_tol*norm, abs_tol);
1272 
1273  prec->iterative_mode = false;
1274 
1275  // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b]
1276  for (it = 0; true; it++)
1277  {
1278  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
1279  if (print_level >= 0)
1280  {
1281  mfem::out << "Newton iteration " << setw(2) << it
1282  << " : ||r|| = " << norm;
1283  if (it > 0)
1284  {
1285  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
1286  }
1287  mfem::out << '\n';
1288  }
1289 
1290  if (norm <= norm_goal)
1291  {
1292  converged = 1;
1293  break;
1294  }
1295 
1296  if (it >= max_iter)
1297  {
1298  converged = 0;
1299  break;
1300  }
1301 
1303 
1304  prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b]
1305 
1306  const double c_scale = ComputeScalingFactor(x, b);
1307  if (c_scale == 0.0)
1308  {
1309  converged = 0;
1310  break;
1311  }
1312  add(x, -c_scale, c, x);
1313 
1314  oper->Mult(x, r);
1315  if (have_b)
1316  {
1317  r -= b;
1318  }
1319  norm = Norm(r);
1320  }
1321 
1322  final_iter = it;
1323  final_norm = norm;
1324 }
1325 
1326 
1327 int aGMRES(const Operator &A, Vector &x, const Vector &b,
1328  const Operator &M, int &max_iter,
1329  int m_max, int m_min, int m_step, double cf,
1330  double &tol, double &atol, int printit)
1331 {
1332  int n = A.Width();
1333 
1334  int m = m_max;
1335 
1336  DenseMatrix H(m+1,m);
1337  Vector s(m+1), cs(m+1), sn(m+1);
1338  Vector w(n), av(n);
1339 
1340  double r1, resid;
1341  int i, j, k;
1342 
1343  M.Mult(b,w);
1344  double normb = w.Norml2(); // normb = ||M b||
1345  if (normb == 0.0)
1346  {
1347  normb = 1;
1348  }
1349 
1350  Vector r(n);
1351  A.Mult(x, r);
1352  subtract(b,r,w);
1353  M.Mult(w, r); // r = M (b - A x)
1354  double beta = r.Norml2(); // beta = ||r||
1355 
1356  resid = beta / normb;
1357 
1358  if (resid * resid <= tol)
1359  {
1360  tol = resid * resid;
1361  max_iter = 0;
1362  return 0;
1363  }
1364 
1365  if (printit)
1366  mfem::out << " Pass : " << setw(2) << 1
1367  << " Iteration : " << setw(3) << 0
1368  << " (r, r) = " << beta*beta << '\n';
1369 
1370  tol *= (normb*normb);
1371  tol = (atol > tol) ? atol : tol;
1372 
1373  m = m_max;
1374  Array<Vector *> v(m+1);
1375  for (i= 0; i<=m; i++)
1376  {
1377  v[i] = new Vector(n);
1378  (*v[i]) = 0.0;
1379  }
1380 
1381  j = 1;
1382  while (j <= max_iter)
1383  {
1384  (*v[0]) = 0.0;
1385  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
1386  s = 0.0; s(0) = beta;
1387 
1388  r1 = beta;
1389 
1390  for (i = 0; i < m && j <= max_iter; i++)
1391  {
1392  A.Mult((*v[i]),av);
1393  M.Mult(av,w); // w = M A v[i]
1394 
1395  for (k = 0; k <= i; k++)
1396  {
1397  H(k,i) = w * (*v[k]); // H(k,i) = w * v[k]
1398  w.Add(-H(k,i), (*v[k])); // w -= H(k,i) * v[k]
1399  }
1400 
1401  H(i+1,i) = w.Norml2(); // H(i+1,i) = ||w||
1402  (*v[i+1]) = 0.0;
1403  v[i+1] -> Add (1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
1404 
1405  for (k = 0; k < i; k++)
1406  {
1407  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
1408  }
1409 
1410  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1411  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1412  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
1413 
1414  resid = fabs(s(i+1));
1415  if (printit)
1416  mfem::out << " Pass : " << setw(2) << j
1417  << " Iteration : " << setw(3) << i+1
1418  << " (r, r) = " << resid*resid << '\n';
1419 
1420  if ( resid*resid < tol)
1421  {
1422  Update(x, i, H, s, v);
1423  tol = resid * resid;
1424  max_iter = j;
1425  for (i= 0; i<=m; i++)
1426  {
1427  delete v[i];
1428  }
1429  return 0;
1430  }
1431  }
1432 
1433  if (printit)
1434  {
1435  mfem::out << "Restarting..." << '\n';
1436  }
1437 
1438  Update(x, i-1, H, s, v);
1439 
1440  A.Mult(x, r);
1441  subtract(b,r,w);
1442  M.Mult(w, r); // r = M (b - A x)
1443  beta = r.Norml2(); // beta = ||r||
1444  if ( resid*resid < tol)
1445  {
1446  tol = resid * resid;
1447  max_iter = j;
1448  for (i= 0; i<=m; i++)
1449  {
1450  delete v[i];
1451  }
1452  return 0;
1453  }
1454 
1455  if (beta/r1 > cf)
1456  {
1457  if (m - m_step >= m_min)
1458  {
1459  m -= m_step;
1460  }
1461  else
1462  {
1463  m = m_max;
1464  }
1465  }
1466 
1467  j++;
1468  }
1469 
1470  tol = resid * resid;
1471  for (i= 0; i<=m; i++)
1472  {
1473  delete v[i];
1474  }
1475  return 1;
1476 }
1477 
1478 
1479 void SLBQPOptimizer::SetBounds(const Vector &_lo, const Vector &_hi)
1480 {
1481  lo.SetDataAndSize(_lo.GetData(), _lo.Size());
1482  hi.SetDataAndSize(_hi.GetData(), _hi.Size());
1483 }
1484 
1486 {
1487  w.SetDataAndSize(_w.GetData(), _w.Size());
1488  a = _a;
1489 }
1490 
1492 {
1493  mfem_error("SLBQPOptimizer::SetPreconditioner() : "
1494  "not meaningful for this solver");
1495 }
1496 
1498 {
1499  mfem_error("SLBQPOptimizer::SetOperator() : "
1500  "not meaningful for this solver");
1501 }
1502 
1503 inline void SLBQPOptimizer::print_iteration(int it, double r, double l) const
1504 {
1505  if (print_level > 1)
1506  mfem::out << "SLBQP iteration " << it << ": residual = " << r
1507  << ", lambda = " << l << '\n';
1508 }
1509 
1510 void SLBQPOptimizer::Mult(const Vector& xt, Vector& x) const
1511 {
1512  // Based on code provided by Denis Ridzal, dridzal@sandia.gov.
1513  // Algorithm adapted from Dai and Fletcher, "New Algorithms for
1514  // Singly Linearly Constrained Quadratic Programs Subject to Lower
1515  // and Upper Bounds", Numerical Analysis Report NA/216, 2003.
1516 
1517  // Set some algorithm-specific constants and temporaries.
1518  int nclip = 0;
1519  double l = 0;
1520  double llow = 0;
1521  double lupp = 0;
1522  double lnew = 0;
1523  double dl = 2;
1524  double r = 0;
1525  double rlow = 0;
1526  double rupp = 0;
1527  double s = 0;
1528 
1529  const double smin = 0.1;
1530 
1531  const double tol = max(abs_tol, rel_tol*a);
1532 
1533  // *** Start bracketing phase of SLBQP ***
1534  if (print_level > 1)
1535  {
1536  mfem::out << "SLBQP bracketing phase" << '\n';
1537  }
1538 
1539  // Solve QP with fixed Lagrange multiplier
1540  r = solve(l,xt,x,nclip);
1541  print_iteration(nclip, r, l);
1542 
1543 
1544  // If x=xt was already within bounds and satisfies the linear
1545  // constraint, then we already have the solution.
1546  if (fabs(r) <= tol)
1547  {
1548  converged = true;
1549  goto slbqp_done;
1550  }
1551 
1552  if (r < 0)
1553  {
1554  llow = l; rlow = r; l = l + dl;
1555 
1556  // Solve QP with fixed Lagrange multiplier
1557  r = solve(l,xt,x,nclip);
1558  print_iteration(nclip, r, l);
1559 
1560  while ((r < 0) && (nclip < max_iter))
1561  {
1562  llow = l;
1563  s = rlow/r - 1.0;
1564  if (s < smin) { s = smin; }
1565  dl = dl + dl/s;
1566  l = l + dl;
1567 
1568  // Solve QP with fixed Lagrange multiplier
1569  r = solve(l,xt,x,nclip);
1570  print_iteration(nclip, r, l);
1571  }
1572 
1573  lupp = l; rupp = r;
1574  }
1575  else
1576  {
1577  lupp = l; rupp = r; l = l - dl;
1578 
1579  // Solve QP with fixed Lagrange multiplier
1580  r = solve(l,xt,x,nclip);
1581  print_iteration(nclip, r, l);
1582 
1583  while ((r > 0) && (nclip < max_iter))
1584  {
1585  lupp = l;
1586  s = rupp/r - 1.0;
1587  if (s < smin) { s = smin; }
1588  dl = dl + dl/s;
1589  l = l - dl;
1590 
1591  // Solve QP with fixed Lagrange multiplier
1592  r = solve(l,xt,x,nclip);
1593  print_iteration(nclip, r, l);
1594  }
1595 
1596  llow = l; rlow = r;
1597  }
1598 
1599  // *** Stop bracketing phase of SLBQP ***
1600 
1601 
1602  // *** Start secant phase of SLBQP ***
1603  if (print_level > 1)
1604  {
1605  mfem::out << "SLBQP secant phase" << '\n';
1606  }
1607 
1608  s = 1.0 - rlow/rupp; dl = dl/s; l = lupp - dl;
1609 
1610  // Solve QP with fixed Lagrange multiplier
1611  r = solve(l,xt,x,nclip);
1612  print_iteration(nclip, r, l);
1613 
1614  while ( (fabs(r) > tol) && (nclip < max_iter) )
1615  {
1616  if (r > 0)
1617  {
1618  if (s <= 2.0)
1619  {
1620  lupp = l; rupp = r; s = 1.0 - rlow/rupp;
1621  dl = (lupp - llow)/s; l = lupp - dl;
1622  }
1623  else
1624  {
1625  s = rupp/r - 1.0;
1626  if (s < smin) { s = smin; }
1627  dl = (lupp - l)/s;
1628  lnew = 0.75*llow + 0.25*l;
1629  if (lnew < l-dl) { lnew = l-dl; }
1630  lupp = l; rupp = r; l = lnew;
1631  s = (lupp - llow)/(lupp - l);
1632  }
1633 
1634  }
1635  else
1636  {
1637  if (s >= 2.0)
1638  {
1639  llow = l; rlow = r; s = 1.0 - rlow/rupp;
1640  dl = (lupp - llow)/s; l = lupp - dl;
1641  }
1642  else
1643  {
1644  s = rlow/r - 1.0;
1645  if (s < smin) { s = smin; }
1646  dl = (l - llow)/s;
1647  lnew = 0.75*lupp + 0.25*l;
1648  if (lnew < l+dl) { lnew = l+dl; }
1649  llow = l; rlow = r; l = lnew;
1650  s = (lupp - llow)/(lupp - l);
1651  }
1652  }
1653 
1654  // Solve QP with fixed Lagrange multiplier
1655  r = solve(l,xt,x,nclip);
1656  print_iteration(nclip, r, l);
1657  }
1658 
1659  // *** Stop secant phase of SLBQP ***
1660 
1661  converged = (fabs(r) <= tol);
1662  if (!converged && print_level >= 0)
1663  {
1664  mfem::err << "SLBQP not converged!" << '\n';
1665  }
1666 
1667 slbqp_done:
1668 
1669  final_iter = nclip;
1670  final_norm = r;
1671 
1672  if (print_level == 1 || (!converged && print_level >= 0))
1673  {
1674  mfem::out << "SLBQP iterations = " << nclip << '\n';
1675  mfem::out << "SLBQP lambda = " << l << '\n';
1676  mfem::out << "SLBQP residual = " << r << '\n';
1677  }
1678 }
1679 
1680 #ifdef MFEM_USE_SUITESPARSE
1681 
1683 {
1684  mat = NULL;
1685  Numeric = NULL;
1686  AI = AJ = NULL;
1687  if (!use_long_ints)
1688  {
1689  umfpack_di_defaults(Control);
1690  }
1691  else
1692  {
1693  umfpack_dl_defaults(Control);
1694  }
1695 }
1696 
1698 {
1699  int *Ap, *Ai;
1700  void *Symbolic;
1701  double *Ax;
1702 
1703  if (Numeric)
1704  {
1705  if (!use_long_ints)
1706  {
1707  umfpack_di_free_numeric(&Numeric);
1708  }
1709  else
1710  {
1711  umfpack_dl_free_numeric(&Numeric);
1712  }
1713  }
1714 
1715  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
1716  MFEM_VERIFY(mat, "not a SparseMatrix");
1717 
1718  // UMFPack requires that the column-indices in mat corresponding to each
1719  // row be sorted.
1720  // Generally, this will modify the ordering of the entries of mat.
1722 
1723  height = mat->Height();
1724  width = mat->Width();
1725  MFEM_VERIFY(width == height, "not a square matrix");
1726 
1727  Ap = mat->GetI();
1728  Ai = mat->GetJ();
1729  Ax = mat->GetData();
1730 
1731  if (!use_long_ints)
1732  {
1733  int status = umfpack_di_symbolic(width, width, Ap, Ai, Ax, &Symbolic,
1734  Control, Info);
1735  if (status < 0)
1736  {
1737  umfpack_di_report_info(Control, Info);
1738  umfpack_di_report_status(Control, status);
1739  mfem_error("UMFPackSolver::SetOperator :"
1740  " umfpack_di_symbolic() failed!");
1741  }
1742 
1743  status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric,
1744  Control, Info);
1745  if (status < 0)
1746  {
1747  umfpack_di_report_info(Control, Info);
1748  umfpack_di_report_status(Control, status);
1749  mfem_error("UMFPackSolver::SetOperator :"
1750  " umfpack_di_numeric() failed!");
1751  }
1752  umfpack_di_free_symbolic(&Symbolic);
1753  }
1754  else
1755  {
1756  SuiteSparse_long status;
1757 
1758  delete [] AJ;
1759  delete [] AI;
1760  AI = new SuiteSparse_long[width + 1];
1761  AJ = new SuiteSparse_long[Ap[width]];
1762  for (int i = 0; i <= width; i++)
1763  {
1764  AI[i] = (SuiteSparse_long)(Ap[i]);
1765  }
1766  for (int i = 0; i < Ap[width]; i++)
1767  {
1768  AJ[i] = (SuiteSparse_long)(Ai[i]);
1769  }
1770 
1771  status = umfpack_dl_symbolic(width, width, AI, AJ, Ax, &Symbolic,
1772  Control, Info);
1773  if (status < 0)
1774  {
1775  umfpack_dl_report_info(Control, Info);
1776  umfpack_dl_report_status(Control, status);
1777  mfem_error("UMFPackSolver::SetOperator :"
1778  " umfpack_dl_symbolic() failed!");
1779  }
1780 
1781  status = umfpack_dl_numeric(AI, AJ, Ax, Symbolic, &Numeric,
1782  Control, Info);
1783  if (status < 0)
1784  {
1785  umfpack_dl_report_info(Control, Info);
1786  umfpack_dl_report_status(Control, status);
1787  mfem_error("UMFPackSolver::SetOperator :"
1788  " umfpack_dl_numeric() failed!");
1789  }
1790  umfpack_dl_free_symbolic(&Symbolic);
1791  }
1792 }
1793 
1794 void UMFPackSolver::Mult(const Vector &b, Vector &x) const
1795 {
1796  if (mat == NULL)
1797  mfem_error("UMFPackSolver::Mult : matrix is not set!"
1798  " Call SetOperator first!");
1799 
1800  if (!use_long_ints)
1801  {
1802  int status =
1803  umfpack_di_solve(UMFPACK_At, mat->GetI(), mat->GetJ(),
1804  mat->GetData(), x, b, Numeric, Control, Info);
1805  umfpack_di_report_info(Control, Info);
1806  if (status < 0)
1807  {
1808  umfpack_di_report_status(Control, status);
1809  mfem_error("UMFPackSolver::Mult : umfpack_di_solve() failed!");
1810  }
1811  }
1812  else
1813  {
1814  SuiteSparse_long status =
1815  umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->GetData(), x, b,
1816  Numeric, Control, Info);
1817  umfpack_dl_report_info(Control, Info);
1818  if (status < 0)
1819  {
1820  umfpack_dl_report_status(Control, status);
1821  mfem_error("UMFPackSolver::Mult : umfpack_dl_solve() failed!");
1822  }
1823  }
1824 }
1825 
1827 {
1828  if (mat == NULL)
1829  mfem_error("UMFPackSolver::MultTranspose : matrix is not set!"
1830  " Call SetOperator first!");
1831 
1832  if (!use_long_ints)
1833  {
1834  int status =
1835  umfpack_di_solve(UMFPACK_A, mat->GetI(), mat->GetJ(),
1836  mat->GetData(), x, b, Numeric, Control, Info);
1837  umfpack_di_report_info(Control, Info);
1838  if (status < 0)
1839  {
1840  umfpack_di_report_status(Control, status);
1841  mfem_error("UMFPackSolver::MultTranspose :"
1842  " umfpack_di_solve() failed!");
1843  }
1844  }
1845  else
1846  {
1847  SuiteSparse_long status =
1848  umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->GetData(), x, b,
1849  Numeric, Control, Info);
1850  umfpack_dl_report_info(Control, Info);
1851  if (status < 0)
1852  {
1853  umfpack_dl_report_status(Control, status);
1854  mfem_error("UMFPackSolver::MultTranspose :"
1855  " umfpack_dl_solve() failed!");
1856  }
1857  }
1858 }
1859 
1861 {
1862  delete [] AJ;
1863  delete [] AI;
1864  if (Numeric)
1865  {
1866  if (!use_long_ints)
1867  {
1868  umfpack_di_free_numeric(&Numeric);
1869  }
1870  else
1871  {
1872  umfpack_dl_free_numeric(&Numeric);
1873  }
1874  }
1875 }
1876 
1878 {
1879  klu_defaults(&Common);
1880 }
1881 
1883 {
1884  if (Numeric)
1885  {
1886  MFEM_ASSERT(Symbolic != 0,
1887  "Had Numeric pointer in KLU, but not Symbolic");
1888  klu_free_symbolic(&Symbolic, &Common);
1889  Symbolic = 0;
1890  klu_free_numeric(&Numeric, &Common);
1891  Numeric = 0;
1892  }
1893 
1894  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
1895  MFEM_VERIFY(mat != NULL, "not a SparseMatrix");
1896 
1897  // KLU requires that the column-indices in mat corresponding to each row be
1898  // sorted. Generally, this will modify the ordering of the entries of mat.
1900 
1901  height = mat->Height();
1902  width = mat->Width();
1903  MFEM_VERIFY(width == height, "not a square matrix");
1904 
1905  int * Ap = mat->GetI();
1906  int * Ai = mat->GetJ();
1907  double * Ax = mat->GetData();
1908 
1909  Symbolic = klu_analyze( height, Ap, Ai, &Common);
1910  Numeric = klu_factor(Ap, Ai, Ax, Symbolic, &Common);
1911 }
1912 
1913 void KLUSolver::Mult(const Vector &b, Vector &x) const
1914 {
1915  MFEM_VERIFY(mat != NULL,
1916  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
1917 
1918  int n = mat->Height();
1919  int numRhs = 1;
1920  // Copy B into X, so we can pass it in and overwrite it.
1921  x = b;
1922  // Solve the transpose, since KLU thinks the matrix is compressed column
1923  // format.
1924  klu_tsolve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
1925 }
1926 
1927 void KLUSolver::MultTranspose(const Vector &b, Vector &x) const
1928 {
1929  MFEM_VERIFY(mat != NULL,
1930  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
1931 
1932  int n = mat->Height();
1933  int numRhs = 1;
1934  // Copy B into X, so we can pass it in and overwrite it.
1935  x = b;
1936  // Solve the regular matrix, not the transpose, since KLU thinks the matrix
1937  // is compressed column format.
1938  klu_solve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
1939 }
1940 
1942 {
1943  klu_free_symbolic (&Symbolic, &Common) ;
1944  klu_free_numeric (&Numeric, &Common) ;
1945  Symbolic = 0;
1946  Numeric = 0;
1947 }
1948 
1949 #endif // MFEM_USE_SUITESPARSE
1950 
1951 }
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:498
int Size() const
Logical size of the array.
Definition: array.hpp:118
Conjugate gradient method.
Definition: solvers.hpp:111
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:1826
double Info[UMFPACK_INFO]
Definition: solvers.hpp:356
SparseMatrix * mat
Definition: solvers.hpp:348
int GetNumIterations() const
Definition: solvers.hpp:66
const Operator * oper
Definition: solvers.hpp:41
SuiteSparse_long * AI
Definition: solvers.hpp:350
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:290
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:400
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:68
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:709
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:146
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:706
int Size() const
Returns the size of the vector.
Definition: vector.hpp:150
int * GetI()
Return the array I.
Definition: sparsemat.hpp:141
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.cpp:1927
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:283
MINRES method.
Definition: solvers.hpp:221
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:159
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:365
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:238
double Norm(const Vector &x) const
Definition: solvers.hpp:52
double GetFinalNorm() const
Definition: solvers.hpp:68
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:151
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:1018
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:112
bool IsFinite(const double &val)
Definition: vector.hpp:365
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:3030
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:526
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1055
void SetLinearConstraint(const Vector &_w, double _a)
Definition: solvers.cpp:1485
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:234
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:67
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:94
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:311
Data type sparse matrix.
Definition: sparsemat.hpp:40
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
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:1213
virtual ~UMFPackSolver()
Definition: solvers.cpp:1860
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:565
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:146
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
SuiteSparse_long * AJ
Definition: solvers.hpp:350
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:156
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:448
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:1503
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:461
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1882
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1041
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:51
int GetConverged() const
Definition: solvers.hpp:67
klu_symbolic * Symbolic
Definition: solvers.hpp:389
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1250
Stationary linear iteration: x &lt;- x + B (b - A x)
Definition: solvers.hpp:79
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:575
void SetAbsTol(double atol)
Definition: solvers.hpp:62
void SetRelTol(double rtol)
Definition: solvers.hpp:61
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:355
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:69
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:889
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1913
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:385
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:618
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1497
GMRES method.
Definition: solvers.hpp:143
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
Definition: solvers.cpp:505
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:125
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:205
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
void SetBounds(const Vector &_lo, const Vector &_hi)
Definition: solvers.cpp:1479
void UpdateVectors()
Definition: solvers.cpp:106
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:253
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:190
klu_common Common
Definition: solvers.hpp:410
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:284
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:476
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:204
virtual void Mult(const Vector &xt, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1510
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:125
const double alpha
Definition: ex15.cpp:339
Vector data type.
Definition: vector.hpp:48
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:88
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:1327
virtual ~KLUSolver()
Definition: solvers.cpp:1941
void UpdateVectors()
Definition: solvers.cpp:283
BiCGSTAB method.
Definition: solvers.hpp:190
SparseMatrix * mat
Definition: solvers.hpp:388
Base class for solvers.
Definition: operator.hpp:279
virtual void SetPreconditioner(Solver &pr)
These are not currently meaningful for this solver and will error out.
Definition: solvers.cpp:1491
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:64
Abstract operator.
Definition: operator.hpp:21
klu_numeric * Numeric
Definition: solvers.hpp:390
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1794
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:93
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:255
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:853
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1239
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:1697