solver_umfpack.cc
Go to the documentation of this file.
// direct solver UMFPACK, seq implementations
//
// Note : why a dis implementation based on umfpack ?
// Because when dis_ext_nnz == 0, then the matrix is block diagonal.
// in that case the umfpack could be better than mumps that initialize stuff
// for the distributed case.
// Is could appends e.g. for block-diagonal mass matrix "Pkd"
// This also occurs when nproc==1.
//
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_UMFPACK
#include "solver_umfpack.h"
#include "rheolef/rheostream.h" // itos
namespace rheolef {
using namespace std;
// =========================================================================
// umfpack utilities
// =========================================================================
static
std::string
umfpack_message (int status) {
switch (status) {
case UMFPACK_OK: return "ok";
case UMFPACK_WARNING_singular_matrix: return "singular matrix";
case UMFPACK_WARNING_determinant_underflow: return "determinant underflow";
case UMFPACK_WARNING_determinant_overflow: return "determinant overflow";
case UMFPACK_ERROR_out_of_memory: return "out of memory";
case UMFPACK_ERROR_invalid_Numeric_object: return "invalid Numeric object";
case UMFPACK_ERROR_invalid_Symbolic_object: return "invalid Symbolic object";
case UMFPACK_ERROR_argument_missing: return "argument missing";
case UMFPACK_ERROR_n_nonpositive: return "size is less or equal to zero";
case UMFPACK_ERROR_invalid_matrix: return "invalid matrix";
case UMFPACK_ERROR_different_pattern: return "different pattern";
case UMFPACK_ERROR_invalid_system: return "invalid system";
case UMFPACK_ERROR_invalid_permutation: return "invalid permutation";
case UMFPACK_ERROR_internal_error: return "internal error";
case UMFPACK_ERROR_file_IO : return "file i/o error";
default: return "unexpected umfpack error status = " + itos(status);
}
}
static
void
umfpack_check_error (int status) {
if (status == UMFPACK_OK) return;
warning_macro (umfpack_message(status));
}
// =========================================================================
// the class interface
// =========================================================================
template<class T, class M>
solver_umfpack_rep<T,M>::solver_umfpack_rep ()
: solver_abstract_rep<T,M>(solver_option()),
_n(0),
_numeric(0),
_ia_p(0),
_ja_p(0),
_a_p(0),
_det(),
_control(),
_info()
{
}
template<class T, class M>
solver_umfpack_rep<T,M>::solver_umfpack_rep (const csr<T,M>& a, const solver_option& opt)
: solver_abstract_rep<T,M>(opt),
_n(0),
_numeric(0),
_ia_p(0),
_ja_p(0),
_a_p(0),
_det(),
_control(),
_info()
{
_set_control();
update_values (a);
}
template<class T, class M>
solver_umfpack_rep<T,M>::solver_umfpack_rep (const solver_umfpack_rep<T,M>&)
: solver_abstract_rep<T,M>(solver_option()),
_n(0),
_numeric(0),
_ia_p(0),
_ja_p(0),
_a_p(0),
_det(),
_control(),
_info()
{
// -Weff_c++ -Werror requires it, because has pointers, but should not happened
error_macro ("solver_umfpack_rep(const solver_umfpack_rep&) : should not happened");
}
template<class T, class M>
solver_umfpack_rep<T,M>&
solver_umfpack_rep<T,M>::operator= (const solver_umfpack_rep<T,M>&)
{
// -Weff_c++ -Werror requires it, because has pointers, but should not happened
error_macro ("solver_umfpack_rep::op=(const solver_umfpack_rep&) : should not happened");
return *this;
}
template<class T, class M>
void
solver_umfpack_rep<T,M>::_set_control()
{
umfpack_di_defaults (_control);
_control [UMFPACK_IRSTEP] = base::option().n_refinement;
}
template<class T, class M>
solver_umfpack_rep<T,M>::~solver_umfpack_rep ()
{
_destroy();
}
template<class T, class M>
void
solver_umfpack_rep<T,M>::_destroy ()
{
if (_numeric) umfpack_di_free_numeric (&_numeric);
if (_ia_p) delete_tab_macro (_ia_p);
if (_ja_p) delete_tab_macro (_ja_p);
if (_a_p) delete_tab_macro (_a_p);
_n = 0;
}
template<class T, class M>
void
solver_umfpack_rep<T,M>::update_values (const csr<T,M>& a)
{
check_macro (base::option().force_seq || a.dis_ext_nnz() == 0, "unexpected non-zero dis_ext_nnz="<<a.dis_nnz());
_destroy(); // TODO: check if the sparse structure is unchanded: could be reused
if (a.nrow() == 0 || a.nnz() == 0) return; // empty matrix
_n = int(a.nrow());
_ia_p = new_tab_macro (int, _n+1);
_ja_p = new_tab_macro (int, a.nnz());
_a_p = new_tab_macro (T, a.nnz());
typename csr<T,M>::const_iterator ia = a.begin();
typedef typename csr<T,M>::size_type size_type;
_ia_p [0] = 0;
for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
_ia_p [i+1] = ia[i+1] - ia[0];
for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p, ++q) {
_ja_p [q] = (*p).first;
_a_p [q] = (*p).second;
}
}
void *symbolic;
umfpack_di_symbolic (_n, _n, _ia_p, _ja_p, _a_p, &symbolic, _control, _info);
umfpack_check_error (int(_info [UMFPACK_STATUS]));
umfpack_di_numeric ( _ia_p, _ja_p, _a_p, symbolic, &_numeric, _control, _info);
umfpack_check_error (int(_info [UMFPACK_STATUS]));
umfpack_di_free_symbolic (&symbolic);
if (base::option().compute_determinant) {
_det.base = 10;
int status = umfpack_di_get_determinant (&_det.mantissa, &_det.exponant, _numeric, _info);
if (status != 0) umfpack_check_error (int(_info [UMFPACK_STATUS]));
}
}
template<class T, class M>
void
solver_umfpack_rep<T,M>::_solve (int transpose_flag, const vec<T,M>& b, vec<T,M>& x) const
{
umfpack_di_solve (transpose_flag, _ia_p, _ja_p, _a_p, x.begin().operator->(), b.begin().operator->(), _numeric, _control, _info);
umfpack_check_error (int(_info [UMFPACK_STATUS]));
}
template<class T, class M>
vec<T,M>
solver_umfpack_rep<T,M>::solve (const vec<T,M>& b) const
{
if (_n == 0) return b; // empty matrix
vec<T,M> x(b.ownership());
_solve (UMFPACK_At, b, x); // umfpack uses csc while rheolef uses csr
return x;
}
template<class T, class M>
vec<T,M>
solver_umfpack_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
if (_n == 0) return b; // empty matrix
vec<T,M> x(b.ownership());
_solve (UMFPACK_A, b, x); // umfpack uses csc while rheolef uses csr
return x;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
// TODO: code is only valid here for T=double
template class solver_umfpack_rep<double,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_umfpack_rep<double,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // HAVE_UMFPACK
rheolef::csr< T, M >
rheolef::solver_abstract_rep
Definition: solver.h:191
rheolef::solver_abstract_rep::size_type
csr< T, M >::size_type size_type
Definition: solver.h:193
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155
rheolef::solver_umfpack_rep
Definition: solver_umfpack.h:40
rheolef::solver_umfpack_rep::_destroy
void _destroy()
Definition: solver_umfpack.cc:137
rheolef::solver_umfpack_rep::_set_control
void _set_control()
Definition: solver_umfpack.cc:125
rheolef::solver_umfpack_rep::update_values
void update_values(const csr< T, M > &a)
Definition: solver_umfpack.cc:147
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::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition: solver_umfpack.cc:196
rheolef::solver_umfpack_rep::~solver_umfpack_rep
~solver_umfpack_rep()
Definition: solver_umfpack.cc:131
rheolef::solver_umfpack_rep::solver_umfpack_rep
solver_umfpack_rep()
Definition: solver_umfpack.cc:73
rheolef::solver_umfpack_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition: solver_umfpack.cc:187
rheolef::vec< T, M >
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
solver_umfpack.h
a
Definition: diffusion_isotropic.h:25
p
Definition: sphere.icc:25
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416