30#include "rheolef/config.h"
31#ifdef _RHEOLEF_HAVE_UMFPACK
42umfpack_message (
int 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);
64umfpack_check_error (
int status) {
65 if (status == UMFPACK_OK)
return;
71template<
class T,
class M>
84template<
class T,
class M>
99template<
class T,
class M>
112 error_macro (
"solver_umfpack_rep(const solver_umfpack_rep&) : should not happened");
114template<
class T,
class M>
115solver_umfpack_rep<T,M>&
116solver_umfpack_rep<T,M>::operator= (
const solver_umfpack_rep<T,M>&)
119 error_macro (
"solver_umfpack_rep::op=(const solver_umfpack_rep&) : should not happened");
122template<
class T,
class M>
126 umfpack_di_defaults (_control);
127 _control [UMFPACK_IRSTEP] = base::option().n_refinement;
129template<
class T,
class M>
134template<
class T,
class M>
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);
144template<
class T,
class M>
148 check_macro (base::option().force_seq || a.dis_ext_nnz() == 0,
"unexpected non-zero dis_ext_nnz="<<a.dis_nnz());
150 if (a.nrow() == 0 || a.nnz() == 0)
return;
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());
158 for (
size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
159 _ia_p [i+1] = ia[i+1] - ia[0];
161 _ja_p [q] = (*p).first;
162 _a_p [q] = (*p).second;
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) {
173 int status = umfpack_di_get_determinant (&_det.mantissa, &_det.exponant, _numeric, _info);
174 if (status != 0) umfpack_check_error (
int(_info [UMFPACK_STATUS]));
177template<
class T,
class M>
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]));
184template<
class T,
class M>
188 if (_n == 0)
return b;
190 _solve (UMFPACK_At, b, x);
193template<
class T,
class M>
197 if (_n == 0)
return b;
199 _solve (UMFPACK_A, b, x);
209#ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the csr page for the full documentation
see the solver_option page for the full documentation
void update_values(const csr< T, M > &a)
base::size_type size_type
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
#define error_macro(message)
#define warning_macro(message)
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.