Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_umfpack.cc
Go to the documentation of this file.
1
21// 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
34namespace rheolef {
35using namespace std;
36
37// =========================================================================
38// umfpack utilities
39// =========================================================================
40static
41std::string
42umfpack_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}
62static
63void
64umfpack_check_error (int status) {
65 if (status == UMFPACK_OK) return;
66 warning_macro (umfpack_message(status));
67}
68// =========================================================================
69// the class interface
70// =========================================================================
71template<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}
84template<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{
97 update_values (a);
98}
99template<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}
114template<class T, class M>
115solver_umfpack_rep<T,M>&
116solver_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}
122template<class T, class M>
123void
125{
126 umfpack_di_defaults (_control);
127 _control [UMFPACK_IRSTEP] = base::option().n_refinement;
128}
129template<class T, class M>
134template<class T, class M>
135void
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}
144template<class T, class M>
145void
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}
177template<class T, class M>
178void
179solver_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}
184template<class T, class 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}
193template<class T, class 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
field::size_type size_type
Definition branch.cc:430
see the csr page for the full documentation
Definition csr.h:317
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
see the vec page for the full documentation
Definition vec.h:79
#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.
STL namespace.
Definition sphere.icc:25
Expr1::memory_type M