rheolef  7.0
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&
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) &&
82  (opt.prefered_library == "umfpack"
83 # if !(defined( _RHEOLEF_HAVE_MUMPS_WITH_METIS) || \
84  defined( _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS) || \
85  defined( _RHEOLEF_HAVE_MUMPS_WITH_SCOTCH) || \
87  || opt.prefered_library == ""
88 #endif // prefer umfpack when mumps has poor ordering
89  ));
90 #else // !_RHEOLEF_HAVE_MUMPS
91  bool prefer_cholmod = (a.dis_ext_nnz() == 0) && a.is_symmetric() && a.is_definite_positive(); // when sdp and diagonal per process
92  bool prefer_umfpack = (a.dis_ext_nnz() == 0) // when (unsym or (sym & undef)) & diag per proc
93  || (opt.prefered_library == "umfpack"); // explicitely prefer (more robust on ill-conditioned systems)
94 #endif // ! _RHEOLEF_HAVE_MUMPS
95 
96 #ifdef _RHEOLEF_HAVE_CHOLMOD
97  if (prefer_cholmod) {
98  // diagonal bloc per process => use cholmod
99  typedef solver_cholmod_rep<T,M> rep;
100  if (_ptr) delete_macro (_ptr);
101  _ptr = new_macro (rep(a,opt));
102  return;
103  }
104 #endif // _RHEOLEF_HAVE_CHOLMOD
105 
106 #ifdef _RHEOLEF_HAVE_UMFPACK
107  if (prefer_umfpack) {
108  trace_macro ("use umfpack");
109  // diagonal bloc per process => use umfpack
110  typedef solver_umfpack_rep<T,M> rep;
111  if (_ptr) delete_macro (_ptr);
112  _ptr = new_macro (rep(a,opt));
113  return;
114  }
115 #endif // _RHEOLEF_HAVE_UMFPACK
116 
117 #ifdef _RHEOLEF_HAVE_MUMPS
118  trace_macro ("use mumps");
119  typedef solver_mumps_rep<T,M> rep;
120 #elif defined(_RHEOLEF_HAVE_PASTIX)
121  typedef solver_pastix_rep<T,M> rep;
122  opt.iterative = 0;
123 #else // ! _RHEOLEF_HAVE_TRILINOS && !_RHEOLEF_HAVE_PASTIX
125 #endif // ! _RHEOLEF_HAVE_MUMPS
126 
127  if (_ptr) delete_macro (_ptr);
128  _ptr = new_macro (rep(a,opt));
129 }
130 template<class T, class M>
131 void
133 {
134 #ifdef _RHEOLEF_HAVE_TRILINOS
135  typedef solver_trilinos_ifpack_rep<T,M> rep;
136 #elif defined(_RHEOLEF_HAVE_PASTIX)
137  typedef solver_pastix_rep<T,M> rep;
138  opt.iterative = 1;
139 #else // ! _RHEOLEF_HAVE_TRILINOS && !_RHEOLEF_HAVE_PASTIX
141 #endif
142  if (_ptr) delete_macro (_ptr);
143  _ptr = new_macro (rep(a,opt));
144 }
145 template<class T, class M>
147 {
148  if (_ptr) delete_macro (_ptr);
149  _ptr = 0;
150 }
151 template<class T, class M>
152 void
154 {
155  _ptr->update_values (a);
156 }
157 template<class T, class M>
158 vec<T,M>
160 {
161  return _ptr->solve (b);
162 }
163 template<class T, class M>
164 vec<T,M>
166 {
167  return _ptr->trans_solve (b);
168 }
169 #define _RHEOLEF_instanciation(T,M) \
170 template class solver_abstract_rep<T,M>; \
171 template class solver_rep<T,M>;
172 
173 _RHEOLEF_instanciation(Float,sequential)
174 #ifdef _RHEOLEF_HAVE_MPI
175 _RHEOLEF_instanciation(Float,distributed)
176 #endif // _RHEOLEF_HAVE_MPI
177 #undef _RHEOLEF_instanciation
178 
179 } // namespace rheolef
virtual determinant_type det() const
Definition: solver.cc:31
void build_lu(const csr< T, M > &a, const solver_option &opt)
Definition: solver.cc:74
#define new_macro(obj)
Definition: compiler.h:361
#define fatal_macro(message)
Definition: compiler.h:98
virtual void update_values(const csr< T, M > &a)=0
solver_abstract_rep< T, M > rep
Definition: solver.h:48
#define delete_macro(ptr)
Definition: compiler.h:363
vec< T, M > solve(const vec< T, M > &b) const
Definition: solver.cc:159
void update_values(const csr< T, M > &a)
Definition: solver.cc:153
eye - the identity matrix
Definition: eye.h:22
const solver_option & option() const
Definition: solver.cc:59
irheostream, orheostream - large data streams
Definition: compiler.h:7
void build_ilu(const csr< T, M > &a, const solver_option &opt)
Definition: solver.cc:132
#define warning_macro(message)
Definition: compiler.h:102
#define check_macro(ok_condition, message)
Definition: compiler.h:104
double Float
Definition: compiler.h:160
solver - direct or interative solver interface
Definition: solver.h:8
#define _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
Definition: config.h:233
void build_eye()
Definition: solver.cc:66
#define trace_macro(message)
Definition: compiler.h:120
determinant_type det() const
Definition: solver.cc:39
virtual vec< T, M > solve(const vec< T, M > &b) const =0
virtual vec< T, M > trans_solve(const vec< T, M > &b) const =0
rheolef::std type
vec - vector in distributed environment
Definition: vec.h:48
#define _RHEOLEF_instanciation(T, M)
Definition: solver.cc:169
vec< T, M > trans_solve(const vec< T, M > &b) const
Definition: solver.cc:165
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
std::string prefered_library