GetFEM++  5.2
getfem_model_solvers.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2004-2017 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"
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  virtual void operator ()(const MAT &, VECT &, const VECT &,
83  gmm::iteration &) const = 0;
84  virtual ~abstract_linear_solver() {}
85  };
86 
87  template <typename MAT, typename VECT>
88  struct linear_solver_cg_preconditioned_ildlt
89  : public abstract_linear_solver<MAT, VECT> {
90  void operator ()(const MAT &M, VECT &x, const VECT &b,
91  gmm::iteration &iter) const {
93  gmm::cg(M, x, b, P, iter);
94  if (!iter.converged()) GMM_WARNING2("cg did not converge!");
95  }
96  };
97 
98  template <typename MAT, typename VECT>
99  struct linear_solver_gmres_preconditioned_ilu
100  : public abstract_linear_solver<MAT, VECT> {
101  void operator ()(const MAT &M, VECT &x, const VECT &b,
102  gmm::iteration &iter) const {
104  gmm::gmres(M, x, b, P, 500, iter);
105  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
106  }
107  };
108 
109  template <typename MAT, typename VECT>
110  struct linear_solver_gmres_unpreconditioned
111  : public abstract_linear_solver<MAT, VECT> {
112  void operator ()(const MAT &M, VECT &x, const VECT &b,
113  gmm::iteration &iter) const {
114  gmm::identity_matrix P;
115  gmm::gmres(M, x, b, P, 500, iter);
116  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
117  }
118  };
119 
120  template <typename MAT, typename VECT>
121  struct linear_solver_gmres_preconditioned_ilut
122  : public abstract_linear_solver<MAT, VECT> {
123  void operator ()(const MAT &M, VECT &x, const VECT &b,
124  gmm::iteration &iter) const {
125  gmm::ilut_precond<MAT> P(M, 40, 1E-7);
126  gmm::gmres(M, x, b, P, 500, iter);
127  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
128  }
129  };
130 
131  template <typename MAT, typename VECT>
132  struct linear_solver_gmres_preconditioned_ilutp
133  : public abstract_linear_solver<MAT, VECT> {
134  void operator ()(const MAT &M, VECT &x, const VECT &b,
135  gmm::iteration &iter) const {
136  gmm::ilutp_precond<MAT> P(M, 20, 1E-7);
137  gmm::gmres(M, x, b, P, 500, iter);
138  if (!iter.converged()) GMM_WARNING2("gmres did not converge!");
139  }
140  };
141 
142  template <typename MAT, typename VECT>
143  struct linear_solver_superlu
144  : public abstract_linear_solver<MAT, VECT> {
145  void operator ()(const MAT &M, VECT &x, const VECT &b,
146  gmm::iteration &iter) const {
147  double rcond;
148  /*gmm::HarwellBoeing_IO::write("test.hb", M);
149  std::fstream f("bbb", std::ios::out);
150  for (unsigned i=0; i < gmm::vect_size(b); ++i) f << b[i] << "\n";*/
151  int info = SuperLU_solve(M, x, b, rcond);
152  iter.enforce_converged(info == 0);
153  if (iter.get_noisy()) cout << "condition number: " << 1.0/rcond<< endl;
154  }
155  };
156 
157  template <typename MAT, typename VECT>
158  struct linear_solver_dense_lu : public abstract_linear_solver<MAT, VECT> {
159  void operator ()(const MAT &M, VECT &x, const VECT &b,
160  gmm::iteration &iter) const {
161  typedef typename gmm::linalg_traits<MAT>::value_type T;
162  gmm::dense_matrix<T> MM(gmm::mat_nrows(M),gmm::mat_ncols(M));
163  gmm::copy(M, MM);
164  gmm::lu_solve(MM, x, b);
165  iter.enforce_converged(true);
166  }
167  };
168 
169 #ifdef GMM_USES_MUMPS
170  template <typename MAT, typename VECT>
171  struct linear_solver_mumps : public abstract_linear_solver<MAT, VECT> {
172  void operator ()(const MAT &M, VECT &x, const VECT &b,
173  gmm::iteration &iter) const {
174  bool ok = gmm::MUMPS_solve(M, x, b, false);
175  iter.enforce_converged(ok);
176  }
177  };
178  template <typename MAT, typename VECT>
179  struct linear_solver_mumps_sym : public abstract_linear_solver<MAT, VECT> {
180  void operator ()(const MAT &M, VECT &x, const VECT &b,
181  gmm::iteration &iter) const {
182  bool ok = gmm::MUMPS_solve(M, x, b, true);
183  iter.enforce_converged(ok);
184  }
185  };
186 #endif
187 
188 #if GETFEM_PARA_LEVEL > 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
189  template <typename MAT, typename VECT>
190  struct linear_solver_distributed_mumps
191  : public abstract_linear_solver<MAT, VECT> {
192  void operator ()(const MAT &M, VECT &x, const VECT &b,
193  gmm::iteration &iter) const {
194  double tt_ref=MPI_Wtime();
195  bool ok = MUMPS_distributed_matrix_solve(M, x, b, false);
196  iter.enforce_converged(ok);
197  if (MPI_IS_MASTER()) cout<<"UNSYMMETRIC MUMPS time "<< MPI_Wtime() - tt_ref<<endl;
198  }
199  };
200 
201  template <typename MAT, typename VECT>
202  struct linear_solver_distributed_mumps_sym
203  : public abstract_linear_solver<MAT, VECT> {
204  void operator ()(const MAT &M, VECT &x, const VECT &b,
205  gmm::iteration &iter) const {
206  double tt_ref=MPI_Wtime();
207  bool ok = MUMPS_distributed_matrix_solve(M, x, b, true);
208  iter.enforce_converged(ok);
209  if (MPI_IS_MASTER()) cout<<"SYMMETRIC MUMPS time "<< MPI_Wtime() - tt_ref<<endl;
210  }
211  };
212 #endif
213 
214 
215  /* ***************************************************************** */
216  /* Newton Line search definition */
217  /* ***************************************************************** */
218 
219  struct abstract_newton_line_search {
220  double conv_alpha, conv_r;
221  size_t it, itmax, glob_it;
222  // size_t tot_it;
223  virtual void init_search(double r, size_t git, double R0 = 0.0) = 0;
224  virtual double next_try(void) = 0;
225  virtual bool is_converged(double, double R1 = 0.0) = 0;
226  virtual double converged_value(void) {
227  // tot_it += it; cout << "tot_it = " << tot_it << endl; it = 0;
228  return conv_alpha;
229  };
230  virtual double converged_residual(void) { return conv_r; };
231  virtual ~abstract_newton_line_search() { }
232  };
233 
234 
235  struct quadratic_newton_line_search : public abstract_newton_line_search {
236  double R0_, R1_;
237 
238  double alpha, alpha_mult, first_res, alpha_max_ratio, alpha_min;
239  virtual void init_search(double r, size_t git, double R0 = 0.0) {
240  GMM_ASSERT1(R0 != 0.0, "You have to specify R0");
241  glob_it = git;
242  conv_alpha = alpha = double(1); conv_r = first_res = r; it = 0;
243  R0_ = R0;
244  }
245  virtual double next_try(void) {
246  ++it;
247  if (it == 1) return double(1);
248  GMM_ASSERT1(R1_ != 0.0, "You have to specify R1");
249  double a = R0_ / R1_;
250  return (a < 0) ? (a*0.5 + sqrt(a*a*0.25-a)) : a*0.5;
251  }
252  virtual bool is_converged(double r, double R1 = 0.0) {
253  conv_r = r;
254  R1_ = R1; return ((gmm::abs(R1_) < gmm::abs(R0_*0.5)) || it >= itmax);
255  }
256  quadratic_newton_line_search(size_t imax = size_t(-1)) { itmax = imax; }
257  };
258 
259 
260  struct simplest_newton_line_search : public abstract_newton_line_search {
261  double alpha, alpha_mult, first_res, alpha_max_ratio, alpha_min,
262  alpha_threshold_res;
263  virtual void init_search(double r, size_t git, double = 0.0) {
264  glob_it = git;
265  conv_alpha = alpha = double(1); conv_r = first_res = r; it = 0;
266  }
267  virtual double next_try(void)
268  { conv_alpha = alpha; alpha *= alpha_mult; ++it; return conv_alpha; }
269  virtual bool is_converged(double r, double = 0.0) {
270  conv_r = r;
271  return ((it <= 1 && r < first_res)
272  || (r <= first_res * alpha_max_ratio && r <= alpha_threshold_res)
273  || (conv_alpha <= alpha_min && r < first_res * 1e5)
274  || it >= itmax);
275  }
276  simplest_newton_line_search
277  (size_t imax = size_t(-1), double a_max_ratio = 6.0/5.0,
278  double a_min = 1.0/1000.0, double a_mult = 3.0/5.0,
279  double a_threshold_res = 1e50)
280  : alpha_mult(a_mult), alpha_max_ratio(a_max_ratio), alpha_min(a_min),
281  alpha_threshold_res(a_threshold_res)
282  { itmax = imax; }
283  };
284 
285  struct default_newton_line_search : public abstract_newton_line_search {
286  // This line search try to detect where is the minimum, dividing the step
287  // by a factor two each time.
288  // - it stops if the residual is less than the previous residual
289  // times alpha_min_ratio (= 0.9).
290  // - if the minimal step is reached with a residual greater than
291  // the previous residual times alpha_min_ratio then it decides
292  // between two options :
293  // - return with the step corresponding to the smallest residual
294  // - return with a greater residual corresponding to a residual
295  // less than the previous residual times alpha_max_ratio.
296  // the decision is taken regarding the previous iterations.
297  // - in order to shorten the line search, the process stops when
298  // the residual increases three times consecutively.
299  // possible improvment : detect the curvature at the origin
300  // (only one evaluation) and take it into account.
301  // Fitted to some experiments in contrib/tests_newton
302 
303  double alpha, alpha_old, alpha_mult, first_res, alpha_max_ratio;
304  double alpha_min_ratio, alpha_min;
305  size_type count, count_pat;
306  bool max_ratio_reached;
307  double alpha_max_ratio_reached, r_max_ratio_reached;
308  size_type it_max_ratio_reached;
309 
310  virtual void init_search(double r, size_t git, double = 0.0);
311  virtual double next_try(void);
312  virtual bool is_converged(double r, double = 0.0);
313  default_newton_line_search(void) { count_pat = 0; }
314  };
315 
316  /* the former default_newton_line_search */
317  struct basic_newton_line_search : public abstract_newton_line_search {
318  double alpha, alpha_mult, first_res, alpha_max_ratio;
319  double alpha_min, prev_res, alpha_max_augment;
320  virtual void init_search(double r, size_t git, double = 0.0) {
321  glob_it = git;
322  conv_alpha = alpha = double(1);
323  prev_res = conv_r = first_res = r; it = 0;
324  }
325  virtual double next_try(void)
326  { conv_alpha = alpha; alpha *= alpha_mult; ++it; return conv_alpha; }
327  virtual bool is_converged(double r, double = 0.0) {
328  if (glob_it == 0 || (r < first_res / double(2))
329  || (conv_alpha <= alpha_min && r < first_res * alpha_max_augment)
330  || it >= itmax)
331  { conv_r = r; return true; }
332  if (it > 1 && r > prev_res && prev_res < alpha_max_ratio * first_res)
333  return true;
334  prev_res = conv_r = r;
335  return false;
336  }
337  basic_newton_line_search
338  (size_t imax = size_t(-1),
339  double a_max_ratio = 5.0/3.0,
340  double a_min = 1.0/1000.0, double a_mult = 3.0/5.0, double a_augm = 2.0)
341  : alpha_mult(a_mult), alpha_max_ratio(a_max_ratio),
342  alpha_min(a_min), alpha_max_augment(a_augm) { itmax = imax; }
343  };
344 
345 
346  struct systematic_newton_line_search : public abstract_newton_line_search {
347  double alpha, alpha_mult, first_res;
348  double alpha_min, prev_res;
349  bool first;
350  virtual void init_search(double r, size_t git, double = 0.0) {
351  glob_it = git;
352  conv_alpha = alpha = double(1);
353  prev_res = conv_r = first_res = r; it = 0; first = true;
354  }
355  virtual double next_try(void)
356  { double a = alpha; alpha *= alpha_mult; ++it; return a; }
357  virtual bool is_converged(double r, double = 0.0) {
358  // cout << "a = " << alpha / alpha_mult << " r = " << r << endl;
359  if (r < conv_r || first)
360  { conv_r = r; conv_alpha = alpha / alpha_mult; first = false; }
361  if ((alpha <= alpha_min*alpha_mult) || it >= itmax) return true;
362  return false;
363  }
364  systematic_newton_line_search
365  (size_t imax = size_t(-1),
366  double a_min = 1.0/10000.0, double a_mult = 3.0/5.0)
367  : alpha_mult(a_mult), alpha_min(a_min) { itmax = imax; }
368  };
369 
370 
371 
372 
373 
374  /* ***************************************************************** */
375  /* Newton algorithm. */
376  /* ***************************************************************** */
377 
378  template <typename PB>
379  void classical_Newton(PB &pb, gmm::iteration &iter,
380  const abstract_linear_solver<typename PB::MATRIX,
381  typename PB::VECTOR> &linear_solver) {
382  typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T;
383  typedef typename gmm::number_traits<T>::magnitude_type R;
384  gmm::iteration iter_linsolv0 = iter;
385  iter_linsolv0.reduce_noisy();
386  iter_linsolv0.set_resmax(iter.get_resmax()/20.0);
387  iter_linsolv0.set_maxiter(10000); // arbitrary
388 
389  pb.compute_residual();
390 
391  typename PB::VECTOR dr(gmm::vect_size(pb.residual()));
392  typename PB::VECTOR b(gmm::vect_size(pb.residual()));
393 
394  scalar_type crit = pb.residual_norm()
395  / std::max(1E-25, pb.approx_external_load_norm());
396  while (!iter.finished(crit)) {
397  gmm::iteration iter_linsolv = iter_linsolv0;
398  if (iter.get_noisy() > 1)
399  cout << "starting computing tangent matrix" << endl;
400 
401  int is_singular = 1;
402  while (is_singular) {
403  pb.compute_tangent_matrix();
404  gmm::clear(dr);
405  gmm::copy(gmm::scaled(pb.residual(), pb.scale_residual()), b);
406  if (iter.get_noisy() > 1) cout << "starting linear solver" << endl;
407  iter_linsolv.init();
408  linear_solver(pb.tangent_matrix(), dr, b, iter_linsolv);
409  if (!iter_linsolv.converged()) {
410  is_singular++;
411  if (is_singular <= 4) {
412  if (iter.get_noisy())
413  cout << "Singular tangent matrix:"
414  " perturbation of the state vector." << endl;
415  pb.perturbation();
416  pb.compute_residual();
417  } else {
418  if (iter.get_noisy())
419  cout << "Singular tangent matrix: perturbation failed, aborting."
420  << endl;
421  return;
422  }
423  }
424  else is_singular = 0;
425  }
426 
427  if (iter.get_noisy() > 1) cout << "linear solver done" << endl;
428  R alpha = pb.line_search(dr, iter); //it is assumed that the linesearch
429  //executes a pb.compute_residual();
430  if (iter.get_noisy()) cout << "alpha = " << std::setw(6) << alpha << " ";
431  ++iter;
432  crit = std::min(pb.residual_norm()
433  / std::max(1E-25, pb.approx_external_load_norm()),
434  gmm::vect_norm1(dr) / std::max(1E-25, pb.state_norm()));
435  }
436  }
437 
438 
439  /* ***************************************************************** */
440  /* Intermediary structure for Newton algorithms with getfem::model. */
441  /* ***************************************************************** */
442 
443  #define TRACE_SOL 0
444 
445  template <typename MAT, typename VEC>
446  struct model_pb {
447 
448  typedef MAT MATRIX;
449  typedef VEC VECTOR;
450  typedef typename gmm::linalg_traits<VECTOR>::value_type T;
451  typedef typename gmm::number_traits<T>::magnitude_type R;
452 
453  model &md;
454  abstract_newton_line_search &ls;
455  VECTOR stateinit, &state;
456  const VECTOR &rhs;
457  const MATRIX &K;
458 
459  void compute_tangent_matrix(void) {
460  md.to_variables(state);
461  md.assembly(model::BUILD_MATRIX);
462  }
463 
464  const MATRIX &tangent_matrix(void) { return K; }
465 
466  inline T scale_residual(void) const { return T(1); }
467 
468  void compute_residual(void) {
469  md.to_variables(state);
470  md.assembly(model::BUILD_RHS);
471  }
472 
473  void perturbation(void) {
474  R res = gmm::vect_norm2(state), ampl = std::max(res*R(1E-20), R(1E-50));
475  std::vector<R> V(gmm::vect_size(state));
476  gmm::fill_random(V);
477  gmm::add(gmm::scaled(V, ampl), state);
478  }
479 
480  const VECTOR &residual(void) { return rhs; }
481 
482  R state_norm(void) const
483  { return gmm::vect_norm1(state); }
484 
485  R approx_external_load_norm(void)
486  { return md.approx_external_load(); }
487 
488  R residual_norm(void) { // A norm1 seems to be better than a norm2
489  // at least for contact problems.
490  return gmm::vect_norm1(rhs);
491  }
492 
493  R compute_res(bool comp = true) {
494  if (comp) compute_residual();
495  return residual_norm();
496  }
497 
498 
499  R line_search(VECTOR &dr, const gmm::iteration &iter) {
500  size_type nit = 0;
501  gmm::resize(stateinit, md.nb_dof());
502  gmm::copy(state, stateinit);
503  R alpha(1), res, /* res_init, */ R0;
504 
505  /* res_init = */ res = compute_res(false);
506  // cout << "first residual = " << residual() << endl << endl;
507  R0 = gmm::real(gmm::vect_sp(dr, rhs));
508 
509 #if TRACE_SOL
510  static int trace_number = 0;
511  int trace_iter = 0;
512  {
513  std::stringstream trace_name;
514  trace_name << "line_search_state" << std::setfill('0')
515  << std::setw(3) << trace_number << "_000_init";
516  gmm::vecsave(trace_name.str(),stateinit);
517  }
518  trace_number++;
519 #endif
520 
521  ls.init_search(res, iter.get_iteration(), R0);
522  do {
523  alpha = ls.next_try();
524  gmm::add(stateinit, gmm::scaled(dr, alpha), state);
525 #if TRACE_SOL
526  {
527  trace_iter++;
528  std::stringstream trace_name;
529  trace_name << "line_search_state" << std::setfill('0')
530  << std::setw(3) << trace_number << "_"
531  << std::setfill('0') << std::setw(3) << trace_iter;
532  gmm::vecsave(trace_name.str(), state);
533  }
534 #endif
535  res = compute_res();
536  // cout << "residual = " << residual() << endl << endl;
537  R0 = gmm::real(gmm::vect_sp(dr, rhs));
538 
539  ++ nit;
540  } while (!ls.is_converged(res, R0));
541 
542  if (alpha != ls.converged_value()) {
543  alpha = ls.converged_value();
544  gmm::add(stateinit, gmm::scaled(dr, alpha), state);
545  res = ls.converged_residual();
546  compute_residual();
547  }
548 
549  return alpha;
550  }
551 
552  model_pb(model &m, abstract_newton_line_search &ls_, VECTOR &st,
553  const VECTOR &rhs_, const MATRIX &K_)
554  : md(m), ls(ls_), state(st), rhs(rhs_), K(K_) {}
555 
556  };
557 
558  //---------------------------------------------------------------------
559  // Default linear solver.
560  //---------------------------------------------------------------------
561 
562  typedef abstract_linear_solver<model_real_sparse_matrix,
563  model_real_plain_vector> rmodel_linear_solver;
564  typedef std::shared_ptr<rmodel_linear_solver> rmodel_plsolver_type;
565  typedef abstract_linear_solver<model_complex_sparse_matrix,
566  model_complex_plain_vector>
567  cmodel_linear_solver;
568  typedef std::shared_ptr<cmodel_linear_solver> cmodel_plsolver_type;
569 
570 
571  template<typename MATRIX, typename VECTOR>
572  std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>
573  default_linear_solver(const model &md) {
574 
575 #if GETFEM_PARA_LEVEL == 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
576  if (md.is_symmetric())
577  return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>();
578  else
579  return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
580 #elif GETFEM_PARA_LEVEL > 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
581  if (md.is_symmetric())
582  return std::make_shared
583  <linear_solver_distributed_mumps_sym<MATRIX, VECTOR>>();
584  else
585  return std::make_shared
586  <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
587 #else
588  size_type ndof = md.nb_dof(), max3d = 15000, dim = md.leading_dimension();
589 # ifdef GMM_USES_MUMPS
590  max3d = 250000;
591 # endif
592  if ((ndof<300000 && dim<=2) || (ndof<max3d && dim<=3) || (ndof<1000)) {
593 # ifdef GMM_USES_MUMPS
594  if (md.is_symmetric())
595  return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>();
596  else
597  return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
598 # else
599  return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>();
600 # endif
601  }
602  else {
603  if (md.is_coercive())
604  return std::make_shared
605  <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
606  else {
607  if (dim <= 2)
608  return std::make_shared
609  <linear_solver_gmres_preconditioned_ilut<MATRIX,VECTOR>>();
610  else
611  return std::make_shared
612  <linear_solver_gmres_preconditioned_ilu<MATRIX,VECTOR>>();
613  }
614  }
615 #endif
616  return std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>();
617  }
618 
619  template <typename MATRIX, typename VECTOR>
620  std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>
621  select_linear_solver(const model &md, const std::string &name) {
622  std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>> p;
623  if (bgeot::casecmp(name, "superlu") == 0)
624  return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>();
625  else if (bgeot::casecmp(name, "dense_lu") == 0)
626  return std::make_shared<linear_solver_dense_lu<MATRIX, VECTOR>>();
627  else if (bgeot::casecmp(name, "mumps") == 0) {
628 #ifdef GMM_USES_MUMPS
629 # if GETFEM_PARA_LEVEL <= 1
630  return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
631 # else
632  return std::make_shared
633  <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
634 # endif
635 #else
636  GMM_ASSERT1(false, "Mumps is not interfaced");
637 #endif
638  }
639  else if (bgeot::casecmp(name, "cg/ildlt") == 0)
640  return std::make_shared
641  <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
642  else if (bgeot::casecmp(name, "gmres/ilu") == 0)
643  return std::make_shared
644  <linear_solver_gmres_preconditioned_ilu<MATRIX, VECTOR>>();
645  else if (bgeot::casecmp(name, "gmres/ilut") == 0)
646  return std::make_shared
647  <linear_solver_gmres_preconditioned_ilut<MATRIX, VECTOR>>();
648  else if (bgeot::casecmp(name, "gmres/ilutp") == 0)
649  return std::make_shared
650  <linear_solver_gmres_preconditioned_ilutp<MATRIX, VECTOR>>();
651  else if (bgeot::casecmp(name, "auto") == 0)
652  return default_linear_solver<MATRIX, VECTOR>(md);
653  else
654  GMM_ASSERT1(false, "Unknown linear solver");
655  return std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>();
656  }
657 
658  inline rmodel_plsolver_type rselect_linear_solver(const model &md,
659  const std::string &name) {
660  return select_linear_solver<model_real_sparse_matrix,
661  model_real_plain_vector>(md, name);
662  }
663 
664  inline cmodel_plsolver_type cselect_linear_solver(const model &md,
665  const std::string &name) {
666  return select_linear_solver<model_complex_sparse_matrix,
667  model_complex_plain_vector>(md, name);
668  }
669 
670  //---------------------------------------------------------------------
671  // Standard solve.
672  //---------------------------------------------------------------------
673 
674 
675  /** A default solver for the model brick system.
676  Of course it could be not very well suited for a particular
677  problem, so it could be copied and adapted to change solvers,
678  add a special traitement on the problem, etc ... This is in
679  fact a model for your own solver.
680 
681  For small problems, a direct solver is used
682  (getfem::SuperLU_solve), for larger problems, a conjugate
683  gradient gmm::cg (if the problem is coercive) or a gmm::gmres is
684  used (preconditioned with an incomplete factorization).
685 
686  When MPI/METIS is enabled, a partition is done via METIS, and a parallel
687  solver can be used.
688 
689  Note that it is possible to disable some variables
690  (see the md.disable_variable(varname) method) in order to
691  solve the problem only with respect to a subset of variables (the
692  disabled variables are the considered as data) for instance to
693  replace the global Newton strategy with a fixed point one.
694 
695  @ingroup bricks
696  */
697  void standard_solve(model &md, gmm::iteration &iter,
698  rmodel_plsolver_type lsolver,
699  abstract_newton_line_search &ls);
700 
701  void standard_solve(model &md, gmm::iteration &iter,
702  cmodel_plsolver_type lsolver,
703  abstract_newton_line_search &ls);
704 
705  void standard_solve(model &md, gmm::iteration &iter,
706  rmodel_plsolver_type lsolver);
707 
708  void standard_solve(model &md, gmm::iteration &iter,
709  cmodel_plsolver_type lsolver);
710 
711  void standard_solve(model &md, gmm::iteration &iter);
712 
713 } /* end of namespace getfem. */
714 
715 
716 #endif /* GETFEM_MODEL_SOLVERS_H__ */
Incomplete LU without fill-in Preconditioner.
The Iteration object calculates whether the solution has reached the desired accuracy, or whether the maximum number of iterations has been reached.
Definition: gmm_iter.h:53
Dense QR factorization.
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.
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
Model representation in Getfem.
Interface with MUMPS (LU direct solver for sparse matrices).
Iteration object.
Include standard gmm iterative solvers (cg, gmres, ...)
ILUTP: Incomplete LU with threshold and K fill-in Preconditioner and column pivoting.
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:46
Incomplete LU with threshold and K fill-in Preconditioner.
Incomplete Level 0 LDLT Preconditioner.