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