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