rheolef  6.6
solver.cc
Go to the documentation of this file.
1 #include "rheolef/solver.h"
2 #include "rheolef/eye.h"
3 #include "solver_wrapper.h"
4 #include "solver_ldlt_builtin.h"
5 #include "solver_cholmod.h"
6 #include "solver_umfpack.h"
7 #include "solver_mumps.h"
10 #include "solver_pastix.h"
11 #include <typeinfo>
12 
13 namespace rheolef {
14 
15 template <class T, class M>
16 void
18 {
19  fatal_macro ("undefined set_preconditionner for this solver");
20 }
21 template <class T, class M>
22 void
24 {
25  check_macro (_ptr != 0, "solver::set_preconditionner(): uninitialized solver");
26  _ptr->set_preconditionner (sa);
27 }
28 
29 template <class T, class M>
32 {
33  std::string type = typeid(*this).name();
34  warning_macro ("undefined computation of the determinant for this solver="<<type);
35  return determinant_type();
36 }
37 template <class T, class M>
40 {
41  check_macro (_ptr != 0, "solver::det(): uninitialized solver");
42  return _ptr->det();
43 }
44 
45 template<class T, class M>
47  : _ptr(0)
48 {
49 }
50 template<class T, class M>
52  : _ptr(0)
53 {
55  _ptr = new_macro (rep(a,opt));
56 }
57 template<class T, class M>
58 const solver_option_type&
60 {
61  check_macro (_ptr != 0, "solver::option(): uninitialized solver");
62  return _ptr -> option();
63 }
64 template<class T, class M>
65 void
67 {
68  typedef eye_rep<T,M> rep;
69  if (_ptr) delete_macro (_ptr);
70  _ptr = new_macro (rep());
71 }
72 template<class T, class M>
73 void
75 {
76 #ifdef _RHEOLEF_HAVE_MUMPS
77  //bool prefer_cholmod = (a.dis_ext_nnz() == 0) && a.is_symmetric() && a.is_definite_positive(); // when sdp and diagonal per process
78  //bool prefer_umfpack = opt.force_seq || (a.dis_ext_nnz() == 0); // when diagonal per process
79  bool prefer_cholmod = false; // prefer mumps that seems faster than cholmod, even in seq case
80  bool prefer_umfpack = opt.force_seq // idem, but when building schur complement, mumps is still buggy
81  || ((a.dis_ext_nnz() == 0) && opt.prefered_library == "umfpack"); // explicitely prefer (more robust on ill-conditioned systems)
82 #else // _RHEOLEF_HAVE_MUMPS
83  bool prefer_cholmod = (a.dis_ext_nnz() == 0) && a.is_symmetric() && a.is_definite_positive(); // when sdp and diagonal per process
84  bool prefer_umfpack = (a.dis_ext_nnz() == 0) // when (unsym or (sym & undef)) & diag per proc
85  || (opt.prefered_library == "umfpack"); // explicitely prefer (more robust on ill-conditioned systems)
86 #endif // _RHEOLEF_HAVE_MUMPS
87 
88 #ifdef _RHEOLEF_HAVE_CHOLMOD
89  if (prefer_cholmod) {
90  // diagonal bloc per process => use cholmod
91  typedef solver_cholmod_rep<T,M> rep;
92  if (_ptr) delete_macro (_ptr);
93  _ptr = new_macro (rep(a,opt));
94  return;
95  }
96 #endif // _RHEOLEF_HAVE_CHOLMOD
97 
98 #ifdef _RHEOLEF_HAVE_UMFPACK
99  if (prefer_umfpack) {
100  // diagonal bloc per process => use umfpack
101  typedef solver_umfpack_rep<T,M> rep;
102  if (_ptr) delete_macro (_ptr);
103  _ptr = new_macro (rep(a,opt));
104  return;
105  }
106 #endif // _RHEOLEF_HAVE_UMFPACK
107 
108 #ifdef _RHEOLEF_HAVE_MUMPS
109  typedef solver_mumps_rep<T,M> rep;
110 #elif defined(_RHEOLEF_HAVE_PASTIX)
111  typedef solver_pastix_rep<T,M> rep;
112  opt.iterative = 0;
113 #else // ! _RHEOLEF_HAVE_TRILINOS && !_RHEOLEF_HAVE_PASTIX
114  typedef solver_ldlt_builtin_rep<T,M> rep;
115 #endif // ! _RHEOLEF_HAVE_MUMPS
116 
117  if (_ptr) delete_macro (_ptr);
118  _ptr = new_macro (rep(a,opt));
119 }
120 template<class T, class M>
121 void
123 {
124 #ifdef _RHEOLEF_HAVE_TRILINOS
125  typedef solver_trilinos_ifpack_rep<T,M> rep;
126 #elif defined(_RHEOLEF_HAVE_PASTIX)
127  typedef solver_pastix_rep<T,M> rep;
128  opt.iterative = 1;
129 #else // ! _RHEOLEF_HAVE_TRILINOS && !_RHEOLEF_HAVE_PASTIX
131 #endif
132  if (_ptr) delete_macro (_ptr);
133  _ptr = new_macro (rep(a,opt));
134 }
135 template<class T, class M>
137 {
138  if (_ptr) delete_macro (_ptr);
139  _ptr = 0;
140 }
141 template<class T, class M>
142 void
144 {
145  _ptr->update_values (a);
146 }
147 template<class T, class M>
148 vec<T,M>
150 {
151  return _ptr->solve (b);
152 }
153 template<class T, class M>
154 vec<T,M>
156 {
157  return _ptr->trans_solve (b);
158 }
159 #define _RHEOLEF_instanciation(T,M) \
160 template class solver_abstract_rep<T,M>; \
161 template class solver_rep<T,M>;
162 
163 _RHEOLEF_instanciation(Float,sequential)
164 #ifdef _RHEOLEF_HAVE_MPI
165 _RHEOLEF_instanciation(Float,distributed)
166 #endif // _RHEOLEF_HAVE_MPI
167 #undef _RHEOLEF_instanciation
168 
169 } // namespace rheolef
vec< T, M > trans_solve(const vec< T, M > &b) const
Definition: solver.cc:155
#define new_macro(obj)
Definition: compiler.h:364
#define fatal_macro(message)
Definition: compiler.h:98
solver_option_type - options for direct or interative solvers
solver_abstract_rep< T, M > rep
Definition: solver.h:49
#define delete_macro(ptr)
Definition: compiler.h:366
void update_values(const csr< T, M > &a)
Definition: solver.cc:143
eye - the identity matrix
Definition: eye.h:22
virtual determinant_type det() const
Definition: solver.cc:31
irheostream, orheostream - large data streams
Definition: compiler.h:7
const solver_option_type & option() const
Definition: solver.cc:59
#define warning_macro(message)
Definition: compiler.h:102
determinant_type det() const
Definition: solver.cc:39
#define check_macro(ok_condition, message)
Definition: compiler.h:104
double Float
Definition: compiler.h:177
void build_eye()
Definition: solver.cc:66
void build_ilu(const csr< T, M > &a, const solver_option_type &opt)
Definition: solver.cc:122
void build_lu(const csr< T, M > &a, const solver_option_type &opt)
Definition: solver.cc:74
vec - vector in distributed environment
Definition: vec.h:44
solver - direct or interative solver interface
Definition: solver.h:8
#define _RHEOLEF_instanciation(T, M)
Definition: solver.cc:159
vec< T, M > solve(const vec< T, M > &b) const
Definition: solver.cc:149
virtual void set_preconditionner(const solver_basic< T, M > &)
Definition: solver.cc:17
void set_preconditionner(const solver_basic< T, M > &)
Definition: solver.cc:23
csr - compressed sparse row matrix
Definition: asr.h:8