Rheolef  7.1
an efficient C++ finite element environment
solver_umfpack.cc
Go to the documentation of this file.
1 // direct solver UMFPACK, seq implementations
22 //
23 // Note : why a dis implementation based on umfpack ?
24 // Because when dis_ext_nnz == 0, then the matrix is block diagonal.
25 // in that case the umfpack could be better than mumps that initialize stuff
26 // for the distributed case.
27 // Is could appends e.g. for block-diagonal mass matrix "Pkd"
28 // This also occurs when nproc==1.
29 //
30 #include "rheolef/config.h"
31 #ifdef _RHEOLEF_HAVE_UMFPACK
32 #include "solver_umfpack.h"
33 #include "rheolef/rheostream.h" // itos
34 
35 namespace rheolef {
36 using namespace std;
37 
38 // =========================================================================
39 // umfpack utilities
40 // =========================================================================
41 static
42 std::string
43 umfpack_message (int status) {
44  switch (status) {
45  case UMFPACK_OK: return "ok";
46  case UMFPACK_WARNING_singular_matrix: return "singular matrix";
47  case UMFPACK_WARNING_determinant_underflow: return "determinant underflow";
48  case UMFPACK_WARNING_determinant_overflow: return "determinant overflow";
49  case UMFPACK_ERROR_out_of_memory: return "out of memory";
50  case UMFPACK_ERROR_invalid_Numeric_object: return "invalid Numeric object";
51  case UMFPACK_ERROR_invalid_Symbolic_object: return "invalid Symbolic object";
52  case UMFPACK_ERROR_argument_missing: return "argument missing";
53  case UMFPACK_ERROR_n_nonpositive: return "size is less or equal to zero";
54  case UMFPACK_ERROR_invalid_matrix: return "invalid matrix";
55  case UMFPACK_ERROR_different_pattern: return "different pattern";
56  case UMFPACK_ERROR_invalid_system: return "invalid system";
57  case UMFPACK_ERROR_invalid_permutation: return "invalid permutation";
58  case UMFPACK_ERROR_internal_error: return "internal error";
59  case UMFPACK_ERROR_file_IO : return "file i/o error";
60  default: return "unexpected umfpack error status = " + itos(status);
61  }
62 }
63 static
64 void
65 umfpack_check_error (int status) {
66  if (status == UMFPACK_OK) return;
67  warning_macro (umfpack_message(status));
68 }
69 // =========================================================================
70 // the class interface
71 // =========================================================================
72 template<class T, class M>
75  _n(0),
76  _numeric(0),
77  _ia_p(0),
78  _ja_p(0),
79  _a_p(0),
80  _det(),
81  _control(),
82  _info()
83 {
84 }
85 template<class T, class M>
87  : solver_abstract_rep<T,M>(opt),
88  _n(0),
89  _numeric(0),
90  _ia_p(0),
91  _ja_p(0),
92  _a_p(0),
93  _det(),
94  _control(),
95  _info()
96 {
97  _set_control();
98  update_values (a);
99 }
100 template<class T, class M>
103  _n(0),
104  _numeric(0),
105  _ia_p(0),
106  _ja_p(0),
107  _a_p(0),
108  _det(),
109  _control(),
110  _info()
111 {
112  // -Weff_c++ -Werror requires it, because has pointers, but should not happened
113  error_macro ("solver_umfpack_rep(const solver_umfpack_rep&) : should not happened");
114 }
115 template<class T, class M>
116 solver_umfpack_rep<T,M>&
117 solver_umfpack_rep<T,M>::operator= (const solver_umfpack_rep<T,M>&)
118 {
119  // -Weff_c++ -Werror requires it, because has pointers, but should not happened
120  error_macro ("solver_umfpack_rep::op=(const solver_umfpack_rep&) : should not happened");
121  return *this;
122 }
123 template<class T, class M>
124 void
126 {
127  umfpack_di_defaults (_control);
128  _control [UMFPACK_IRSTEP] = base::option().n_refinement;
129 }
130 template<class T, class M>
132 {
133  _destroy();
134 }
135 template<class T, class M>
136 void
138 {
139  if (_numeric) umfpack_di_free_numeric (&_numeric);
140  if (_ia_p) delete_tab_macro (_ia_p);
141  if (_ja_p) delete_tab_macro (_ja_p);
142  if (_a_p) delete_tab_macro (_a_p);
143  _n = 0;
144 }
145 template<class T, class M>
146 void
148 {
149  check_macro (base::option().force_seq || a.dis_ext_nnz() == 0, "unexpected non-zero dis_ext_nnz="<<a.dis_nnz());
150  _destroy(); // TODO: check if the sparse structure is unchanded: could be reused
151  if (a.nrow() == 0 || a.nnz() == 0) return; // empty matrix
152  _n = int(a.nrow());
153  _ia_p = new_tab_macro (int, _n+1);
154  _ja_p = new_tab_macro (int, a.nnz());
155  _a_p = new_tab_macro (T, a.nnz());
156  typename csr<T,M>::const_iterator ia = a.begin();
157  typedef typename csr<T,M>::size_type size_type;
158  _ia_p [0] = 0;
159  for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
160  _ia_p [i+1] = ia[i+1] - ia[0];
161  for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p, ++q) {
162  _ja_p [q] = (*p).first;
163  _a_p [q] = (*p).second;
164  }
165  }
166  void *symbolic;
167  umfpack_di_symbolic (_n, _n, _ia_p, _ja_p, _a_p, &symbolic, _control, _info);
168  umfpack_check_error (int(_info [UMFPACK_STATUS]));
169  umfpack_di_numeric ( _ia_p, _ja_p, _a_p, symbolic, &_numeric, _control, _info);
170  umfpack_check_error (int(_info [UMFPACK_STATUS]));
171  umfpack_di_free_symbolic (&symbolic);
172  if (base::option().compute_determinant) {
173  _det.base = 10;
174  int status = umfpack_di_get_determinant (&_det.mantissa, &_det.exponant, _numeric, _info);
175  if (status != 0) umfpack_check_error (int(_info [UMFPACK_STATUS]));
176  }
177 }
178 template<class T, class M>
179 void
180 solver_umfpack_rep<T,M>::_solve (int transpose_flag, const vec<T,M>& b, vec<T,M>& x) const
181 {
182  umfpack_di_solve (transpose_flag, _ia_p, _ja_p, _a_p, x.begin().operator->(), b.begin().operator->(), _numeric, _control, _info);
183  umfpack_check_error (int(_info [UMFPACK_STATUS]));
184 }
185 template<class T, class M>
186 vec<T,M>
188 {
189  if (_n == 0) return b; // empty matrix
190  vec<T,M> x(b.ownership());
191  _solve (UMFPACK_At, b, x); // umfpack uses csc while rheolef uses csr
192  return x;
193 }
194 template<class T, class M>
195 vec<T,M>
197 {
198  if (_n == 0) return b; // empty matrix
199  vec<T,M> x(b.ownership());
200  _solve (UMFPACK_A, b, x); // umfpack uses csc while rheolef uses csr
201  return x;
202 }
203 // ----------------------------------------------------------------------------
204 // instanciation in library
205 // ----------------------------------------------------------------------------
206 // TODO: code is only valid here for T=double
207 
209 
210 #ifdef _RHEOLEF_HAVE_MPI
212 #endif // _RHEOLEF_HAVE_MPI
213 
214 } // namespace rheolef
215 #endif // HAVE_UMFPACK
solver_umfpack.h
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::solver_umfpack_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition: solver_umfpack.cc:187
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
p
Definition: sphere.icc:25
rheolef::solver_umfpack_rep::_destroy
void _destroy()
Definition: solver_umfpack.cc:137
rheolef::solver_abstract_rep
Definition: solver.h:191
a
Definition: diffusion_isotropic.h:25
rheolef::csr
see the csr page for the full documentation
Definition: asr.h:31
rheolef::solver_umfpack_rep::_set_control
void _set_control()
Definition: solver_umfpack.cc:125
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::solver_abstract_rep::size_type
csr< T, M >::size_type size_type
Definition: solver.h:193
rheolef::solver_umfpack_rep::~solver_umfpack_rep
~solver_umfpack_rep()
Definition: solver_umfpack.cc:131
rheolef::solver_umfpack_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition: solver_umfpack.cc:196
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::solver_umfpack_rep
Definition: solver_umfpack.h:40
rheolef::solver_umfpack_rep::_solve
void _solve(int transpose_flag, const vec< T, M > &b, vec< T, M > &x) const
Definition: solver_umfpack.cc:180
rheolef::solver_umfpack_rep::solver_umfpack_rep
solver_umfpack_rep()
Definition: solver_umfpack.cc:73
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::solver_umfpack_rep::update_values
void update_values(const csr< T, M > &a)
Definition: solver_umfpack.cc:147
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
M
Expr1::memory_type M
Definition: vec_expr_v2.h:385
T
Expr1::float_type T
Definition: field_expr.h:218
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155