GetFEM  5.4.2
getfem_model_solvers.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2004-2020 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**
33  @file getfem_model_solvers.h
34  @author Yves Renard <Yves.Renard@insa-lyon.fr>
35  @date June 15, 2004.
36  @brief Standard solvers for model bricks
37  @see getfem_modeling.h
38 */
39 
40 #ifndef GETFEM_MODEL_SOLVERS_H__
41 #define GETFEM_MODEL_SOLVERS_H__
42 #include "getfem_models.h"
43 #include "gmm/gmm_MUMPS_interface.h"
44 #include "gmm/gmm_iter.h"
45 #include "gmm/gmm_iter_solvers.h"
46 #include "gmm/gmm_dense_qr.h"
47 
48 //#include "gmm/gmm_inoutput.h"
49 
50 
51 namespace getfem {
52 
53  template <typename T> struct sort_abs_val_
54  { bool operator()(T x, T y) { return (gmm::abs(x) < gmm::abs(y)); } };
55 
56  template <typename MAT> void print_eigval(const MAT &M) {
57  // just to test a stiffness matrix if needing
58  typedef typename gmm::linalg_traits<MAT>::value_type T;
59  size_type n = gmm::mat_nrows(M);
60  gmm::dense_matrix<T> MM(n, n), Q(n, n);
61  std::vector<T> eigval(n);
62  gmm::copy(M, MM);
63  // gmm::symmetric_qr_algorithm(MM, eigval, Q);
64  gmm::implicit_qr_algorithm(MM, eigval, Q);
65  std::sort(eigval.begin(), eigval.end(), sort_abs_val_<T>());
66  cout << "eival = " << eigval << endl;
67 // cout << "vectp : " << gmm::mat_col(Q, n-1) << endl;
68 // cout << "vectp : " << gmm::mat_col(Q, n-2) << endl;
69 // double emax, emin;
70 // cout << "condition number" << condition_number(MM,emax,emin) << endl;
71 // cout << "emin = " << emin << endl;
72 // cout << "emax = " << emax << endl;
73  }
74 
75 
76  /* ***************************************************************** */
77  /* Linear solvers definition */
78  /* ***************************************************************** */
79 
80  template <typename MAT, typename VECT>
81  struct abstract_linear_solver {
82  typedef MAT MATRIX;
83  typedef VECT VECTOR;
84  virtual void operator ()(const MAT &, VECT &, const VECT &,
85  gmm::iteration &) const = 0;
86  virtual ~abstract_linear_solver() {}
87  };
88 
89  template <typename MAT, typename VECT>
90  struct linear_solver_cg_preconditioned_ildlt
91  : public abstract_linear_solver<MAT, VECT> {
92  void operator ()(const MAT &M, VECT &x, const VECT &b,
93  gmm::iteration &iter) const {
94  gmm::ildlt_precond<MAT> P(M);
95  gmm::cg(M, x, b, P, iter);
96  if (!iter.converged()) GMM_WARNING2("cg did not converge!");
97  }
98  };
99 
100  template <typename MAT, typename VECT>
101  struct linear_solver_gmres_preconditioned_ilu
102  : public abstract_linear_solver<MAT, VECT> {
103  void operator ()(const MAT &M, VECT &x, const VECT &b,
104  gmm::iteration &iter) const {
105  gmm::ilu_precond<MAT> P(M);
106  gmm::gmres(M, x, b, P, 500, iter);
107  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
108  }
109  };
110 
111  template <typename MAT, typename VECT>
112  struct linear_solver_gmres_unpreconditioned
113  : public abstract_linear_solver<MAT, VECT> {
114  void operator ()(const MAT &M, VECT &x, const VECT &b,
115  gmm::iteration &iter) const {
116  gmm::identity_matrix P;
117  gmm::gmres(M, x, b, P, 500, iter);
118  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
119  }
120  };
121 
122  template <typename MAT, typename VECT>
123  struct linear_solver_gmres_preconditioned_ilut
124  : public abstract_linear_solver<MAT, VECT> {
125  void operator ()(const MAT &M, VECT &x, const VECT &b,
126  gmm::iteration &iter) const {
127  gmm::ilut_precond<MAT> P(M, 40, 1E-7);
128  gmm::gmres(M, x, b, P, 500, iter);
129  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
130  }
131  };
132 
133  template <typename MAT, typename VECT>
134  struct linear_solver_gmres_preconditioned_ilutp
135  : public abstract_linear_solver<MAT, VECT> {
136  void operator ()(const MAT &M, VECT &x, const VECT &b,
137  gmm::iteration &iter) const {
138  gmm::ilutp_precond<MAT> P(M, 20, 1E-7);
139  gmm::gmres(M, x, b, P, 500, iter);
140  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
141  }
142  };
143 
144  template <typename MAT, typename VECT>
145  struct linear_solver_superlu
146  : public abstract_linear_solver<MAT, VECT> {
147  void operator ()(const MAT &M, VECT &x, const VECT &b,
148  gmm::iteration &iter) const {
149  double rcond;
150  /*gmm::HarwellBoeing_IO::write("test.hb", M);
151  std::fstream f("bbb", std::ios::out);
152  for (unsigned i=0; i < gmm::vect_size(b); ++i) f << b[i] << "\n";*/
153  int info = SuperLU_solve(M, x, b, rcond);
154  iter.enforce_converged(info == 0);
155  if (iter.get_noisy()) cout << "condition number: " << 1.0/rcond<< endl;
156  }
157  };
158 
159  template <typename MAT, typename VECT>
160  struct linear_solver_dense_lu : public abstract_linear_solver<MAT, VECT> {
161  void operator ()(const MAT &M, VECT &x, const VECT &b,
162  gmm::iteration &iter) const {
163  typedef typename gmm::linalg_traits<MAT>::value_type T;
164  gmm::dense_matrix<T> MM(gmm::mat_nrows(M),gmm::mat_ncols(M));
165  gmm::copy(M, MM);
166  gmm::lu_solve(MM, x, b);
167  iter.enforce_converged(true);
168  }
169  };
170 
171 #ifdef GMM_USES_MUMPS
172  template <typename MAT, typename VECT>
173  struct linear_solver_mumps : public abstract_linear_solver<MAT, VECT> {
174  void operator ()(const MAT &M, VECT &x, const VECT &b,
175  gmm::iteration &iter) const {
176  bool ok = gmm::MUMPS_solve(M, x, b, false);
177  iter.enforce_converged(ok);
178  }
179  };
180  template <typename MAT, typename VECT>
181  struct linear_solver_mumps_sym : public abstract_linear_solver<MAT, VECT> {
182  void operator ()(const MAT &M, VECT &x, const VECT &b,
183  gmm::iteration &iter) const {
184  bool ok = gmm::MUMPS_solve(M, x, b, true);
185  iter.enforce_converged(ok);
186  }
187  };
188 #endif
189 
190 #if GETFEM_PARA_LEVEL > 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
191  template <typename MAT, typename VECT>
192  struct linear_solver_distributed_mumps
193  : public abstract_linear_solver<MAT, VECT> {
194  void operator ()(const MAT &M, VECT &x, const VECT &b,
195  gmm::iteration &iter) const {
196  double tt_ref=MPI_Wtime();
197  bool ok = MUMPS_distributed_matrix_solve(M, x, b, false);
198  iter.enforce_converged(ok);
199  if (MPI_IS_MASTER()) cout<<"UNSYMMETRIC MUMPS time "<< MPI_Wtime() - tt_ref<<endl;
200  }
201  };
202 
203  template <typename MAT, typename VECT>
204  struct linear_solver_distributed_mumps_sym
205  : public abstract_linear_solver<MAT, VECT> {
206  void operator ()(const MAT &M, VECT &x, const VECT &b,
207  gmm::iteration &iter) const {
208  double tt_ref=MPI_Wtime();
209  bool ok = MUMPS_distributed_matrix_solve(M, x, b, true);
210  iter.enforce_converged(ok);
211  if (MPI_IS_MASTER()) cout<<"SYMMETRIC MUMPS time "<< MPI_Wtime() - tt_ref<<endl;
212  }
213  };
214 #endif
215 
216 
217  /* ***************************************************************** */
218  /* Newton Line search definition */
219  /* ***************************************************************** */
220 
221  struct abstract_newton_line_search {
222  double conv_alpha, conv_r;
223  size_t it, itmax, glob_it;
224  // size_t tot_it;
225  virtual void init_search(double r, size_t git, double R0 = 0.0) = 0;
226  virtual double next_try(void) = 0;
227  virtual bool is_converged(double, double R1 = 0.0) = 0;
228  virtual double converged_value(void) {
229  // tot_it += it; cout << "tot_it = " << tot_it << endl; it = 0;
230  return conv_alpha;
231  };
232  virtual double converged_residual(void) { return conv_r; };
233  virtual ~abstract_newton_line_search() { }
234  };
235 
236  // Dummy linesearch for the newton with step control
237  struct newton_search_with_step_control : public abstract_newton_line_search {
238 
239  virtual void init_search(double /*r*/, size_t /*git*/, double /*R0*/ = 0.0)
240  { GMM_ASSERT1(false, "Not to be used"); }
241 
242  virtual double next_try(void)
243  { GMM_ASSERT1(false, "Not to be used"); }
244 
245  virtual bool is_converged(double /*r*/, double /*R1*/ = 0.0)
246  { GMM_ASSERT1(false, "Not to be used"); }
247 
248  newton_search_with_step_control() {}
249  };
250 
251 
252  struct quadratic_newton_line_search : public abstract_newton_line_search {
253  double R0_, R1_;
254 
255  double alpha, alpha_mult, first_res, alpha_max_ratio, alpha_min;
256  virtual void init_search(double r, size_t git, double R0 = 0.0) {
257  GMM_ASSERT1(R0 != 0.0, "You have to specify R0");
258  glob_it = git;
259  conv_alpha = alpha = double(1); conv_r = first_res = r; it = 0;
260  R0_ = R0;
261  }
262  virtual double next_try(void) {
263  ++it;
264  if (it == 1) return double(1);
265  GMM_ASSERT1(R1_ != 0.0, "You have to specify R1");
266  double a = R0_ / R1_;
267  return (a < 0) ? (a*0.5 + sqrt(a*a*0.25-a)) : a*0.5;
268  }
269  virtual bool is_converged(double r, double R1 = 0.0) {
270  conv_r = r;
271  R1_ = R1; return ((gmm::abs(R1_) < gmm::abs(R0_*0.5)) || it >= itmax);
272  }
273  quadratic_newton_line_search(size_t imax = size_t(-1)) { itmax = imax; }
274  };
275 
276 
277  struct simplest_newton_line_search : public abstract_newton_line_search {
278  double alpha, alpha_mult, first_res, alpha_max_ratio, alpha_min,
279  alpha_threshold_res;
280  virtual void init_search(double r, size_t git, double = 0.0) {
281  glob_it = git;
282  conv_alpha = alpha = double(1); conv_r = first_res = r; it = 0;
283  }
284  virtual double next_try(void)
285  { conv_alpha = alpha; alpha *= alpha_mult; ++it; return conv_alpha; }
286  virtual bool is_converged(double r, double = 0.0) {
287  conv_r = r;
288  return ((it <= 1 && r < first_res)
289  || (r <= first_res * alpha_max_ratio && r <= alpha_threshold_res)
290  || (conv_alpha <= alpha_min && r < first_res * 1e5)
291  || it >= itmax);
292  }
293  simplest_newton_line_search
294  (size_t imax = size_t(-1), double a_max_ratio = 6.0/5.0,
295  double a_min = 1.0/1000.0, double a_mult = 3.0/5.0,
296  double a_threshold_res = 1e50)
297  : alpha_mult(a_mult), alpha_max_ratio(a_max_ratio), alpha_min(a_min),
298  alpha_threshold_res(a_threshold_res)
299  { itmax = imax; }
300  };
301 
302  struct default_newton_line_search : public abstract_newton_line_search {
303  // This line search try to detect where is the minimum, dividing the step
304  // by a factor two each time.
305  // - it stops if the residual is less than the previous residual
306  // times alpha_min_ratio (= 0.9).
307  // - if the minimal step is reached with a residual greater than
308  // the previous residual times alpha_min_ratio then it decides
309  // between two options :
310  // - return with the step corresponding to the smallest residual
311  // - return with a greater residual corresponding to a residual
312  // less than the previous residual times alpha_max_ratio.
313  // the decision is taken regarding the previous iterations.
314  // - in order to shorten the line search, the process stops when
315  // the residual increases three times consecutively.
316  // possible improvment : detect the curvature at the origin
317  // (only one evaluation) and take it into account.
318  // Fitted to some experiments in contrib/tests_newton
319 
320  double alpha, alpha_old, alpha_mult, first_res, alpha_max_ratio;
321  double alpha_min_ratio, alpha_min;
322  size_type count, count_pat;
323  bool max_ratio_reached;
324  double alpha_max_ratio_reached, r_max_ratio_reached;
325  size_type it_max_ratio_reached;
326 
327  virtual void init_search(double r, size_t git, double = 0.0);
328  virtual double next_try(void);
329  virtual bool is_converged(double r, double = 0.0);
330  default_newton_line_search(void) { count_pat = 0; }
331  };
332 
333  /* the former default_newton_line_search */
334  struct basic_newton_line_search : public abstract_newton_line_search {
335  double alpha, alpha_mult, first_res, alpha_max_ratio;
336  double alpha_min, prev_res, alpha_max_augment;
337  virtual void init_search(double r, size_t git, double = 0.0) {
338  glob_it = git;
339  conv_alpha = alpha = double(1);
340  prev_res = conv_r = first_res = r; it = 0;
341  }
342  virtual double next_try(void)
343  { conv_alpha = alpha; alpha *= alpha_mult; ++it; return conv_alpha; }
344  virtual bool is_converged(double r, double = 0.0) {
345  if (glob_it == 0 || (r < first_res / double(2))
346  || (conv_alpha <= alpha_min && r < first_res * alpha_max_augment)
347  || it >= itmax)
348  { conv_r = r; return true; }
349  if (it > 1 && r > prev_res && prev_res < alpha_max_ratio * first_res)
350  return true;
351  prev_res = conv_r = r;
352  return false;
353  }
354  basic_newton_line_search
355  (size_t imax = size_t(-1),
356  double a_max_ratio = 5.0/3.0,
357  double a_min = 1.0/1000.0, double a_mult = 3.0/5.0, double a_augm = 2.0)
358  : alpha_mult(a_mult), alpha_max_ratio(a_max_ratio),
359  alpha_min(a_min), alpha_max_augment(a_augm) { itmax = imax; }
360  };
361 
362 
363  struct systematic_newton_line_search : public abstract_newton_line_search {
364  double alpha, alpha_mult, first_res;
365  double alpha_min, prev_res;
366  bool first;
367  virtual void init_search(double r, size_t git, double = 0.0) {
368  glob_it = git;
369  conv_alpha = alpha = double(1);
370  prev_res = conv_r = first_res = r; it = 0; first = true;
371  }
372  virtual double next_try(void)
373  { double a = alpha; alpha *= alpha_mult; ++it; return a; }
374  virtual bool is_converged(double r, double = 0.0) {
375  // cout << "a = " << alpha / alpha_mult << " r = " << r << endl;
376  if (r < conv_r || first)
377  { conv_r = r; conv_alpha = alpha / alpha_mult; first = false; }
378  if ((alpha <= alpha_min*alpha_mult) || it >= itmax) return true;
379  return false;
380  }
381  systematic_newton_line_search
382  (size_t imax = size_t(-1),
383  double a_min = 1.0/10000.0, double a_mult = 3.0/5.0)
384  : alpha_mult(a_mult), alpha_min(a_min) { itmax = imax; }
385  };
386 
387 
388  /* ***************************************************************** */
389  /* Newton(-Raphson) algorithm with step control. */
390  /* ***************************************************************** */
391  /* Still solves a problem F(X) = 0 starting at X0 but setting */
392  /* B0 = F(X0) and solving */
393  /* F(X) = (1-alpha)B0 (1) */
394  /* with alpha growing from 0 to 1. */
395  /* A very basic line search is applied. */
396  /* */
397  /* Step 0 : set alpha0 = 0, alpha = 1, X0 given and B0 = F(X0). */
398  /* Step 1 : Set Ri = F(Xi) - (1-alpha)B0 */
399  /* If ||Ri|| < rho, Xi+1 = Xi and go to step 2 */
400  /* Perform Newton step on problem (1) */
401  /* If the Newton do not converge (stagnation) */
402  /* alpha <- max(alpha0+1E-4, (alpha+alpha0)/2) */
403  /* Loop on step 1 with Xi unchanged */
404  /* Step 2 : if alpha=1 stop */
405  /* else alpha0 <- alpha, */
406  /* alpha <- min(1,alpha0+2*(alpha-alpha0)), */
407  /* Go to step 1 with Xi+1 */
408  /* being the result of Newton iterations of step1. */
409  /* */
410  /*********************************************************************/
411 
412  template <typename PB>
413  void Newton_with_step_control(PB &pb, gmm::iteration &iter)
414  {
415  typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T;
416  typedef typename gmm::number_traits<T>::magnitude_type R;
417  gmm::iteration iter_linsolv0 = iter;
418  iter_linsolv0.reduce_noisy();
419  iter_linsolv0.set_resmax(iter.get_resmax()/20.0);
420  iter_linsolv0.set_maxiter(10000); // arbitrary
421 
422  pb.compute_residual();
423  R approx_eln = pb.approx_external_load_norm();
424  if (approx_eln == R(0)) approx_eln = R(1);
425 
426  typename PB::VECTOR b0(gmm::vect_size(pb.residual()));
427  gmm::copy(pb.residual(), b0);
428  typename PB::VECTOR Xi(gmm::vect_size(pb.residual()));
429  gmm::copy(pb.state_vector(), Xi);
430 
431  typename PB::VECTOR dr(gmm::vect_size(pb.residual()));
432 
433  R alpha0(0), alpha(1), res0(gmm::vect_norm1(b0)), minres(res0);
434  // const newton_search_with_step_control *ls
435  // = dynamic_cast<const newton_search_with_step_control *>(&(pb.ls));
436  // GMM_ASSERT1(ls, "Internal error");
437  size_type nit = 0, stagn = 0;
438  R coeff = R(2);
439 
440  scalar_type crit = pb.residual_norm() / approx_eln;
441  if (iter.finished(crit)) return;
442  for(;;) {
443 
444  crit = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha))
445  / approx_eln;
446  if (!iter.converged(crit)) {
447  gmm::iteration iter_linsolv = iter_linsolv0;
448 
449  int is_singular = 1;
450  while (is_singular) { // Linear system solve
451  gmm::clear(dr);
452  pb.add_to_residual(b0, alpha-R(1)); // canceled at next compute_residual
453  iter_linsolv.init();
454  if (iter.get_noisy() > 1)
455  cout << "starting tangent matrix computation" << endl;
456  pb.compute_tangent_matrix();
457  if (iter.get_noisy() > 1)
458  cout << "starting linear solver" << endl;
459  pb.linear_solve(dr, iter_linsolv);
460  if (!iter_linsolv.converged()) {
461  is_singular++;
462  if (is_singular <= 4) {
463  if (iter.get_noisy())
464  cout << "Singular tangent matrix:"
465  " perturbation of the state vector." << endl;
466  pb.perturbation();
467  pb.compute_residual(); // cancels add_to_residual
468  } else {
469  if (iter.get_noisy())
470  cout << "Singular tangent matrix: perturbation failed, "
471  << "aborting." << endl;
472  return;
473  }
474  }
475  else is_singular = 0;
476  }
477  if (iter.get_noisy() > 1) cout << "linear solver done" << endl;
478 
479  gmm::add(dr, pb.state_vector());
480  pb.compute_residual(); // cancels add_to_residual
481  R res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
482  R dec = R(1)/R(2), coeff2 = coeff * R(1.5);
483 
484  while (dec > R(1E-5) && res >= res0 * coeff) {
485  gmm::add(gmm::scaled(dr, -dec), pb.state_vector());
486  pb.compute_residual();
487  R res2 = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
488  if (res2 < res*R(0.95) || res2 >= res0 * coeff2) {
489  dec /= R(2); res = res2; coeff2 *= R(1.5);
490  } else {
491  gmm::add(gmm::scaled(dr, dec), pb.state_vector());
492  break;
493  }
494  }
495  dec *= R(2);
496 
497  nit++;
498  coeff = std::max(R(1.05), coeff*R(0.93));
499  bool near_end = (iter.get_iteration() > iter.get_maxiter()/2);
500  bool cut = (alpha < R(1)) && near_end;
501  if ((res > minres && nit > 4) || cut) {
502  stagn++;
503 
504  if ((stagn > 10 && alpha > alpha0 + R(5E-2)) || cut) {
505  alpha = (alpha + alpha0) / R(2);
506  if (near_end) alpha = R(1);
507  gmm::copy(Xi, pb.state_vector());
508  pb.compute_residual();
509  res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
510  nit = 0;
511  stagn = 0; coeff = R(2);
512  }
513  }
514  if (res < minres || (alpha == R(1) && nit < 10)) stagn = 0;
515  res0 = res;
516  if (nit < 5) minres = res0; else minres = std::min(minres, res0);
517 
518  if (iter.get_noisy())
519  cout << "step control [" << std::setw(8) << alpha0 << ","
520  << std::setw(8) << alpha << "," << std::setw(10) << dec << "]";
521  ++iter;
522  // crit = std::min(res / approx_eln,
523  // gmm::vect_norm1(dr) / std::max(1E-25,pb.state_norm()));
524  crit = res / approx_eln;
525  }
526 
527  if (iter.finished(crit)) {
528  if (iter.converged() && alpha < R(1)) {
529  R a = alpha;
530  alpha = std::min(R(1), alpha*R(3) - alpha0*R(2));
531  alpha0 = a;
532  gmm::copy(pb.state_vector(), Xi);
533  nit = 0; stagn = 0; coeff = R(2);
534  } else return;
535  }
536  }
537  }
538 
539 
540 
541  /* ***************************************************************** */
542  /* Classical Newton(-Raphson) algorithm. */
543  /* ***************************************************************** */
544 
545  template <typename PB>
546  void classical_Newton(PB &pb, gmm::iteration &iter)
547  {
548  typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T;
549  typedef typename gmm::number_traits<T>::magnitude_type R;
550  gmm::iteration iter_linsolv0 = iter;
551  iter_linsolv0.reduce_noisy();
552  iter_linsolv0.set_resmax(iter.get_resmax()/20.0);
553  iter_linsolv0.set_maxiter(10000); // arbitrary
554 
555  pb.compute_residual();
556  R approx_eln = pb.approx_external_load_norm();
557  if (approx_eln == R(0)) approx_eln = R(1);
558 
559  typename PB::VECTOR dr(gmm::vect_size(pb.residual()));
560 
561  scalar_type crit = pb.residual_norm() / approx_eln;
562  while (!iter.finished(crit)) {
563  gmm::iteration iter_linsolv = iter_linsolv0;
564 
565  int is_singular = 1;
566  while (is_singular) {
567  gmm::clear(dr);
568  iter_linsolv.init();
569  if (iter.get_noisy() > 1)
570  cout << "starting computing tangent matrix" << endl;
571  pb.compute_tangent_matrix();
572  if (iter.get_noisy() > 1)
573  cout << "starting linear solver" << endl;
574  pb.linear_solve(dr, iter_linsolv);
575  if (!iter_linsolv.converged()) {
576  is_singular++;
577  if (is_singular <= 4) {
578  if (iter.get_noisy())
579  cout << "Singular tangent matrix:"
580  " perturbation of the state vector." << endl;
581  pb.perturbation();
582  pb.compute_residual();
583  } else {
584  if (iter.get_noisy())
585  cout << "Singular tangent matrix: perturbation failed, aborting."
586  << endl;
587  return;
588  }
589  }
590  else is_singular = 0;
591  }
592 
593  if (iter.get_noisy() > 1) cout << "linear solver done" << endl;
594  R alpha = pb.line_search(dr, iter); //it is assumed that the linesearch
595  //executes a pb.compute_residual();
596  if (iter.get_noisy()) cout << "alpha = " << std::setw(6) << alpha << " ";
597  ++iter;
598  crit = std::min(pb.residual_norm() / approx_eln,
599  gmm::vect_norm1(dr) / std::max(1E-25, pb.state_norm()));
600  }
601  }
602 
603 
604 
605  //---------------------------------------------------------------------
606  // Default linear solver.
607  //---------------------------------------------------------------------
608 
609  typedef std::shared_ptr<abstract_linear_solver<model_real_sparse_matrix,
610  model_real_plain_vector> >
611  rmodel_plsolver_type;
612  typedef std::shared_ptr<abstract_linear_solver<model_complex_sparse_matrix,
613  model_complex_plain_vector> >
614  cmodel_plsolver_type;
615 
616  template<typename MATRIX, typename VECTOR>
617  std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>
618  default_linear_solver(const model &md) {
619 
620 #if GETFEM_PARA_LEVEL == 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
621  if (md.is_symmetric())
622  return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>();
623  else
624  return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
625 #elif GETFEM_PARA_LEVEL > 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
626  if (md.is_symmetric())
627  return std::make_shared
628  <linear_solver_distributed_mumps_sym<MATRIX, VECTOR>>();
629  else
630  return std::make_shared
631  <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
632 #else
633  size_type ndof = md.nb_dof(), max3d = 15000, dim = md.leading_dimension();
634 # ifdef GMM_USES_MUMPS
635  max3d = 250000;
636 # endif
637  if ((ndof<300000 && dim<=2) || (ndof<max3d && dim<=3) || (ndof<1000)) {
638 # ifdef GMM_USES_MUMPS
639  if (md.is_symmetric())
640  return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>();
641  else
642  return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
643 # else
644  return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>();
645 # endif
646  }
647  else {
648  if (md.is_coercive())
649  return std::make_shared
650  <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
651  else {
652  if (dim <= 2)
653  return std::make_shared
654  <linear_solver_gmres_preconditioned_ilut<MATRIX,VECTOR>>();
655  else
656  return std::make_shared
657  <linear_solver_gmres_preconditioned_ilu<MATRIX,VECTOR>>();
658  }
659  }
660 #endif
661  return std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>();
662  }
663 
664  template <typename MATRIX, typename VECTOR>
665  std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>
666  select_linear_solver(const model &md, const std::string &name) {
667  std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>> p;
668  if (bgeot::casecmp(name, "superlu") == 0)
669  return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>();
670  else if (bgeot::casecmp(name, "dense_lu") == 0)
671  return std::make_shared<linear_solver_dense_lu<MATRIX, VECTOR>>();
672  else if (bgeot::casecmp(name, "mumps") == 0) {
673 #ifdef GMM_USES_MUMPS
674 # if GETFEM_PARA_LEVEL <= 1
675  return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
676 # else
677  return std::make_shared
678  <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
679 # endif
680 #else
681  GMM_ASSERT1(false, "Mumps is not interfaced");
682 #endif
683  }
684  else if (bgeot::casecmp(name, "cg/ildlt") == 0)
685  return std::make_shared
686  <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
687  else if (bgeot::casecmp(name, "gmres/ilu") == 0)
688  return std::make_shared
689  <linear_solver_gmres_preconditioned_ilu<MATRIX, VECTOR>>();
690  else if (bgeot::casecmp(name, "gmres/ilut") == 0)
691  return std::make_shared
692  <linear_solver_gmres_preconditioned_ilut<MATRIX, VECTOR>>();
693  else if (bgeot::casecmp(name, "gmres/ilutp") == 0)
694  return std::make_shared
695  <linear_solver_gmres_preconditioned_ilutp<MATRIX, VECTOR>>();
696  else if (bgeot::casecmp(name, "auto") == 0)
697  return default_linear_solver<MATRIX, VECTOR>(md);
698  else
699  GMM_ASSERT1(false, "Unknown linear solver");
700  return std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>();
701  }
702 
703  inline rmodel_plsolver_type rselect_linear_solver(const model &md,
704  const std::string &name) {
705  return select_linear_solver<model_real_sparse_matrix,
706  model_real_plain_vector>(md, name);
707  }
708 
709  inline cmodel_plsolver_type cselect_linear_solver(const model &md,
710  const std::string &name) {
711  return select_linear_solver<model_complex_sparse_matrix,
712  model_complex_plain_vector>(md, name);
713  }
714 
715  //---------------------------------------------------------------------
716  // Standard solve.
717  //---------------------------------------------------------------------
718 
719 
720  /** A default solver for the model brick system.
721  Of course it could be not very well suited for a particular
722  problem, so it could be copied and adapted to change solvers,
723  add a special traitement on the problem, etc ... This is in
724  fact a model for your own solver.
725 
726  For small problems, a direct solver is used
727  (getfem::SuperLU_solve), for larger problems, a conjugate
728  gradient gmm::cg (if the problem is coercive) or a gmm::gmres is
729  used (preconditioned with an incomplete factorization).
730 
731  When MPI/METIS is enabled, a partition is done via METIS, and a parallel
732  solver can be used.
733 
734  Note that it is possible to disable some variables
735  (see the md.disable_variable(varname) method) in order to
736  solve the problem only with respect to a subset of variables (the
737  disabled variables are the considered as data) for instance to
738  replace the global Newton strategy with a fixed point one.
739 
740  @ingroup bricks
741  */
742  void standard_solve(model &md, gmm::iteration &iter,
743  rmodel_plsolver_type lsolver,
744  abstract_newton_line_search &ls);
745 
746  void standard_solve(model &md, gmm::iteration &iter,
747  cmodel_plsolver_type lsolver,
748  abstract_newton_line_search &ls);
749 
750  void standard_solve(model &md, gmm::iteration &iter,
751  rmodel_plsolver_type lsolver);
752 
753  void standard_solve(model &md, gmm::iteration &iter,
754  cmodel_plsolver_type lsolver);
755 
756  void standard_solve(model &md, gmm::iteration &iter);
757 
758 } /* end of namespace getfem. */
759 
760 
761 #endif /* GETFEM_MODEL_SOLVERS_H__ */
gmm::vect_dist1
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist1(const V1 &v1, const V2 &v2)
1-distance between two vectors
Definition: gmm_blas.h:658
gmm::ildlt_precond
Incomplete Level 0 LDLT Preconditioner.
Definition: gmm_precond_ildlt.h:93
gmm::ilu_precond
Incomplete LU without fill-in Preconditioner.
Definition: gmm_precond_ilu.h:86
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm_MUMPS_interface.h
Interface with MUMPS (LU direct solver for sparse matrices).
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
gmm_iter_solvers.h
Include standard gmm iterative solvers (cg, gmres, ...)
getfem_models.h
Model representation in Getfem.
gmm::ilutp_precond
ILUTP: Incomplete LU with threshold and K fill-in Preconditioner and column pivoting.
Definition: gmm_precond_ilutp.h:57
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
gmm_dense_qr.h
Dense QR factorization.
getfem::standard_solve
void standard_solve(model &md, gmm::iteration &iter, rmodel_plsolver_type lsolver, abstract_newton_line_search &ls)
A default solver for the model brick system.
Definition: getfem_model_solvers.cc:410
gmm::ilut_precond
Incomplete LU with threshold and K fill-in Preconditioner.
Definition: gmm_precond_ilut.h:102
gmm_iter.h
Iteration object.