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"
 
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 = " + std::to_string(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