MFEM  v4.1.0 Finite element discretization library
solvers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "linalg.hpp"
13 #include "../general/forall.hpp"
14 #include "../general/globals.hpp"
15 #include "../fem/bilinearform.hpp"
16 #include <iostream>
17 #include <iomanip>
18 #include <algorithm>
19 #include <cmath>
20 #include <set>
21
22 namespace mfem
23 {
24
25 using namespace std;
26
28  : Solver(0, true)
29 {
30  oper = NULL;
31  prec = NULL;
32  max_iter = 10;
33  print_level = -1;
34  rel_tol = abs_tol = 0.0;
35 #ifdef MFEM_USE_MPI
36  dot_prod_type = 0;
37 #endif
38 }
39
40 #ifdef MFEM_USE_MPI
42  : Solver(0, true)
43 {
44  oper = NULL;
45  prec = NULL;
46  max_iter = 10;
47  print_level = -1;
48  rel_tol = abs_tol = 0.0;
49  dot_prod_type = 1;
50  comm = _comm;
51 }
52 #endif
53
54 double IterativeSolver::Dot(const Vector &x, const Vector &y) const
55 {
56 #ifndef MFEM_USE_MPI
57  return (x * y);
58 #else
59  if (dot_prod_type == 0)
60  {
61  return (x * y);
62  }
63  else
64  {
65  return InnerProduct(comm, x, y);
66  }
67 #endif
68 }
69
70 void IterativeSolver::SetPrintLevel(int print_lvl)
71 {
72 #ifndef MFEM_USE_MPI
73  print_level = print_lvl;
74 #else
75  if (dot_prod_type == 0)
76  {
77  print_level = print_lvl;
78  }
79  else
80  {
81  int rank;
82  MPI_Comm_rank(comm, &rank);
83  if (rank == 0)
84  {
85  print_level = print_lvl;
86  }
87  }
88 #endif
89 }
90
92 {
93  prec = &pr;
94  prec->iterative_mode = false;
95 }
96
98 {
99  oper = &op;
100  height = op.Height();
101  width = op.Width();
102  if (prec)
103  {
104  prec->SetOperator(*oper);
105  }
106 }
107
108
110  const Array<int> &ess_tdofs,
111  const double dmpng)
112  :
113  Solver(a.FESpace()->GetTrueVSize()),
114  N(height),
115  dinv(N),
116  damping(dmpng),
117  ess_tdof_list(ess_tdofs),
118  residual(N)
119 {
120  Vector diag(N);
121  a.AssembleDiagonal(diag);
122  oper = &a;
123  Setup(diag);
124 }
125
127  const Array<int> &ess_tdofs,
128  const double dmpng)
129  :
130  Solver(d.Size()),
131  N(d.Size()),
132  dinv(N),
133  damping(dmpng),
134  ess_tdof_list(ess_tdofs),
135  residual(N)
136 {
137  Setup(d);
138 }
139
141 {
142  residual.UseDevice(true);
143  const double delta = damping;
145  auto DI = dinv.Write();
146  MFEM_FORALL(i, N, DI[i] = delta / D[i]; );
148  MFEM_FORALL(i, ess_tdof_list.Size(), DI[I[i]] = delta; );
149 }
150
152 {
153  MFEM_ASSERT(x.Size() == N, "invalid input vector");
154  MFEM_ASSERT(y.Size() == N, "invalid output vector");
155
156  if (iterative_mode && oper)
157  {
158  oper->Mult(y, residual); // r = A x
159  subtract(x, residual, residual); // r = b - A x
160  }
161  else
162  {
163  residual = x;
164  y.UseDevice(true);
165  y = 0.0;
166  }
170  MFEM_FORALL(i, N, Y[i] += DI[i] * R[i]; );
171 }
172
173
175 {
176  r.SetSize(width);
177  z.SetSize(width);
178 }
179
180 void SLISolver::Mult(const Vector &b, Vector &x) const
181 {
182  int i;
183
184  // Optimized preconditioned SLI with fixed number of iterations and given
185  // initial guess
186  if (!rel_tol && iterative_mode && prec)
187  {
188  for (i = 0; i < max_iter; i++)
189  {
190  oper->Mult(x, r); // r = A x
191  subtract(b, r, r); // r = b - A x
192  prec->Mult(r, z); // z = B r
193  add(x, 1.0, z, x); // x = x + B (b - A x)
194  }
195  converged = 1;
196  final_iter = i;
197  return;
198  }
199
200  // Optimized preconditioned SLI with fixed number of iterations and zero
201  // initial guess
202  if (!rel_tol && !iterative_mode && prec)
203  {
204  prec->Mult(b, x); // x = B b (initial guess 0)
205  for (i = 1; i < max_iter; i++)
206  {
207  oper->Mult(x, r); // r = A x
208  subtract(b, r, r); // r = b - A x
209  prec->Mult(r, z); // z = B r
210  add(x, 1.0, z, x); // x = x + B (b - A x)
211  }
212  converged = 1;
213  final_iter = i;
214  return;
215  }
216
217  // General version of SLI with a relative tolerance, optional preconditioner
218  // and optional initial guess
219  double r0, nom, nom0, nomold = 1, cf;
220
221  if (iterative_mode)
222  {
223  oper->Mult(x, r);
224  subtract(b, r, r); // r = b - A x
225  }
226  else
227  {
228  r = b;
229  x = 0.0;
230  }
231
232  if (prec)
233  {
234  prec->Mult(r, z); // z = B r
235  nom0 = nom = Dot(z, r);
236  }
237  else
238  {
239  nom0 = nom = Dot(r, r);
240  }
241
242  if (print_level == 1)
243  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
244  << nom << '\n';
245
246  r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
247  if (nom <= r0)
248  {
249  converged = 1;
250  final_iter = 0;
251  final_norm = sqrt(nom);
252  return;
253  }
254
255  // start iteration
256  converged = 0;
258  for (i = 1; true; )
259  {
260  if (prec) // x = x + B (b - A x)
261  {
263  }
264  else
265  {
267  }
268
269  oper->Mult(x, r);
270  subtract(b, r, r); // r = b - A x
271
272  if (prec)
273  {
274  prec->Mult(r, z); // z = B r
275  nom = Dot(z, r);
276  }
277  else
278  {
279  nom = Dot(r, r);
280  }
281
282  cf = sqrt(nom/nomold);
283  if (print_level == 1)
284  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
285  << nom << "\tConv. rate: " << cf << '\n';
286  nomold = nom;
287
288  if (nom < r0)
289  {
290  if (print_level == 2)
291  mfem::out << "Number of SLI iterations: " << i << '\n'
292  << "Conv. rate: " << cf << '\n';
293  else if (print_level == 3)
294  mfem::out << "(B r_0, r_0) = " << nom0 << '\n'
295  << "(B r_N, r_N) = " << nom << '\n'
296  << "Number of SLI iterations: " << i << '\n';
297  converged = 1;
298  final_iter = i;
299  break;
300  }
301
302  if (++i > max_iter)
303  {
304  break;
305  }
306  }
307
308  if (print_level >= 0 && !converged)
309  {
310  mfem::err << "SLI: No convergence!" << '\n';
311  mfem::out << "(B r_0, r_0) = " << nom0 << '\n'
312  << "(B r_N, r_N) = " << nom << '\n'
313  << "Number of SLI iterations: " << final_iter << '\n';
314  }
315  if (print_level >= 1 || (print_level >= 0 && !converged))
316  {
317  mfem::out << "Average reduction factor = "
318  << pow (nom/nom0, 0.5/final_iter) << '\n';
319  }
320  final_norm = sqrt(nom);
321 }
322
323 void SLI(const Operator &A, const Vector &b, Vector &x,
324  int print_iter, int max_num_iter,
325  double RTOLERANCE, double ATOLERANCE)
326 {
327  SLISolver sli;
328  sli.SetPrintLevel(print_iter);
329  sli.SetMaxIter(max_num_iter);
330  sli.SetRelTol(sqrt(RTOLERANCE));
331  sli.SetAbsTol(sqrt(ATOLERANCE));
332  sli.SetOperator(A);
333  sli.Mult(b, x);
334 }
335
336 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
337  int print_iter, int max_num_iter,
338  double RTOLERANCE, double ATOLERANCE)
339 {
340  SLISolver sli;
341  sli.SetPrintLevel(print_iter);
342  sli.SetMaxIter(max_num_iter);
343  sli.SetRelTol(sqrt(RTOLERANCE));
344  sli.SetAbsTol(sqrt(ATOLERANCE));
345  sli.SetOperator(A);
346  sli.SetPreconditioner(B);
347  sli.Mult(b, x);
348 }
349
350
352 {
353  r.SetSize(width);
354  d.SetSize(width);
355  z.SetSize(width);
356 }
357
358 void CGSolver::Mult(const Vector &b, Vector &x) const
359 {
360  int i;
361  double r0, den, nom, nom0, betanom, alpha, beta;
362
363  if (iterative_mode)
364  {
365  oper->Mult(x, r);
366  subtract(b, r, r); // r = b - A x
367  }
368  else
369  {
370  r = b;
371  x = 0.0;
372  }
373
374  if (prec)
375  {
376  prec->Mult(r, z); // z = B r
377  d = z;
378  }
379  else
380  {
381  d = r;
382  }
383  nom0 = nom = Dot(d, r);
384  MFEM_ASSERT(IsFinite(nom), "nom = " << nom);
385
386  if (print_level == 1 || print_level == 3)
387  {
388  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
389  << nom << (print_level == 3 ? " ...\n" : "\n");
390  }
391
392  r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
393  if (nom <= r0)
394  {
395  converged = 1;
396  final_iter = 0;
397  final_norm = sqrt(nom);
398  return;
399  }
400
401  oper->Mult(d, z); // z = A d
402  den = Dot(z, d);
403  MFEM_ASSERT(IsFinite(den), "den = " << den);
404  if (den <= 0.0)
405  {
406  if (Dot(d, d) > 0.0 && print_level >= 0)
407  {
408  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
409  << den << '\n';
410  }
411  if (den == 0.0)
412  {
413  converged = 0;
414  final_iter = 0;
415  final_norm = sqrt(nom);
416  return;
417  }
418  }
419
420  // start iteration
421  converged = 0;
423  for (i = 1; true; )
424  {
425  alpha = nom/den;
426  add(x, alpha, d, x); // x = x + alpha d
427  add(r, -alpha, z, r); // r = r - alpha A d
428
429  if (prec)
430  {
431  prec->Mult(r, z); // z = B r
432  betanom = Dot(r, z);
433  }
434  else
435  {
436  betanom = Dot(r, r);
437  }
438  MFEM_ASSERT(IsFinite(betanom), "betanom = " << betanom);
439
440  if (print_level == 1)
441  {
442  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
443  << betanom << '\n';
444  }
445
446  if (betanom < r0)
447  {
448  if (print_level == 2)
449  {
450  mfem::out << "Number of PCG iterations: " << i << '\n';
451  }
452  else if (print_level == 3)
453  {
454  mfem::out << " Iteration : " << setw(3) << i << " (B r, r) = "
455  << betanom << '\n';
456  }
457  converged = 1;
458  final_iter = i;
459  break;
460  }
461
462  if (++i > max_iter)
463  {
464  break;
465  }
466
467  beta = betanom/nom;
468  if (prec)
469  {
470  add(z, beta, d, d); // d = z + beta d
471  }
472  else
473  {
475  }
476  oper->Mult(d, z); // z = A d
477  den = Dot(d, z);
478  MFEM_ASSERT(IsFinite(den), "den = " << den);
479  if (den <= 0.0)
480  {
481  if (Dot(d, d) > 0.0 && print_level >= 0)
482  {
483  mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
484  << den << '\n';
485  }
486  if (den == 0.0)
487  {
488  final_iter = i;
489  break;
490  }
491  }
492  nom = betanom;
493  }
494  if (print_level >= 0 && !converged)
495  {
496  if (print_level != 1)
497  {
498  if (print_level != 3)
499  {
500  mfem::out << " Iteration : " << setw(3) << 0 << " (B r, r) = "
501  << nom0 << " ...\n";
502  }
503  mfem::out << " Iteration : " << setw(3) << final_iter << " (B r, r) = "
504  << betanom << '\n';
505  }
506  mfem::out << "PCG: No convergence!" << '\n';
507  }
508  if (print_level >= 1 || (print_level >= 0 && !converged))
509  {
510  mfem::out << "Average reduction factor = "
511  << pow (betanom/nom0, 0.5/final_iter) << '\n';
512  }
513  final_norm = sqrt(betanom);
514 }
515
516 void CG(const Operator &A, const Vector &b, Vector &x,
517  int print_iter, int max_num_iter,
518  double RTOLERANCE, double ATOLERANCE)
519 {
520  CGSolver cg;
521  cg.SetPrintLevel(print_iter);
522  cg.SetMaxIter(max_num_iter);
523  cg.SetRelTol(sqrt(RTOLERANCE));
524  cg.SetAbsTol(sqrt(ATOLERANCE));
525  cg.SetOperator(A);
526  cg.Mult(b, x);
527 }
528
529 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
530  int print_iter, int max_num_iter,
531  double RTOLERANCE, double ATOLERANCE)
532 {
533  CGSolver pcg;
534  pcg.SetPrintLevel(print_iter);
535  pcg.SetMaxIter(max_num_iter);
536  pcg.SetRelTol(sqrt(RTOLERANCE));
537  pcg.SetAbsTol(sqrt(ATOLERANCE));
538  pcg.SetOperator(A);
539  pcg.SetPreconditioner(B);
540  pcg.Mult(b, x);
541 }
542
543
544 inline void GeneratePlaneRotation(double &dx, double &dy,
545  double &cs, double &sn)
546 {
547  if (dy == 0.0)
548  {
549  cs = 1.0;
550  sn = 0.0;
551  }
552  else if (fabs(dy) > fabs(dx))
553  {
554  double temp = dx / dy;
555  sn = 1.0 / sqrt( 1.0 + temp*temp );
556  cs = temp * sn;
557  }
558  else
559  {
560  double temp = dy / dx;
561  cs = 1.0 / sqrt( 1.0 + temp*temp );
562  sn = temp * cs;
563  }
564 }
565
566 inline void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
567 {
568  double temp = cs * dx + sn * dy;
569  dy = -sn * dx + cs * dy;
570  dx = temp;
571 }
572
573 inline void Update(Vector &x, int k, DenseMatrix &h, Vector &s,
574  Array<Vector*> &v)
575 {
576  Vector y(s);
577
578  // Backsolve:
579  for (int i = k; i >= 0; i--)
580  {
581  y(i) /= h(i,i);
582  for (int j = i - 1; j >= 0; j--)
583  {
584  y(j) -= h(j,i) * y(i);
585  }
586  }
587
588  for (int j = 0; j <= k; j++)
589  {
591  }
592 }
593
594 void GMRESSolver::Mult(const Vector &b, Vector &x) const
595 {
596  // Generalized Minimum Residual method following the algorithm
597  // on p. 20 of the SIAM Templates book.
598
599  int n = width;
600
601  DenseMatrix H(m+1, m);
602  Vector s(m+1), cs(m+1), sn(m+1);
603  Vector r(n), w(n);
604  Array<Vector *> v;
605
606  double resid;
607  int i, j, k;
608
609  if (iterative_mode)
610  {
611  oper->Mult(x, r);
612  }
613  else
614  {
615  x = 0.0;
616  }
617
618  if (prec)
619  {
620  if (iterative_mode)
621  {
622  subtract(b, r, w);
623  prec->Mult(w, r); // r = M (b - A x)
624  }
625  else
626  {
627  prec->Mult(b, r);
628  }
629  }
630  else
631  {
632  if (iterative_mode)
633  {
634  subtract(b, r, r);
635  }
636  else
637  {
638  r = b;
639  }
640  }
641  double beta = Norm(r); // beta = ||r||
642  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
643
644  final_norm = std::max(rel_tol*beta, abs_tol);
645
646  if (beta <= final_norm)
647  {
648  final_norm = beta;
649  final_iter = 0;
650  converged = 1;
651  goto finish;
652  }
653
654  if (print_level == 1 || print_level == 3)
655  {
656  mfem::out << " Pass : " << setw(2) << 1
657  << " Iteration : " << setw(3) << 0
658  << " ||B r|| = " << beta << (print_level == 3 ? " ...\n" : "\n");
659  }
660
661  v.SetSize(m+1, NULL);
662
663  for (j = 1; j <= max_iter; )
664  {
665  if (v[0] == NULL) { v[0] = new Vector(n); }
666  v[0]->Set(1.0/beta, r);
667  s = 0.0; s(0) = beta;
668
669  for (i = 0; i < m && j <= max_iter; i++, j++)
670  {
671  if (prec)
672  {
673  oper->Mult(*v[i], r);
674  prec->Mult(r, w); // w = M A v[i]
675  }
676  else
677  {
678  oper->Mult(*v[i], w);
679  }
680
681  for (k = 0; k <= i; k++)
682  {
683  H(k,i) = Dot(w, *v[k]); // H(k,i) = w * v[k]
684  w.Add(-H(k,i), *v[k]); // w -= H(k,i) * v[k]
685  }
686
687  H(i+1,i) = Norm(w); // H(i+1,i) = ||w||
688  MFEM_ASSERT(IsFinite(H(i+1,i)), "Norm(w) = " << H(i+1,i));
689  if (v[i+1] == NULL) { v[i+1] = new Vector(n); }
690  v[i+1]->Set(1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
691
692  for (k = 0; k < i; k++)
693  {
694  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
695  }
696
697  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
698  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
699  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
700
701  resid = fabs(s(i+1));
702  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
703
704  if (resid <= final_norm)
705  {
706  Update(x, i, H, s, v);
707  final_norm = resid;
708  final_iter = j;
709  converged = 1;
710  goto finish;
711  }
712
713  if (print_level == 1)
714  {
715  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
716  << " Iteration : " << setw(3) << j
717  << " ||B r|| = " << resid << '\n';
718  }
719  }
720
721  if (print_level == 1 && j <= max_iter)
722  {
723  mfem::out << "Restarting..." << '\n';
724  }
725
726  Update(x, i-1, H, s, v);
727
728  oper->Mult(x, r);
729  if (prec)
730  {
731  subtract(b, r, w);
732  prec->Mult(w, r); // r = M (b - A x)
733  }
734  else
735  {
736  subtract(b, r, r);
737  }
738  beta = Norm(r); // beta = ||r||
739  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
740  if (beta <= final_norm)
741  {
742  final_norm = beta;
743  final_iter = j;
744  converged = 1;
745  goto finish;
746  }
747  }
748
749  final_norm = beta;
751  converged = 0;
752
753 finish:
754  if (print_level == 1 || print_level == 3)
755  {
756  mfem::out << " Pass : " << setw(2) << (final_iter-1)/m+1
757  << " Iteration : " << setw(3) << final_iter
758  << " ||B r|| = " << final_norm << '\n';
759  }
760  else if (print_level == 2)
761  {
762  mfem::out << "GMRES: Number of iterations: " << final_iter << '\n';
763  }
764  if (print_level >= 0 && !converged)
765  {
766  mfem::out << "GMRES: No convergence!\n";
767  }
768  for (i = 0; i < v.Size(); i++)
769  {
770  delete v[i];
771  }
772 }
773
774 void FGMRESSolver::Mult(const Vector &b, Vector &x) const
775 {
776  DenseMatrix H(m+1,m);
777  Vector s(m+1), cs(m+1), sn(m+1);
778  Vector r(b.Size());
779
780  int i, j, k;
781
782
783  if (iterative_mode)
784  {
785  oper->Mult(x, r);
786  subtract(b,r,r);
787  }
788  else
789  {
790  x = 0.;
791  r = b;
792  }
793  double beta = Norm(r); // beta = ||r||
794  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
795
796  final_norm = std::max(rel_tol*beta, abs_tol);
797
798  if (beta <= final_norm)
799  {
800  final_norm = beta;
801  final_iter = 0;
802  converged = 1;
803  return;
804  }
805
806  if (print_level == 1)
807  {
808  mfem::out << " Pass : " << setw(2) << 1
809  << " Iteration : " << setw(3) << 0
810  << " || r || = " << beta << endl;
811  }
812
813  Array<Vector*> v(m+1);
814  Array<Vector*> z(m+1);
815  for (i= 0; i<=m; i++)
816  {
817  v[i] = NULL;
818  z[i] = NULL;
819  }
820
821  j = 1;
822  while (j <= max_iter)
823  {
824  if (v[0] == NULL) { v[0] = new Vector(b.Size()); }
825  (*v[0]) = 0.0;
826  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
827  s = 0.0; s(0) = beta;
828
829  for (i = 0; i < m && j <= max_iter; i++, j++)
830  {
831
832  if (z[i] == NULL) { z[i] = new Vector(b.Size()); }
833  (*z[i]) = 0.0;
834
835  if (prec)
836  {
837  prec->Mult(*v[i], *z[i]);
838  }
839  else
840  {
841  (*z[i]) = (*v[i]);
842  }
843  oper->Mult(*z[i], r);
844
845  for (k = 0; k <= i; k++)
846  {
847  H(k,i) = Dot( r, *v[k]); // H(k,i) = r * v[k]
848  r.Add(-H(k,i), (*v[k])); // r -= H(k,i) * v[k]
849  }
850
851  H(i+1,i) = Norm(r); // H(i+1,i) = ||r||
852  if (v[i+1] == NULL) { v[i+1] = new Vector(b.Size()); }
853  (*v[i+1]) = 0.0;
854  v[i+1] -> Add (1.0/H(i+1,i), r); // v[i+1] = r / H(i+1,i)
855
856  for (k = 0; k < i; k++)
857  {
858  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
859  }
860
861  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
862  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
863  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
864
865  double resid = fabs(s(i+1));
866  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
867  if (print_level == 1)
868  {
869  mfem::out << " Pass : " << setw(2) << (j-1)/m+1
870  << " Iteration : " << setw(3) << j
871  << " || r || = " << resid << endl;
872  }
873
874  if (resid <= final_norm)
875  {
876  Update(x, i, H, s, z);
877  final_norm = resid;
878  final_iter = j;
879  converged = 1;
880
881  if (print_level == 2)
882  {
883  mfem::out << "Number of FGMRES iterations: " << final_iter << endl;
884  }
885
886  for (i= 0; i<=m; i++)
887  {
888  if (v[i]) { delete v[i]; }
889  if (z[i]) { delete z[i]; }
890  }
891  return;
892  }
893  }
894
895  if (print_level == 1)
896  {
897  mfem::out << "Restarting..." << endl;
898  }
899
900  Update(x, i-1, H, s, z);
901
902  oper->Mult(x, r);
903  subtract(b,r,r);
904  beta = Norm(r);
905  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
906  if (beta <= final_norm)
907  {
908  final_norm = beta;
909  final_iter = j;
910  converged = 1;
911
912  if (print_level == 2)
913  {
914  mfem::out << "Number of FGMRES iterations: " << final_iter << endl;
915  }
916
917  for (i= 0; i<=m; i++)
918  {
919  if (v[i]) { delete v[i]; }
920  if (z[i]) { delete z[i]; }
921  }
922  return;
923  }
924  }
925
926  for (i = 0; i <= m; i++)
927  {
928  if (v[i]) { delete v[i]; }
929  if (z[i]) { delete z[i]; }
930  }
931  converged = 0;
932
933  if (print_level >= 0)
934  {
935  mfem::out << "FGMRES: No convergence!" << endl;
936  }
937
938  return;
939 }
940
941
942 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
943  int &max_iter, int m, double &tol, double atol, int printit)
944 {
945  GMRESSolver gmres;
946  gmres.SetPrintLevel(printit);
947  gmres.SetMaxIter(max_iter);
948  gmres.SetKDim(m);
949  gmres.SetRelTol(sqrt(tol));
950  gmres.SetAbsTol(sqrt(atol));
951  gmres.SetOperator(A);
952  gmres.SetPreconditioner(M);
953  gmres.Mult(b, x);
954  max_iter = gmres.GetNumIterations();
955  tol = gmres.GetFinalNorm()*gmres.GetFinalNorm();
956  return gmres.GetConverged();
957 }
958
959 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
960  int print_iter, int max_num_iter, int m, double rtol, double atol)
961 {
962  GMRES(A, x, b, B, max_num_iter, m, rtol, atol, print_iter);
963 }
964
965
967 {
968  p.SetSize(width);
969  phat.SetSize(width);
970  s.SetSize(width);
971  shat.SetSize(width);
972  t.SetSize(width);
973  v.SetSize(width);
974  r.SetSize(width);
976 }
977
978 void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
979 {
980  // BiConjugate Gradient Stabilized method following the algorithm
981  // on p. 27 of the SIAM Templates book.
982
983  int i;
984  double resid, tol_goal;
985  double rho_1, rho_2=1.0, alpha=1.0, beta, omega=1.0;
986
987  if (iterative_mode)
988  {
989  oper->Mult(x, r);
990  subtract(b, r, r); // r = b - A x
991  }
992  else
993  {
994  x = 0.0;
995  r = b;
996  }
997  rtilde = r;
998
999  resid = Norm(r);
1000  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1001  if (print_level >= 0)
1002  mfem::out << " Iteration : " << setw(3) << 0
1003  << " ||r|| = " << resid << '\n';
1004
1005  tol_goal = std::max(resid*rel_tol, abs_tol);
1006
1007  if (resid <= tol_goal)
1008  {
1009  final_norm = resid;
1010  final_iter = 0;
1011  converged = 1;
1012  return;
1013  }
1014
1015  for (i = 1; i <= max_iter; i++)
1016  {
1017  rho_1 = Dot(rtilde, r);
1018  if (rho_1 == 0)
1019  {
1020  if (print_level >= 0)
1021  mfem::out << " Iteration : " << setw(3) << i
1022  << " ||r|| = " << resid << '\n';
1023  final_norm = resid;
1024  final_iter = i;
1025  converged = 0;
1026  return;
1027  }
1028  if (i == 1)
1029  {
1030  p = r;
1031  }
1032  else
1033  {
1034  beta = (rho_1/rho_2) * (alpha/omega);
1035  add(p, -omega, v, p); // p = p - omega * v
1036  add(r, beta, p, p); // p = r + beta * p
1037  }
1038  if (prec)
1039  {
1040  prec->Mult(p, phat); // phat = M^{-1} * p
1041  }
1042  else
1043  {
1044  phat = p;
1045  }
1046  oper->Mult(phat, v); // v = A * phat
1047  alpha = rho_1 / Dot(rtilde, v);
1048  add(r, -alpha, v, s); // s = r - alpha * v
1049  resid = Norm(s);
1050  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1051  if (resid < tol_goal)
1052  {
1053  x.Add(alpha, phat); // x = x + alpha * phat
1054  if (print_level >= 0)
1055  mfem::out << " Iteration : " << setw(3) << i
1056  << " ||s|| = " << resid << '\n';
1057  final_norm = resid;
1058  final_iter = i;
1059  converged = 1;
1060  return;
1061  }
1062  if (print_level >= 0)
1063  mfem::out << " Iteration : " << setw(3) << i
1064  << " ||s|| = " << resid;
1065  if (prec)
1066  {
1067  prec->Mult(s, shat); // shat = M^{-1} * s
1068  }
1069  else
1070  {
1071  shat = s;
1072  }
1073  oper->Mult(shat, t); // t = A * shat
1074  omega = Dot(t, s) / Dot(t, t);
1075  x.Add(alpha, phat); // x += alpha * phat
1076  x.Add(omega, shat); // x += omega * shat
1077  add(s, -omega, t, r); // r = s - omega * t
1078
1079  rho_2 = rho_1;
1080  resid = Norm(r);
1081  MFEM_ASSERT(IsFinite(resid), "resid = " << resid);
1082  if (print_level >= 0)
1083  {
1084  mfem::out << " ||r|| = " << resid << '\n';
1085  }
1086  if (resid < tol_goal)
1087  {
1088  final_norm = resid;
1089  final_iter = i;
1090  converged = 1;
1091  return;
1092  }
1093  if (omega == 0)
1094  {
1095  final_norm = resid;
1096  final_iter = i;
1097  converged = 0;
1098  return;
1099  }
1100  }
1101
1102  final_norm = resid;
1103  final_iter = max_iter;
1104  converged = 0;
1105 }
1106
1107 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
1108  int &max_iter, double &tol, double atol, int printit)
1109 {
1110  BiCGSTABSolver bicgstab;
1111  bicgstab.SetPrintLevel(printit);
1112  bicgstab.SetMaxIter(max_iter);
1113  bicgstab.SetRelTol(sqrt(tol));
1114  bicgstab.SetAbsTol(sqrt(atol));
1115  bicgstab.SetOperator(A);
1116  bicgstab.SetPreconditioner(M);
1117  bicgstab.Mult(b, x);
1118  max_iter = bicgstab.GetNumIterations();
1119  tol = bicgstab.GetFinalNorm()*bicgstab.GetFinalNorm();
1120  return bicgstab.GetConverged();
1121 }
1122
1123 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
1124  int print_iter, int max_num_iter, double rtol, double atol)
1125 {
1126  BiCGSTAB(A, x, b, B, max_num_iter, rtol, atol, print_iter);
1127 }
1128
1129
1131 {
1133  v0.SetSize(width);
1134  v1.SetSize(width);
1135  w0.SetSize(width);
1136  w1.SetSize(width);
1137  q.SetSize(width);
1138  if (prec)
1139  {
1140  u1.SetSize(width);
1141  }
1142 }
1143
1144 void MINRESSolver::Mult(const Vector &b, Vector &x) const
1145 {
1146  // Based on the MINRES algorithm on p. 86, Fig. 6.9 in
1147  // "Iterative Krylov Methods for Large Linear Systems",
1148  // by Henk A. van der Vorst, 2003.
1149  // Extended to support an SPD preconditioner.
1150
1151  int it;
1152  double beta, eta, gamma0, gamma1, sigma0, sigma1;
1153  double alpha, delta, rho1, rho2, rho3, norm_goal;
1154  Vector *z = (prec) ? &u1 : &v1;
1155
1156  converged = 1;
1157
1158  if (!iterative_mode)
1159  {
1160  v1 = b;
1161  x = 0.;
1162  }
1163  else
1164  {
1165  oper->Mult(x, v1);
1166  subtract(b, v1, v1);
1167  }
1168
1169  if (prec)
1170  {
1171  prec->Mult(v1, u1);
1172  }
1173  eta = beta = sqrt(Dot(*z, v1));
1174  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1175  gamma0 = gamma1 = 1.;
1176  sigma0 = sigma1 = 0.;
1177
1178  norm_goal = std::max(rel_tol*eta, abs_tol);
1179
1180  if (eta <= norm_goal)
1181  {
1182  it = 0;
1183  goto loop_end;
1184  }
1185
1186  if (print_level == 1 || print_level == 3)
1187  {
1188  mfem::out << "MINRES: iteration " << setw(3) << 0 << ": ||r||_B = "
1189  << eta << (print_level == 3 ? " ...\n" : "\n");
1190  }
1191
1192  for (it = 1; it <= max_iter; it++)
1193  {
1194  v1 /= beta;
1195  if (prec)
1196  {
1197  u1 /= beta;
1198  }
1199  oper->Mult(*z, q);
1200  alpha = Dot(*z, q);
1201  MFEM_ASSERT(IsFinite(alpha), "alpha = " << alpha);
1202  if (it > 1) // (v0 == 0) for (it == 1)
1203  {
1205  }
1207
1208  delta = gamma1*alpha - gamma0*sigma1*beta;
1209  rho3 = sigma0*beta;
1210  rho2 = sigma1*alpha + gamma0*gamma1*beta;
1211  if (!prec)
1212  {
1213  beta = Norm(v0);
1214  }
1215  else
1216  {
1217  prec->Mult(v0, q);
1218  beta = sqrt(Dot(v0, q));
1219  }
1220  MFEM_ASSERT(IsFinite(beta), "beta = " << beta);
1221  rho1 = hypot(delta, beta);
1222
1223  if (it == 1)
1224  {
1225  w0.Set(1./rho1, *z); // (w0 == 0) and (w1 == 0)
1226  }
1227  else if (it == 2)
1228  {
1229  add(1./rho1, *z, -rho2/rho1, w1, w0); // (w0 == 0)
1230  }
1231  else
1232  {
1233  add(-rho3/rho1, w0, -rho2/rho1, w1, w0);
1235  }
1236
1237  gamma0 = gamma1;
1238  gamma1 = delta/rho1;
1239
1241
1242  sigma0 = sigma1;
1243  sigma1 = beta/rho1;
1244
1245  eta = -sigma1*eta;
1246  MFEM_ASSERT(IsFinite(eta), "eta = " << eta);
1247
1248  if (fabs(eta) <= norm_goal)
1249  {
1250  goto loop_end;
1251  }
1252
1253  if (print_level == 1)
1254  {
1255  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1256  << fabs(eta) << '\n';
1257  }
1258
1259  if (prec)
1260  {
1261  Swap(u1, q);
1262  }
1263  Swap(v0, v1);
1264  Swap(w0, w1);
1265  }
1266  converged = 0;
1267  it--;
1268
1269 loop_end:
1270  final_iter = it;
1271  final_norm = fabs(eta);
1272
1273  if (print_level == 1 || print_level == 3)
1274  {
1275  mfem::out << "MINRES: iteration " << setw(3) << final_iter << ": ||r||_B = "
1276  << final_norm << '\n';
1277  }
1278  else if (print_level == 2)
1279  {
1280  mfem::out << "MINRES: number of iterations: " << final_iter << '\n';
1281  }
1282 #if 0
1283  if (print_level >= 1)
1284  {
1285  oper->Mult(x, v1);
1286  subtract(b, v1, v1);
1287  if (prec)
1288  {
1289  prec->Mult(v1, u1);
1290  }
1291  eta = sqrt(Dot(*z, v1));
1292  mfem::out << "MINRES: iteration " << setw(3) << it << ": ||r||_B = "
1293  << eta << " (re-computed)" << '\n';
1294  }
1295 #endif
1296  if (!converged && print_level >= 0)
1297  {
1298  mfem::out << "MINRES: No convergence!\n";
1299  }
1300 }
1301
1302 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it,
1303  int max_it, double rtol, double atol)
1304 {
1305  MINRESSolver minres;
1306  minres.SetPrintLevel(print_it);
1307  minres.SetMaxIter(max_it);
1308  minres.SetRelTol(sqrt(rtol));
1309  minres.SetAbsTol(sqrt(atol));
1310  minres.SetOperator(A);
1311  minres.Mult(b, x);
1312 }
1313
1314 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
1315  int print_it, int max_it, double rtol, double atol)
1316 {
1317  MINRESSolver minres;
1318  minres.SetPrintLevel(print_it);
1319  minres.SetMaxIter(max_it);
1320  minres.SetRelTol(sqrt(rtol));
1321  minres.SetAbsTol(sqrt(atol));
1322  minres.SetOperator(A);
1323  minres.SetPreconditioner(B);
1324  minres.Mult(b, x);
1325 }
1326
1327
1329 {
1330  oper = &op;
1331  height = op.Height();
1332  width = op.Width();
1333  MFEM_ASSERT(height == width, "square Operator is required.");
1334
1335  r.SetSize(width);
1336  c.SetSize(width);
1337 }
1338
1339 void NewtonSolver::Mult(const Vector &b, Vector &x) const
1340 {
1341  MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
1342  MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver).");
1343
1344  int it;
1345  double norm0, norm, norm_goal;
1346  const bool have_b = (b.Size() == Height());
1347
1348  if (!iterative_mode)
1349  {
1350  x = 0.0;
1351  }
1352
1353  oper->Mult(x, r);
1354  if (have_b)
1355  {
1356  r -= b;
1357  }
1358
1359  norm0 = norm = Norm(r);
1360  norm_goal = std::max(rel_tol*norm, abs_tol);
1361
1362  prec->iterative_mode = false;
1363
1364  // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b]
1365  for (it = 0; true; it++)
1366  {
1367  MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
1368  if (print_level >= 0)
1369  {
1370  mfem::out << "Newton iteration " << setw(2) << it
1371  << " : ||r|| = " << norm;
1372  if (it > 0)
1373  {
1374  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
1375  }
1376  mfem::out << '\n';
1377  }
1378
1379  if (norm <= norm_goal)
1380  {
1381  converged = 1;
1382  break;
1383  }
1384
1385  if (it >= max_iter)
1386  {
1387  converged = 0;
1388  break;
1389  }
1390
1392
1393  prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b]
1394
1395  const double c_scale = ComputeScalingFactor(x, b);
1396  if (c_scale == 0.0)
1397  {
1398  converged = 0;
1399  break;
1400  }
1402
1403  ProcessNewState(x);
1404
1405  oper->Mult(x, r);
1406  if (have_b)
1407  {
1408  r -= b;
1409  }
1410  norm = Norm(r);
1411  }
1412
1413  final_iter = it;
1414  final_norm = norm;
1415 }
1416
1417
1418 int aGMRES(const Operator &A, Vector &x, const Vector &b,
1419  const Operator &M, int &max_iter,
1420  int m_max, int m_min, int m_step, double cf,
1421  double &tol, double &atol, int printit)
1422 {
1423  int n = A.Width();
1424
1425  int m = m_max;
1426
1427  DenseMatrix H(m+1,m);
1428  Vector s(m+1), cs(m+1), sn(m+1);
1429  Vector w(n), av(n);
1430
1431  double r1, resid;
1432  int i, j, k;
1433
1434  M.Mult(b,w);
1435  double normb = w.Norml2(); // normb = ||M b||
1436  if (normb == 0.0)
1437  {
1438  normb = 1;
1439  }
1440
1441  Vector r(n);
1442  A.Mult(x, r);
1443  subtract(b,r,w);
1444  M.Mult(w, r); // r = M (b - A x)
1445  double beta = r.Norml2(); // beta = ||r||
1446
1447  resid = beta / normb;
1448
1449  if (resid * resid <= tol)
1450  {
1451  tol = resid * resid;
1452  max_iter = 0;
1453  return 0;
1454  }
1455
1456  if (printit)
1457  mfem::out << " Pass : " << setw(2) << 1
1458  << " Iteration : " << setw(3) << 0
1459  << " (r, r) = " << beta*beta << '\n';
1460
1461  tol *= (normb*normb);
1462  tol = (atol > tol) ? atol : tol;
1463
1464  m = m_max;
1465  Array<Vector *> v(m+1);
1466  for (i= 0; i<=m; i++)
1467  {
1468  v[i] = new Vector(n);
1469  (*v[i]) = 0.0;
1470  }
1471
1472  j = 1;
1473  while (j <= max_iter)
1474  {
1475  (*v[0]) = 0.0;
1476  v[0] -> Add (1.0/beta, r); // v[0] = r / ||r||
1477  s = 0.0; s(0) = beta;
1478
1479  r1 = beta;
1480
1481  for (i = 0; i < m && j <= max_iter; i++)
1482  {
1483  A.Mult((*v[i]),av);
1484  M.Mult(av,w); // w = M A v[i]
1485
1486  for (k = 0; k <= i; k++)
1487  {
1488  H(k,i) = w * (*v[k]); // H(k,i) = w * v[k]
1489  w.Add(-H(k,i), (*v[k])); // w -= H(k,i) * v[k]
1490  }
1491
1492  H(i+1,i) = w.Norml2(); // H(i+1,i) = ||w||
1493  (*v[i+1]) = 0.0;
1494  v[i+1] -> Add (1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
1495
1496  for (k = 0; k < i; k++)
1497  {
1498  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
1499  }
1500
1501  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1502  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
1503  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
1504
1505  resid = fabs(s(i+1));
1506  if (printit)
1507  mfem::out << " Pass : " << setw(2) << j
1508  << " Iteration : " << setw(3) << i+1
1509  << " (r, r) = " << resid*resid << '\n';
1510
1511  if ( resid*resid < tol)
1512  {
1513  Update(x, i, H, s, v);
1514  tol = resid * resid;
1515  max_iter = j;
1516  for (i= 0; i<=m; i++)
1517  {
1518  delete v[i];
1519  }
1520  return 0;
1521  }
1522  }
1523
1524  if (printit)
1525  {
1526  mfem::out << "Restarting..." << '\n';
1527  }
1528
1529  Update(x, i-1, H, s, v);
1530
1531  A.Mult(x, r);
1532  subtract(b,r,w);
1533  M.Mult(w, r); // r = M (b - A x)
1534  beta = r.Norml2(); // beta = ||r||
1535  if ( resid*resid < tol)
1536  {
1537  tol = resid * resid;
1538  max_iter = j;
1539  for (i= 0; i<=m; i++)
1540  {
1541  delete v[i];
1542  }
1543  return 0;
1544  }
1545
1546  if (beta/r1 > cf)
1547  {
1548  if (m - m_step >= m_min)
1549  {
1550  m -= m_step;
1551  }
1552  else
1553  {
1554  m = m_max;
1555  }
1556  }
1557
1558  j++;
1559  }
1560
1561  tol = resid * resid;
1562  for (i= 0; i<=m; i++)
1563  {
1564  delete v[i];
1565  }
1566  return 1;
1567 }
1568
1570  const Operator *C_,
1571  const Operator *D_)
1572  : C(C_), D(D_), c_e(NULL), d_lo(NULL), d_hi(NULL), x_lo(NULL), x_hi(NULL),
1573  input_size(insize)
1574 {
1575  if (C) { MFEM_ASSERT(C->Width() == input_size, "Wrong width of C."); }
1576  if (D) { MFEM_ASSERT(D->Width() == input_size, "Wrong width of D."); }
1577 }
1578
1580 {
1581  MFEM_ASSERT(C, "The C operator is unspecified -- can't set constraints.");
1582  MFEM_ASSERT(c.Size() == C->Height(), "Wrong size of the constraint.");
1583
1584  c_e = &c;
1585 }
1586
1588  const Vector &dh)
1589 {
1590  MFEM_ASSERT(D, "The D operator is unspecified -- can't set constraints.");
1591  MFEM_ASSERT(dl.Size() == D->Height() && dh.Size() == D->Height(),
1592  "Wrong size of the constraint.");
1593
1594  d_lo = &dl; d_hi = &dh;
1595 }
1596
1598 {
1599  MFEM_ASSERT(xl.Size() == input_size && xh.Size() == input_size,
1600  "Wrong size of the constraint.");
1601
1602  x_lo = &xl; x_hi = &xh;
1603 }
1604
1606 {
1607  int m = 0;
1608  if (C) { m += C->Height(); }
1609  if (D) { m += D->Height(); }
1610  return m;
1611 }
1612
1614 {
1615  if (print_level > 1)
1616  {
1617  MFEM_WARNING("Objective functional is ignored as SLBQP always minimizes"
1618  "the l2 norm of (x - x_target).");
1619  }
1620  MFEM_ASSERT(prob.GetC(), "Linear constraint is not set.");
1621  MFEM_ASSERT(prob.GetC()->Height() == 1, "Solver expects scalar constraint.");
1622
1623  problem = &prob;
1624 }
1625
1626 void SLBQPOptimizer::SetBounds(const Vector &_lo, const Vector &_hi)
1627 {
1628  lo.SetDataAndSize(_lo.GetData(), _lo.Size());
1629  hi.SetDataAndSize(_hi.GetData(), _hi.Size());
1630 }
1631
1633 {
1634  w.SetDataAndSize(_w.GetData(), _w.Size());
1635  a = _a;
1636 }
1637
1638 inline void SLBQPOptimizer::print_iteration(int it, double r, double l) const
1639 {
1640  if (print_level > 1)
1641  mfem::out << "SLBQP iteration " << it << ": residual = " << r
1642  << ", lambda = " << l << '\n';
1643 }
1644
1645 void SLBQPOptimizer::Mult(const Vector& xt, Vector& x) const
1646 {
1647  // Based on code provided by Denis Ridzal, dridzal@sandia.gov.
1648  // Algorithm adapted from Dai and Fletcher, "New Algorithms for
1649  // Singly Linearly Constrained Quadratic Programs Subject to Lower
1650  // and Upper Bounds", Numerical Analysis Report NA/216, 2003.
1651
1652  // Set some algorithm-specific constants and temporaries.
1653  int nclip = 0;
1654  double l = 0;
1655  double llow = 0;
1656  double lupp = 0;
1657  double lnew = 0;
1658  double dl = 2;
1659  double r = 0;
1660  double rlow = 0;
1661  double rupp = 0;
1662  double s = 0;
1663
1664  const double smin = 0.1;
1665
1666  const double tol = max(abs_tol, rel_tol*a);
1667
1668  // *** Start bracketing phase of SLBQP ***
1669  if (print_level > 1)
1670  {
1671  mfem::out << "SLBQP bracketing phase" << '\n';
1672  }
1673
1674  // Solve QP with fixed Lagrange multiplier
1675  r = solve(l,xt,x,nclip);
1676  print_iteration(nclip, r, l);
1677
1678
1679  // If x=xt was already within bounds and satisfies the linear
1680  // constraint, then we already have the solution.
1681  if (fabs(r) <= tol)
1682  {
1683  converged = true;
1684  goto slbqp_done;
1685  }
1686
1687  if (r < 0)
1688  {
1689  llow = l; rlow = r; l = l + dl;
1690
1691  // Solve QP with fixed Lagrange multiplier
1692  r = solve(l,xt,x,nclip);
1693  print_iteration(nclip, r, l);
1694
1695  while ((r < 0) && (nclip < max_iter))
1696  {
1697  llow = l;
1698  s = rlow/r - 1.0;
1699  if (s < smin) { s = smin; }
1700  dl = dl + dl/s;
1701  l = l + dl;
1702
1703  // Solve QP with fixed Lagrange multiplier
1704  r = solve(l,xt,x,nclip);
1705  print_iteration(nclip, r, l);
1706  }
1707
1708  lupp = l; rupp = r;
1709  }
1710  else
1711  {
1712  lupp = l; rupp = r; l = l - dl;
1713
1714  // Solve QP with fixed Lagrange multiplier
1715  r = solve(l,xt,x,nclip);
1716  print_iteration(nclip, r, l);
1717
1718  while ((r > 0) && (nclip < max_iter))
1719  {
1720  lupp = l;
1721  s = rupp/r - 1.0;
1722  if (s < smin) { s = smin; }
1723  dl = dl + dl/s;
1724  l = l - dl;
1725
1726  // Solve QP with fixed Lagrange multiplier
1727  r = solve(l,xt,x,nclip);
1728  print_iteration(nclip, r, l);
1729  }
1730
1731  llow = l; rlow = r;
1732  }
1733
1734  // *** Stop bracketing phase of SLBQP ***
1735
1736
1737  // *** Start secant phase of SLBQP ***
1738  if (print_level > 1)
1739  {
1740  mfem::out << "SLBQP secant phase" << '\n';
1741  }
1742
1743  s = 1.0 - rlow/rupp; dl = dl/s; l = lupp - dl;
1744
1745  // Solve QP with fixed Lagrange multiplier
1746  r = solve(l,xt,x,nclip);
1747  print_iteration(nclip, r, l);
1748
1749  while ( (fabs(r) > tol) && (nclip < max_iter) )
1750  {
1751  if (r > 0)
1752  {
1753  if (s <= 2.0)
1754  {
1755  lupp = l; rupp = r; s = 1.0 - rlow/rupp;
1756  dl = (lupp - llow)/s; l = lupp - dl;
1757  }
1758  else
1759  {
1760  s = rupp/r - 1.0;
1761  if (s < smin) { s = smin; }
1762  dl = (lupp - l)/s;
1763  lnew = 0.75*llow + 0.25*l;
1764  if (lnew < l-dl) { lnew = l-dl; }
1765  lupp = l; rupp = r; l = lnew;
1766  s = (lupp - llow)/(lupp - l);
1767  }
1768
1769  }
1770  else
1771  {
1772  if (s >= 2.0)
1773  {
1774  llow = l; rlow = r; s = 1.0 - rlow/rupp;
1775  dl = (lupp - llow)/s; l = lupp - dl;
1776  }
1777  else
1778  {
1779  s = rlow/r - 1.0;
1780  if (s < smin) { s = smin; }
1781  dl = (l - llow)/s;
1782  lnew = 0.75*lupp + 0.25*l;
1783  if (lnew < l+dl) { lnew = l+dl; }
1784  llow = l; rlow = r; l = lnew;
1785  s = (lupp - llow)/(lupp - l);
1786  }
1787  }
1788
1789  // Solve QP with fixed Lagrange multiplier
1790  r = solve(l,xt,x,nclip);
1791  print_iteration(nclip, r, l);
1792  }
1793
1794  // *** Stop secant phase of SLBQP ***
1795
1796  converged = (fabs(r) <= tol);
1797  if (!converged && print_level >= 0)
1798  {
1799  mfem::err << "SLBQP not converged!" << '\n';
1800  }
1801
1802 slbqp_done:
1803
1804  final_iter = nclip;
1805  final_norm = r;
1806
1807  if (print_level == 1 || (!converged && print_level >= 0))
1808  {
1809  mfem::out << "SLBQP iterations = " << nclip << '\n';
1810  mfem::out << "SLBQP lambda = " << l << '\n';
1811  mfem::out << "SLBQP residual = " << r << '\n';
1812  }
1813 }
1814
1815 struct WeightMinHeap
1816 {
1817  const std::vector<double> &w;
1818  std::vector<size_t> c;
1819  std::vector<int> loc;
1820
1821  WeightMinHeap(const std::vector<double> &w_) : w(w_)
1822  {
1823  c.reserve(w.size());
1824  loc.resize(w.size());
1825  for (size_t i=0; i<w.size(); ++i) { push(i); }
1826  }
1827
1828  size_t percolate_up(size_t pos, double val)
1829  {
1830  for (; pos > 0 && w[c[(pos-1)/2]] > val; pos = (pos-1)/2)
1831  {
1832  c[pos] = c[(pos-1)/2];
1833  loc[c[(pos-1)/2]] = pos;
1834  }
1835  return pos;
1836  }
1837
1838  size_t percolate_down(size_t pos, double val)
1839  {
1840  while (2*pos+1 < c.size())
1841  {
1842  size_t left = 2*pos+1;
1843  size_t right = left+1;
1844  size_t tgt;
1845  if (right < c.size() && w[c[right]] < w[c[left]]) { tgt = right; }
1846  else { tgt = left; }
1847  if (w[c[tgt]] < val)
1848  {
1849  c[pos] = c[tgt];
1850  loc[c[tgt]] = pos;
1851  pos = tgt;
1852  }
1853  else
1854  {
1855  break;
1856  }
1857  }
1858  return pos;
1859  }
1860
1861  void push(size_t i)
1862  {
1863  double val = w[i];
1864  c.push_back(0);
1865  size_t pos = c.size()-1;
1866  pos = percolate_up(pos, val);
1867  c[pos] = i;
1868  loc[i] = pos;
1869  }
1870
1871  int pop()
1872  {
1873  size_t i = c[0];
1874  size_t j = c.back();
1875  c.pop_back();
1876  // Mark as removed
1877  loc[i] = -1;
1878  if (c.empty()) { return i; }
1879  double val = w[j];
1880  size_t pos = 0;
1881  pos = percolate_down(pos, val);
1882  c[pos] = j;
1883  loc[j] = pos;
1884  return i;
1885  }
1886
1887  void update(size_t i)
1888  {
1889  size_t pos = loc[i];
1890  double val = w[i];
1891  pos = percolate_up(pos, val);
1892  pos = percolate_down(pos, val);
1893  c[pos] = i;
1894  loc[i] = pos;
1895  }
1896
1897  bool picked(size_t i)
1898  {
1899  return loc[i] < 0;
1900  }
1901 };
1902
1904 {
1905  int n = C.Width();
1906  // Scale rows by reciprocal of diagonal and take absolute value
1907  Vector D;
1908  C.GetDiag(D);
1909  int *I = C.GetI();
1910  int *J = C.GetJ();
1911  double *V = C.GetData();
1912  for (int i=0; i<n; ++i)
1913  {
1914  for (int j=I[i]; j<I[i+1]; ++j)
1915  {
1916  V[j] = abs(V[j]/D[i]);
1917  }
1918  }
1919
1920  std::vector<double> w(n, 0.0);
1921  // Compute the discarded-fill weights
1922  for (int k=0; k<n; ++k)
1923  {
1924  for (int ii=I[k]; ii<I[k+1]; ++ii)
1925  {
1926  double C_ki = V[ii];
1927  for (int jj=I[k]; jj<I[k+1]; ++jj)
1928  {
1929  if (jj == ii) { continue; }
1930  double C_jk = V[jj];
1931  w[k] += pow(C_jk*C_ki, 2);
1932  }
1933  }
1934  w[k] = sqrt(w[k]);
1935  }
1936
1937  WeightMinHeap w_heap(w);
1938
1939  // Compute ordering
1940  p.SetSize(n);
1941  for (int i=0; i<n; ++i)
1942  {
1943  int pi = w_heap.pop();
1944  p[n-1-i] = pi;
1945  w[pi] = -1;
1946  for (int kk=I[pi]; kk<I[pi+1]; ++kk)
1947  {
1948  int k = J[kk];
1949  if (w_heap.picked(k)) { continue; }
1950  // Recompute weight
1951  w[k] = 0.0;
1952  for (int ii=I[k]; ii<I[k+1]; ++ii)
1953  {
1954  if (w_heap.picked(J[ii])) { continue; }
1955  double C_ki = V[ii];
1956  for (int jj=I[k]; jj<I[k+1]; ++jj)
1957  {
1958  if (jj == ii || w_heap.picked(J[jj])) { continue; }
1959  double C_jk = V[jj];
1960  w[k] += pow(C_jk*C_ki, 2);
1961  }
1962  }
1963  w[k] = sqrt(w[k]);
1964  w_heap.update(k);
1965  }
1966  }
1967 }
1968
1969 BlockILU::BlockILU(int block_size_,
1970  Reordering reordering_,
1971  int k_fill_)
1972  : Solver(0),
1973  block_size(block_size_),
1974  k_fill(k_fill_),
1975  reordering(reordering_)
1976 { }
1977
1979  int block_size_,
1980  Reordering reordering_,
1981  int k_fill_)
1982  : BlockILU(block_size_, reordering_, k_fill_)
1983 {
1984  SetOperator(op);
1985 }
1986
1988 {
1989  const SparseMatrix *A = NULL;
1990 #ifdef MFEM_USE_MPI
1991  const HypreParMatrix *A_par = dynamic_cast<const HypreParMatrix *>(&op);
1992  SparseMatrix A_par_diag;
1993  if (A_par != NULL)
1994  {
1995  A_par->GetDiag(A_par_diag);
1996  A = &A_par_diag;
1997  }
1998 #endif
1999  if (A == NULL)
2000  {
2001  A = dynamic_cast<const SparseMatrix *>(&op);
2002  if (A == NULL)
2003  {
2004  MFEM_ABORT("BlockILU must be created with a SparseMatrix or HypreParMatrix");
2005  }
2006  }
2007  height = op.Height();
2008  width = op.Width();
2009  MFEM_ASSERT(A->Finalized(), "Matrix must be finalized.");
2010  CreateBlockPattern(*A);
2011  Factorize();
2012 }
2013
2014 void BlockILU::CreateBlockPattern(const SparseMatrix &A)
2015 {
2016  MFEM_VERIFY(k_fill == 0, "Only block ILU(0) is currently supported.");
2017  if (A.Height() % block_size != 0)
2018  {
2019  MFEM_ABORT("BlockILU: block size must evenly divide the matrix size");
2020  }
2021
2022  int nrows = A.Height();
2023  const int *I = A.GetI();
2024  const int *J = A.GetJ();
2025  const double *V = A.GetData();
2026  int nnz = 0;
2027  int nblockrows = nrows / block_size;
2028
2029  std::vector<std::set<int>> unique_block_cols(nblockrows);
2030
2031  for (int iblock = 0; iblock < nblockrows; ++iblock)
2032  {
2033  for (int bi = 0; bi < block_size; ++bi)
2034  {
2035  int i = iblock * block_size + bi;
2036  for (int k = I[i]; k < I[i + 1]; ++k)
2037  {
2038  unique_block_cols[iblock].insert(J[k] / block_size);
2039  }
2040  }
2041  nnz += unique_block_cols[iblock].size();
2042  }
2043
2044  if (reordering != Reordering::NONE)
2045  {
2046  SparseMatrix C(nblockrows, nblockrows);
2047  for (int iblock = 0; iblock < nblockrows; ++iblock)
2048  {
2049  for (int jblock : unique_block_cols[iblock])
2050  {
2051  for (int bi = 0; bi < block_size; ++bi)
2052  {
2053  int i = iblock * block_size + bi;
2054  for (int k = I[i]; k < I[i + 1]; ++k)
2055  {
2056  int j = J[k];
2057  if (j >= jblock * block_size && j < (jblock + 1) * block_size)
2058  {
2060  }
2061  }
2062  }
2063  }
2064  }
2065  C.Finalize(false);
2066  double *CV = C.GetData();
2067  for (int i=0; i<C.NumNonZeroElems(); ++i)
2068  {
2069  CV[i] = sqrt(CV[i]);
2070  }
2071
2072  switch (reordering)
2073  {
2076  break;
2077  default:
2078  MFEM_ABORT("BlockILU: unknown reordering")
2079  }
2080  }
2081  else
2082  {
2083  // No reordering: permutation is identity
2084  P.SetSize(nblockrows);
2085  for (int i=0; i<nblockrows; ++i)
2086  {
2087  P[i] = i;
2088  }
2089  }
2090
2091  // Compute inverse permutation
2092  Pinv.SetSize(nblockrows);
2093  for (int i=0; i<nblockrows; ++i)
2094  {
2095  Pinv[P[i]] = i;
2096  }
2097
2098  // Permute columns
2099  std::vector<std::vector<int>> unique_block_cols_perminv(nblockrows);
2100  for (int i=0; i<nblockrows; ++i)
2101  {
2102  std::vector<int> &cols = unique_block_cols_perminv[i];
2103  for (int j : unique_block_cols[P[i]])
2104  {
2105  cols.push_back(Pinv[j]);
2106  }
2107  std::sort(cols.begin(), cols.end());
2108  }
2109
2110  ID.SetSize(nblockrows);
2111  IB.SetSize(nblockrows + 1);
2112  IB[0] = 0;
2113  JB.SetSize(nnz);
2114  AB.SetSize(block_size, block_size, nnz);
2115  DB.SetSize(block_size, block_size, nblockrows);
2116  AB = 0.0;
2117  DB = 0.0;
2118  ipiv.SetSize(block_size*nblockrows);
2119  int counter = 0;
2120
2121  for (int iblock = 0; iblock < nblockrows; ++iblock)
2122  {
2123  int iblock_perm = P[iblock];
2124  for (int jblock : unique_block_cols_perminv[iblock])
2125  {
2126  int jblock_perm = P[jblock];
2127  if (iblock == jblock)
2128  {
2129  ID[iblock] = counter;
2130  }
2131  JB[counter] = jblock;
2132  for (int bi = 0; bi < block_size; ++bi)
2133  {
2134  int i = iblock_perm*block_size + bi;
2135  for (int k = I[i]; k < I[i + 1]; ++k)
2136  {
2137  int j = J[k];
2138  if (j >= jblock_perm*block_size && j < (jblock_perm + 1)*block_size)
2139  {
2140  int bj = j - jblock_perm*block_size;
2141  double val = V[k];
2142  AB(bi, bj, counter) = val;
2143  // Extract the diagonal
2144  if (iblock == jblock)
2145  {
2146  DB(bi, bj, iblock) = val;
2147  }
2148  }
2149  }
2150  }
2151  ++counter;
2152  }
2153  IB[iblock + 1] = counter;
2154  }
2155 }
2156
2157 void BlockILU::Factorize()
2158 {
2159  int nblockrows = Height()/block_size;
2160
2161  // Precompute LU factorization of diagonal blocks
2162  for (int i=0; i<nblockrows; ++i)
2163  {
2164  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2165  factorization.Factor(block_size);
2166  }
2167
2168  // Note: we use UseExternalData to extract submatrices from the tensor AB
2169  // instead of the DenseTensor call operator, because the call operator does
2170  // not allow for two simultaneous submatrix views into the same tensor
2171  DenseMatrix A_ik, A_ij, A_kj;
2172  // Loop over block rows (starting with second block row)
2173  for (int i=1; i<nblockrows; ++i)
2174  {
2175  // Find all nonzeros to the left of the diagonal in row i
2176  for (int kk=IB[i]; kk<IB[i+1]; ++kk)
2177  {
2178  int k = JB[kk];
2179  // Make sure we're still to the left of the diagonal
2180  if (k == i) { break; }
2181  if (k > i)
2182  {
2183  MFEM_ABORT("Matrix must be sorted with nonzero diagonal");
2184  }
2185  LUFactors A_kk_inv(DB.GetData(k), &ipiv[k*block_size]);
2186  A_ik.UseExternalData(&AB(0,0,kk), block_size, block_size);
2187  // A_ik = A_ik * A_kk^{-1}
2188  A_kk_inv.RightSolve(block_size, block_size, A_ik.GetData());
2189  // Modify everything to the right of k in row i
2190  for (int jj=kk+1; jj<IB[i+1]; ++jj)
2191  {
2192  int j = JB[jj];
2193  if (j <= k) { continue; } // Superfluous because JB is sorted?
2194  A_ij.UseExternalData(&AB(0,0,jj), block_size, block_size);
2195  for (int ll=IB[k]; ll<IB[k+1]; ++ll)
2196  {
2197  int l = JB[ll];
2198  if (l == j)
2199  {
2200  A_kj.UseExternalData(&AB(0,0,ll), block_size, block_size);
2201  // A_ij = A_ij - A_ik*A_kj;
2203  // If we need to, update diagonal factorization
2204  if (j == i)
2205  {
2206  DB(i) = A_ij;
2207  LUFactors factorization(DB.GetData(i), &ipiv[i*block_size]);
2208  factorization.Factor(block_size);
2209  }
2210  break;
2211  }
2212  }
2213  }
2214  }
2215  }
2216 }
2217
2218 void BlockILU::Mult(const Vector &b, Vector &x) const
2219 {
2220  MFEM_ASSERT(height > 0, "BlockILU(0) preconditioner is not constructed");
2221  int nblockrows = Height()/block_size;
2222  y.SetSize(Height());
2223
2224  DenseMatrix B;
2225  Vector yi, yj, xi, xj;
2226  Vector tmp(block_size);
2227  // Forward substitute to solve Ly = b
2228  // Implicitly, L has identity on the diagonal
2229  y = 0.0;
2230  for (int i=0; i<nblockrows; ++i)
2231  {
2232  yi.SetDataAndSize(&y[i*block_size], block_size);
2233  for (int ib=0; ib<block_size; ++ib)
2234  {
2235  yi[ib] = b[ib + P[i]*block_size];
2236  }
2237  for (int k=IB[i]; k<ID[i]; ++k)
2238  {
2239  int j = JB[k];
2240  const DenseMatrix &L_ij = AB(k);
2241  yj.SetDataAndSize(&y[j*block_size], block_size);
2242  // y_i = y_i - L_ij*y_j
2244  }
2245  }
2246  // Backward substitution to solve Ux = y
2247  for (int i=nblockrows-1; i >= 0; --i)
2248  {
2249  xi.SetDataAndSize(&x[P[i]*block_size], block_size);
2250  for (int ib=0; ib<block_size; ++ib)
2251  {
2252  xi[ib] = y[ib + i*block_size];
2253  }
2254  for (int k=ID[i]+1; k<IB[i+1]; ++k)
2255  {
2256  int j = JB[k];
2257  const DenseMatrix &U_ij = AB(k);
2258  xj.SetDataAndSize(&x[P[j]*block_size], block_size);
2259  // x_i = x_i - U_ij*x_j
2261  }
2262  LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
2263  // x_i = D_ii^{-1} x_i
2264  A_ii_inv.Solve(block_size, 1, xi);
2265  }
2266 }
2267
2268 #ifdef MFEM_USE_SUITESPARSE
2269
2271 {
2272  mat = NULL;
2273  Numeric = NULL;
2274  AI = AJ = NULL;
2275  if (!use_long_ints)
2276  {
2277  umfpack_di_defaults(Control);
2278  }
2279  else
2280  {
2281  umfpack_dl_defaults(Control);
2282  }
2283 }
2284
2286 {
2287  int *Ap, *Ai;
2288  void *Symbolic;
2289  double *Ax;
2290
2291  if (Numeric)
2292  {
2293  if (!use_long_ints)
2294  {
2295  umfpack_di_free_numeric(&Numeric);
2296  }
2297  else
2298  {
2299  umfpack_dl_free_numeric(&Numeric);
2300  }
2301  }
2302
2303  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
2304  MFEM_VERIFY(mat, "not a SparseMatrix");
2305
2306  // UMFPack requires that the column-indices in mat corresponding to each
2307  // row be sorted.
2308  // Generally, this will modify the ordering of the entries of mat.
2310
2311  height = mat->Height();
2312  width = mat->Width();
2313  MFEM_VERIFY(width == height, "not a square matrix");
2314
2315  Ap = mat->GetI();
2316  Ai = mat->GetJ();
2317  Ax = mat->GetData();
2318
2319  if (!use_long_ints)
2320  {
2321  int status = umfpack_di_symbolic(width, width, Ap, Ai, Ax, &Symbolic,
2322  Control, Info);
2323  if (status < 0)
2324  {
2325  umfpack_di_report_info(Control, Info);
2326  umfpack_di_report_status(Control, status);
2327  mfem_error("UMFPackSolver::SetOperator :"
2328  " umfpack_di_symbolic() failed!");
2329  }
2330
2331  status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric,
2332  Control, Info);
2333  if (status < 0)
2334  {
2335  umfpack_di_report_info(Control, Info);
2336  umfpack_di_report_status(Control, status);
2337  mfem_error("UMFPackSolver::SetOperator :"
2338  " umfpack_di_numeric() failed!");
2339  }
2340  umfpack_di_free_symbolic(&Symbolic);
2341  }
2342  else
2343  {
2344  SuiteSparse_long status;
2345
2346  delete [] AJ;
2347  delete [] AI;
2348  AI = new SuiteSparse_long[width + 1];
2349  AJ = new SuiteSparse_long[Ap[width]];
2350  for (int i = 0; i <= width; i++)
2351  {
2352  AI[i] = (SuiteSparse_long)(Ap[i]);
2353  }
2354  for (int i = 0; i < Ap[width]; i++)
2355  {
2356  AJ[i] = (SuiteSparse_long)(Ai[i]);
2357  }
2358
2359  status = umfpack_dl_symbolic(width, width, AI, AJ, Ax, &Symbolic,
2360  Control, Info);
2361  if (status < 0)
2362  {
2363  umfpack_dl_report_info(Control, Info);
2364  umfpack_dl_report_status(Control, status);
2365  mfem_error("UMFPackSolver::SetOperator :"
2366  " umfpack_dl_symbolic() failed!");
2367  }
2368
2369  status = umfpack_dl_numeric(AI, AJ, Ax, Symbolic, &Numeric,
2370  Control, Info);
2371  if (status < 0)
2372  {
2373  umfpack_dl_report_info(Control, Info);
2374  umfpack_dl_report_status(Control, status);
2375  mfem_error("UMFPackSolver::SetOperator :"
2376  " umfpack_dl_numeric() failed!");
2377  }
2378  umfpack_dl_free_symbolic(&Symbolic);
2379  }
2380 }
2381
2382 void UMFPackSolver::Mult(const Vector &b, Vector &x) const
2383 {
2384  if (mat == NULL)
2385  mfem_error("UMFPackSolver::Mult : matrix is not set!"
2386  " Call SetOperator first!");
2387
2388  if (!use_long_ints)
2389  {
2390  int status =
2391  umfpack_di_solve(UMFPACK_At, mat->GetI(), mat->GetJ(),
2392  mat->GetData(), x, b, Numeric, Control, Info);
2393  umfpack_di_report_info(Control, Info);
2394  if (status < 0)
2395  {
2396  umfpack_di_report_status(Control, status);
2397  mfem_error("UMFPackSolver::Mult : umfpack_di_solve() failed!");
2398  }
2399  }
2400  else
2401  {
2402  SuiteSparse_long status =
2403  umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->GetData(), x, b,
2404  Numeric, Control, Info);
2405  umfpack_dl_report_info(Control, Info);
2406  if (status < 0)
2407  {
2408  umfpack_dl_report_status(Control, status);
2409  mfem_error("UMFPackSolver::Mult : umfpack_dl_solve() failed!");
2410  }
2411  }
2412 }
2413
2415 {
2416  if (mat == NULL)
2417  mfem_error("UMFPackSolver::MultTranspose : matrix is not set!"
2418  " Call SetOperator first!");
2419
2420  if (!use_long_ints)
2421  {
2422  int status =
2423  umfpack_di_solve(UMFPACK_A, mat->GetI(), mat->GetJ(),
2424  mat->GetData(), x, b, Numeric, Control, Info);
2425  umfpack_di_report_info(Control, Info);
2426  if (status < 0)
2427  {
2428  umfpack_di_report_status(Control, status);
2429  mfem_error("UMFPackSolver::MultTranspose :"
2430  " umfpack_di_solve() failed!");
2431  }
2432  }
2433  else
2434  {
2435  SuiteSparse_long status =
2436  umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->GetData(), x, b,
2437  Numeric, Control, Info);
2438  umfpack_dl_report_info(Control, Info);
2439  if (status < 0)
2440  {
2441  umfpack_dl_report_status(Control, status);
2442  mfem_error("UMFPackSolver::MultTranspose :"
2443  " umfpack_dl_solve() failed!");
2444  }
2445  }
2446 }
2447
2449 {
2450  delete [] AJ;
2451  delete [] AI;
2452  if (Numeric)
2453  {
2454  if (!use_long_ints)
2455  {
2456  umfpack_di_free_numeric(&Numeric);
2457  }
2458  else
2459  {
2460  umfpack_dl_free_numeric(&Numeric);
2461  }
2462  }
2463 }
2464
2466 {
2467  klu_defaults(&Common);
2468 }
2469
2471 {
2472  if (Numeric)
2473  {
2474  MFEM_ASSERT(Symbolic != 0,
2475  "Had Numeric pointer in KLU, but not Symbolic");
2476  klu_free_symbolic(&Symbolic, &Common);
2477  Symbolic = 0;
2478  klu_free_numeric(&Numeric, &Common);
2479  Numeric = 0;
2480  }
2481
2482  mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
2483  MFEM_VERIFY(mat != NULL, "not a SparseMatrix");
2484
2485  // KLU requires that the column-indices in mat corresponding to each row be
2486  // sorted. Generally, this will modify the ordering of the entries of mat.
2488
2489  height = mat->Height();
2490  width = mat->Width();
2491  MFEM_VERIFY(width == height, "not a square matrix");
2492
2493  int * Ap = mat->GetI();
2494  int * Ai = mat->GetJ();
2495  double * Ax = mat->GetData();
2496
2497  Symbolic = klu_analyze( height, Ap, Ai, &Common);
2498  Numeric = klu_factor(Ap, Ai, Ax, Symbolic, &Common);
2499 }
2500
2501 void KLUSolver::Mult(const Vector &b, Vector &x) const
2502 {
2503  MFEM_VERIFY(mat != NULL,
2504  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
2505
2506  int n = mat->Height();
2507  int numRhs = 1;
2508  // Copy B into X, so we can pass it in and overwrite it.
2509  x = b;
2510  // Solve the transpose, since KLU thinks the matrix is compressed column
2511  // format.
2512  klu_tsolve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
2513 }
2514
2515 void KLUSolver::MultTranspose(const Vector &b, Vector &x) const
2516 {
2517  MFEM_VERIFY(mat != NULL,
2518  "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
2519
2520  int n = mat->Height();
2521  int numRhs = 1;
2522  // Copy B into X, so we can pass it in and overwrite it.
2523  x = b;
2524  // Solve the regular matrix, not the transpose, since KLU thinks the matrix
2525  // is compressed column format.
2526  klu_solve( Symbolic, Numeric, n, numRhs, x.GetData(), &Common);
2527 }
2528
2530 {
2531  klu_free_symbolic (&Symbolic, &Common) ;
2532  klu_free_numeric (&Numeric, &Common) ;
2533  Symbolic = 0;
2534  Numeric = 0;
2535 }
2536
2537 #endif // MFEM_USE_SUITESPARSE
2538
2539 }
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:566
int Size() const
Logical size of the array.
Definition: array.hpp:124
Definition: solvers.hpp:152
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:2414
double Info[UMFPACK_INFO]
Definition: solvers.hpp:592
void Mult(const Vector &b, Vector &x) const
Solve the system LUx = b, where L and U are the block ILU factors.
Definition: solvers.cpp:2218
SparseMatrix * mat
Definition: solvers.hpp:584
int GetNumIterations() const
Definition: solvers.hpp:68
const Operator * oper
Definition: solvers.hpp:43
SuiteSparse_long * AI
Definition: solvers.hpp:586
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:358
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
virtual void ProcessNewState(const Vector &x) const
This method can be overloaded in derived classes to perform computations that need knowledge of the n...
Definition: solvers.hpp:330
void MinimumDiscardedFillOrdering(SparseMatrix &C, Array< int > &p)
Definition: solvers.cpp:1903
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:89
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:711
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:146
void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:86
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:151
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:774
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
int * GetI()
Return the array I.
Definition: sparsemat.hpp:141
void SetSize(int i, int j, int k)
Definition: densemat.hpp:768
void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag.
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:2515
void SetOperator(const Operator &op)
Definition: solvers.cpp:1987
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:504
MINRES method.
Definition: solvers.hpp:262
const Operator * GetC() const
Definition: solvers.hpp:382
double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:341
Reordering
The reordering method used by the BlockILU factorization.
Definition: solvers.hpp:497
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:166
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:368
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:238
double Norm(const Vector &x) const
Definition: solvers.hpp:54
const Vector * d_hi
Definition: solvers.hpp:364
double GetFinalNorm() const
Definition: solvers.hpp:70
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:151
int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, double &tol, double atol, int printit)
BiCGSTAB method. (tolerances are squared)
Definition: solvers.cpp:1107
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:180
bool IsFinite(const double &val)
Definition: vector.hpp:372
const Vector * x_lo
Definition: solvers.hpp:364
const Vector * x_hi
Definition: solvers.hpp:364
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1926
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:594
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1144
void SetLinearConstraint(const Vector &_w, double _a)
Definition: solvers.cpp:1632
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:275
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
Definition: vector.hpp:349
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:97
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:433
int GetNumConstraints() const
Definition: solvers.cpp:1605
Data type sparse matrix.
Definition: sparsemat.hpp:40
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
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:1302
virtual ~UMFPackSolver()
Definition: solvers.cpp:2448
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:542
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
double b
Definition: lissajous.cpp:42
void SetMaxIter(int max_it)
Definition: solvers.hpp:65
double delta
Definition: lissajous.cpp:43
SuiteSparse_long * AJ
Definition: solvers.hpp:586
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:197
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:516
const T * Read(bool on_dev=true) const
Definition: array.hpp:273
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:1638
double * GetData(int k)
Definition: densemat.hpp:813
void Setup(const Vector &diag)
Definition: solvers.cpp:140
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:529
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:939
void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3019
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:2470
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1130
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:54
int GetConverged() const
Definition: solvers.hpp:69
klu_symbolic * Symbolic
Definition: solvers.hpp:625
BlockILU(int block_size_, Reordering reordering_=Reordering::MINIMUM_DISCARDED_FILL, int k_fill_=0)
Definition: solvers.cpp:1969
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1339
OperatorJacobiSmoother(const BilinearForm &a, const Array< int > &ess_tdof_list, const double damping=1.0)
Definition: solvers.cpp:109
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:428
Stationary linear iteration: x &lt;- x + B (b - A x)
Definition: solvers.hpp:120
void SetInequalityConstraint(const Vector &dl, const Vector &dh)
Definition: solvers.cpp:1587
void SetEqualityConstraint(const Vector &c)
Definition: solvers.cpp:1579
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:592
void SetAbsTol(double atol)
Definition: solvers.hpp:64
void SetRelTol(double rtol)
Definition: solvers.hpp:63
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:591
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:978
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:2501
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:385
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
GMRES method.
Definition: solvers.hpp:184
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
Definition: solvers.cpp:573
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:125
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:205
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.cpp:1613
void SetBounds(const Vector &_lo, const Vector &_hi)
Definition: solvers.cpp:1626
double a
Definition: lissajous.cpp:41
void UpdateVectors()
Definition: solvers.cpp:174
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:246
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:190
klu_common Common
Definition: solvers.hpp:646
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:489
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:325
OptimizationProblem(int insize, const Operator *C_, const Operator *D_)
In parallel, insize is the number of the local true dofs.
Definition: solvers.cpp:1569
void SetSolutionBounds(const Vector &xl, const Vector &xh)
Definition: solvers.cpp:1597
const double * Read(bool on_dev=true) const
Definition: vector.hpp:333
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: solvers.cpp:544
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:245
virtual void Mult(const Vector &xt, Vector &x) const
Definition: solvers.cpp:1645
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:166
const double alpha
Definition: ex15.cpp:336
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:260
Vector data type.
Definition: vector.hpp:48
const Vector * d_lo
Definition: solvers.hpp:364
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
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:1418
virtual ~KLUSolver()
Definition: solvers.cpp:2529
void UpdateVectors()
Definition: solvers.cpp:351
BiCGSTAB method.
Definition: solvers.hpp:231
SparseMatrix * mat
Definition: solvers.hpp:624
Base class for solvers.
Definition: operator.hpp:500
const OptimizationProblem * problem
Definition: solvers.hpp:397
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:187
void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
Definition: densemat.cpp:2024
const Operator * D
Definition: solvers.hpp:363
const Vector * c_e
Definition: solvers.hpp:364
klu_numeric * Numeric
Definition: solvers.hpp:626
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:2382
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:134
const Operator * C
Not owned, some can remain unused (NULL).
Definition: solvers.hpp:363
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:323
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:942
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1328
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:2285