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