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