solver_mumps.cc

Go to the documentation of this file.

// direct solver MUMPS, mpi implementations
//
// Note : why no seq implementation based on mumps ?
//
// 1. mumps has a pure seq library implementation.
//   Both seq and mpi libraries cannot be linked together because they define
//   both the same dmumps() function.
// 2. but, when using mpi, the mpi implementation is unable to solve efficiently
//   in // some independant sparse linear systems (eg. sparse mass matrix local
//   on an element) because the right-hand side may be merged on the proc=0 
//   for an obscure raison.
//
// Thus rheolef uses mumps for mpi-only implementation of the direct solver
// and use another library (eg umfpack) for the seq direct solver.
//
// The present 'seq' implementation has been tested and is invalid :
// the code is maintained (but no more instanciated) for two reasons:
// 1. it can be used when mpi is not present, as a seq solver with the pure seq impl of mumps
// 2. mumps will propose in the future to provide a distributed RHS:
//      https://listes.ens-lyon.fr/sympa/arc/mumps-users
//      From: Alfredo Buttari <alfredo.buttari@enseeiht.fr>
//      Date: Mon, 31 Jan 2011 19:07:03 +0100
//      Subject: Re: [mumps-users] Scalability of MUMPS solution phase
//      Jack,
//      yes, adding support for distributed RHS is on our TODO list but I
//      can't say at this moment when it will be available.
//
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_MUMPS
#include "solver_mumps.h"
 
#undef  DEBUG_MUMPS_SCOTCH
#ifdef  DEBUG_MUMPS_SCOTCH
#undef _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH // has failed on some linear systems: prefer seq-scotch when available
#endif // DEBUG_MUMPS
 
namespace rheolef {
using namespace std;
 
// =========================================================================
// 1. local functions
// =========================================================================
// count non-zero entries of the upper part of A
// 1) diagonal block case
template<class T, class M>
static
typename csr<T,M>::size_type
nnz_dia_upper (const csr<T,M>& a)
{
  typedef typename csr<T,M>::size_type size_type;
  size_type dia_nnz = 0;
  typename csr<T,M>::const_iterator dia_ia = a.begin();
  for (size_type i = 0, n = a.nrow(); i < n; i++) {
    for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
      size_type j = (*p).first;
      if (j >= i) dia_nnz++;
    }
  }
  return dia_nnz;
}
// 2) extra-diagonal block case, for distributed csr
template<class T, class M>
static
typename csr<T,M>::size_type
nnz_ext_upper (const csr<T,M>& a)
{
  return 0;
}
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
static
typename csr<T,distributed>::size_type
nnz_ext_upper (const csr<T,distributed>& a)
{
  typedef typename csr<T,distributed>::size_type size_type;
  size_type first_i = a.row_ownership().first_index();
  size_type ext_nnz = 0;
  typename csr<T,distributed>::const_iterator ext_ia = a.ext_begin();
  for (size_type i = 0, n = a.nrow(); i < n; i++) {
    for (typename csr<T,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
      size_type dis_i = first_i + i;
      size_type dis_j = a.jext2dis_j ((*p).first);
      if (dis_j >= dis_i) ext_nnz++;
    }
  }
  return ext_nnz;
}
#endif // _RHEOLEF_HAVE_MPI
template<class T, class M>
static
void
csr2mumps_struct_seq (const csr<T,M>& a, vector<int>& row, vector<int>& col, bool use_symmetry)
{
  typedef typename csr<T,M>::size_type size_type;
  typename csr<T,M>::const_iterator dia_ia = a.begin();
  for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
    for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
      size_type j = (*p).first; 
      if (!use_symmetry || j >= i) {
        row[q] = i + 1;
        col[q] = j + 1;
        q++;
      }
    }
  }
}
template<class T, class M>
static
void
csr2mumps_values_seq (const csr<T,M>& a, vector<T>& val, bool use_symmetry)
{
  typedef typename csr<T,M>::size_type size_type;
  typename csr<T,M>::const_iterator dia_ia = a.begin();
  for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
    for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
      size_type j = (*p).first; 
      if (!use_symmetry || j >= i) {
        val[q] = (*p).second;
        q++;
      }
    }
  }
}
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
static
void
csr2mumps_struct_dis (const csr<T,sequential>& a, vector<int>& row, vector<int>& col, bool use_symmetry, bool)
{
  // dummy case: for the compiler to be happy
}
template<class T>
static
void
csr2mumps_struct_dis (const csr<T,distributed>& a, vector<int>& row, vector<int>& col, bool use_symmetry, bool drop_ext_nnz)
{
  typedef typename csr<T,distributed>::size_type size_type;
  typename csr<T,distributed>::const_iterator dia_ia = a.begin();
  typename csr<T,distributed>::const_iterator ext_ia = a.ext_begin();
  size_type dis_nr = a.dis_nrow();
  size_type dis_nc = a.dis_ncol();
  size_type nc = a.ncol();
  size_type first_i = a.row_ownership().first_index();
  for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
    for (typename csr<T,distributed>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
      size_type j = (*p).first; 
      size_type dis_i = first_i + i;
      size_type dis_j = first_i + j;    
      assert_macro (j < nc, "invalid matrix");
      assert_macro (dis_i < dis_nr && dis_j < dis_nc, "invalid matrix");
      if (!use_symmetry || dis_j >= dis_i) {
        row[q] = dis_i + 1;
        col[q] = dis_j + 1;
        q++;
      }
    }
    if (drop_ext_nnz) continue;
    for (typename csr<T,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
      size_type dis_i = first_i + i;
      size_type dis_j = a.jext2dis_j ((*p).first);
      assert_macro (dis_i < dis_nr && dis_j < dis_nc, "invalid matrix");
      if (!use_symmetry || dis_j >= dis_i) {
        row[q] = dis_i + 1;
        col[q] = dis_j + 1;
        q++;
      }
    }
  }
}
template<class T>
static
void
csr2mumps_values_dis (const csr<T,sequential>& a, vector<T>& val, bool use_symmetry, bool)
{
  // dummy case: for the compiler to be happy
}
template<class T>
static
void
csr2mumps_values_dis (const csr<T,distributed>& a, vector<T>& val, bool use_symmetry, bool drop_ext_nnz)
{
#ifdef TO_CLEAN
  std:fill (val.begin(), val.end(), std::numeric_limits<T>::max());
#endif // TO_CLEAN
 
  typedef typename csr<T,distributed>::size_type size_type;
  size_type first_i = a.row_ownership().first_index();
  typename csr<T,distributed>::const_iterator dia_ia = a.begin();
  typename csr<T,distributed>::const_iterator ext_ia = a.ext_begin();
  for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
    for (typename csr<T,distributed>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
      size_type dis_i = first_i + i;
      size_type dis_j = first_i + (*p).first;   
      if (!use_symmetry || dis_j >= dis_i) {
#ifdef TO_CLEAN
        check_macro(q < val.size(), "invalid index");
#endif // TO_CLEAN
        val[q] = (*p).second;
        q++;
      }
    }
    if (drop_ext_nnz) continue;
    for (typename csr<T,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
      size_type dis_i = first_i + i;
      size_type dis_j = a.jext2dis_j ((*p).first);
      if (!use_symmetry || dis_j >= dis_i) {
#ifdef TO_CLEAN
        check_macro(q < val.size(), "invalid index");
#endif // TO_CLEAN
        val[q] = (*p).second;
        q++;
      }
    }
  }
#ifdef TO_CLEAN
  T val_min = std::numeric_limits<T>::max(), val_max = 0;
  for (size_t q = 0; q < val.size(); q++) {
    val_min = std::min (val_min, abs(val[q]));
    val_max = std::max (val_max, abs(val[q]));
  }
  warning_macro ("val_min="<<val_min<<", val_max="<<val_max);
#endif // TO_CLEAN
}
#endif // _RHEOLEF_HAVE_MPI
// =========================================================================
// 2. the class interface
// =========================================================================
template<class T, class M>
solver_mumps_rep<T,M>::solver_mumps_rep (const csr<T,M>& a, const solver_option& opt)
 : solver_abstract_rep<T,M>(opt),
   _has_mumps_instance(false),
   _drop_ext_nnz(false),
   _mumps_par(),
   _row(),
   _col(),
   _val(),
   _a00(0),
   _det()
{
  _drop_ext_nnz = is_distributed<M>::value && opt.force_seq;
  update_values (a);
}
template<class T, class M>
void
solver_mumps_rep<T,M>::update_values (const csr<T,M>& a)
{
  size_type nproc = communicator().size();
  if (a.dis_nrow() <= 1) {
    // skip 0 or 1 dis_nrow, where mumps core dump
    if (a.nrow() == 1) {
      _a00 = (*(a.begin()[0])).second;
    }
    _has_mumps_instance = false;
    return;
  }
  _has_mumps_instance = true;
  // ----------------------------------
  //  step 0 : init a mumps instance
  // ----------------------------------
  bool use_symmetry          = a.is_symmetric();
  bool use_definite_positive = a.is_definite_positive();
  _mumps_par.par          = 1; // use parallel
#ifdef _RHEOLEF_HAVE_MPI
  _mumps_par.comm_fortran = MPI_Comm_c2f (a.row_ownership().comm());
#endif // _RHEOLEF_HAVE_MPI
  if (use_symmetry && !use_definite_positive) {
    // sym invertible but could be undef: could have either < 0 or > 0 eigenvalues
    // => leads to similar cpu with s.d.p. but without the asumption of positive or negative :
    _mumps_par.sym        = 2;
  } else if (use_symmetry && use_definite_positive) {
#if 0
    // sym def pos
    // a) standard approach
    //    Note: could fail in scalapack when sym def negative
    _mumps_par.sym        = 1;
    // b) second approach: from AmeButGue-2011-mumps, page 15:
    // Another approach to suppress numerical pivoting which works with ScaLAPACK
    // for both positive definite and negative definite matrices consists in setting SYM=2 and
    // CNTL(1)=0.0D0 (recommended strategy).
    _mumps_par.sym        = 2;
    _mumps_par.cntl[1-1]  = 0.0;
#endif // 0
    // c) third approach: do not use the knowledge of positive or negative, since
    //    this knowledge is contained in bilinear forms defined by the user
    //    use only the general symetry:
    //    tests show CPU is very similar than two previous approaches
    //    and the user do not need know about positive_definite or negative_definite
    _mumps_par.sym        = 2;
  } else {
    // unsym invertible
    _mumps_par.sym        = 0;
  }
  _mumps_par.job = -1;
  dmumps_c(&_mumps_par);
  check_macro (_mumps_par.infog[1-1] == 0, "mumps: an error has occurred");
  // ----------------------------------
  //  step 1 : symbolic factorization
  // ----------------------------------
  if (base::option().verbose_level != 0) {
    _mumps_par.icntl[1-1] = 1; // error output
    _mumps_par.icntl[2-1] = 1; // verbose output
    _mumps_par.icntl[3-1] = 1; // global infos
    _mumps_par.icntl[4-1] = base::option().verbose_level; // verbosity level
    // strcpy (_mumps_par.write_problem, "dump_mumps"); // dump sparse structure
  } else {
    _mumps_par.icntl[1-1] = -1;
    _mumps_par.icntl[2-1] = -1;
    _mumps_par.icntl[3-1] = -1;
    _mumps_par.icntl[4-1] =  0; 
  }
  if (base::option().compute_determinant) {
    _mumps_par.icntl[33-1] = 1; 
  }
  _mumps_par.icntl[5-1] = 0; // matrix is in assembled format
#if defined(_RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_SCOTCH)
  _mumps_par.icntl[7-1] = 3; // sequential ordering: 3==scotch
#elif defined(_RHEOLEF_HAVE_MUMPS_WITH_PARMETIS) || defined(_RHEOLEF_HAVE_MUMPS_WITH_METIS)
  _mumps_par.icntl[7-1] = 5; // sequential ordering: 5==metis
#else // _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS
  _mumps_par.icntl[7-1] = 7; // ordering: 7=auto
#endif // _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
  _mumps_par.icntl[14-1] = 40; // default 20 seems to be insufficient in some cases
  _mumps_par.icntl[29-1] = 0; // 0: auto; 1: ptscotch; 2: parmetis
  _mumps_par.icntl[22-1] = 0; // 0: in-core ; 1: out-of-core TODO
  if (is_distributed<M>::value) {
    _mumps_par.icntl[18-1] = 3; // distributed
 
#if defined(_RHEOLEF_HAVE_MUMPS_WITH_SCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_METIS)
    _mumps_par.icntl[28-1] = 1; // sequential ordering (preferred, more robust)
#elif defined(_RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_PARMETIS)
    _mumps_par.icntl[28-1] = 2; // parallel ordering (less robust, but faster)
#endif // _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH || _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS
 
#ifdef _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
    _mumps_par.icntl[29-1] = 1; // ptscotch=1
#elif _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS
    _mumps_par.icntl[29-1] = 2; // parmetis=2
#else
    _mumps_par.icntl[29-1] = 0; // automatic choice
#endif // _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
  } else {
    _mumps_par.icntl[18-1] = 0; // sequential
  }
  // load the matrix in mumps
  _mumps_par.n = a.dis_nrow();
  size_type dia_nnz = ((use_symmetry) ? nnz_dia_upper(a) : a.nnz());
  size_type ext_nnz = ((use_symmetry) ? nnz_ext_upper(a) : a.ext_nnz());
  size_type dis_nnz = dia_nnz + (_drop_ext_nnz ? 0 : ext_nnz);
  if (is_distributed<M>::value) {
#ifdef _RHEOLEF_HAVE_MPI
    _row.resize (dis_nnz);
    _col.resize (dis_nnz);
    csr2mumps_struct_dis (a, _row, _col, use_symmetry, _drop_ext_nnz);
    _mumps_par.nz_loc  = dis_nnz;
    _mumps_par.irn_loc = _row.begin().operator->();
    _mumps_par.jcn_loc = _col.begin().operator->();
#endif // _RHEOLEF_HAVE_MPI
  } else { // sequential
    _row.resize (dia_nnz);
    _col.resize (dia_nnz);
    csr2mumps_struct_seq (a, _row, _col, use_symmetry);
    _mumps_par.nz      = dia_nnz;
#if (_RHEOLEF_MUMPS_VERSION_MAJOR >= 5)
    _mumps_par.nnz     = dia_nnz;
#endif
    _mumps_par.irn     = _row.begin().operator->();
    _mumps_par.jcn     = _col.begin().operator->();
  }
  _mumps_par.job = 1;
  dmumps_c(&_mumps_par);
  trace_macro ("symbolic: ordering effectively used = " << _mumps_par.infog[7-1]); // 3=scoth, etc
  check_macro (_mumps_par.infog[1-1] == 0, "mumps: error infog(1)="<<_mumps_par.infog[1-1]<<" has occurred (see mumps/infog(1) documentation; infog(2)="<<_mumps_par.infog[2-1]<<")");
  if (a.dis_nrow() <= 1) {
    // skip 0 or 1 dis_nrow, where mumps core dump
    if (a.nrow() == 1) {
      _a00 = (*(a.begin()[0])).second;
    }
    return;
  }
  // ----------------------------------
  //  step 2 : numeric factorization
  // ----------------------------------
  if (is_distributed<M>::value) {
#ifdef _RHEOLEF_HAVE_MPI
    _val.resize (dis_nnz);
    csr2mumps_values_dis (a, _val, use_symmetry, _drop_ext_nnz);
    _mumps_par.a_loc   = _val.begin().operator->();
#endif // _RHEOLEF_HAVE_MPI
  } else { // sequential
    _val.resize (dia_nnz);
    csr2mumps_values_seq (a, _val, use_symmetry);
    _mumps_par.a       = _val.begin().operator->();
  }
  _mumps_par.job = 2;
  bool finished = false;
  while (!finished && _mumps_par.icntl[14-1] <= 2000) {
    dmumps_c(&_mumps_par);
    if ((_mumps_par.infog[1-1] == -8 || _mumps_par.infog[1-1] == -9) && 
         _mumps_par.infog[2-1] >= 0) {
      // not enought working space:
      // default 20 seems to be insufficient in some cases (periodic,surfacic,p-laplacian)
      _mumps_par.icntl[14-1] += 10;
      if (_mumps_par.icntl[14-1] > 100) {
        dis_warning_macro ("numeric factorization: increase working space of "<<_mumps_par.icntl[14-1]<< "% and retry...");
      }
    } else {
      finished = true;
    }
  }
  check_macro (_mumps_par.infog[1-1] == 0,
        "solver failed: mumps error infog(1)=" <<_mumps_par.infog[1-1]
                              << ", infog(2)=" <<_mumps_par.infog[2-1]
                              << ", icntl(14)="<<_mumps_par.icntl[14-1]
                              << ", icntl(23)="<<_mumps_par.icntl[23-1]);
  if (_mumps_par.icntl[33-1] != 0) {
    // get determinant from infos:
    _det.mantissa = _mumps_par.rinfog[12-1];
    //    imag part _mumps_par.rinfog[13-1] is zero as matrix is real here
    _det.exponant = _mumps_par.infog[34-1];
    _det.base     = 2;
  }
}
template<class T, class M>
vec<T,M>
solver_mumps_rep<T,M>::solve (const vec<T,M>& b) const
{
  if (b.dis_size() <= 1) {
    // skip 0 or 1 dis_nrow, where mumps core dump
    if (b.size() == 1) {
      return vec<T,M>(b.ownership(), b[0]/_a00);
    } else {
      return vec<T,M>(b.ownership());
    }
  }
  // ----------------------------------
  //  step 3 : define rhs & solve
  // ----------------------------------
  vec<T,M> x;
  vector<T>   sol;
  vector<int> perm;
  size_type sol_size = 0;
  vec<T,M> b_host;
  T      dummy; // mumps faild when zero sizes and 0 pointers...
  int    idummy;
  if (is_distributed<M>::value) {
    // 3.1a. merge the rhs b on proc=0 as required by distributed mumps (why ?)
    _mumps_par.icntl[21-1] = 1; // 1: solution is distributed
    _mumps_par.icntl[10-1] = 0; // nstep of iterative refinement (not in distributed version?)
    communicator comm = b.ownership().comm();
    distributor host_ownership (b.dis_size(), comm, comm.rank() == 0 ? b.dis_size() : 0);
    b_host.resize (host_ownership);
    size_type first_i = b.ownership().first_index();
    for (size_type i = 0, n = b.size(); i < n; i++) {
      size_type dis_i = first_i + i;
      b_host.dis_entry(dis_i) = b[i];
    }
    b_host.dis_entry_assembly();
    if (comm.rank() == 0) {
      _mumps_par.nrhs = 1;
      _mumps_par.rhs = b_host.begin().operator->();
    }
    sol_size = _mumps_par.info[23-1];
    sol.resize  (sol_size);
    perm.resize (sol_size);
    _mumps_par.lsol_loc = sol_size;
    _mumps_par.sol_loc  = (sol_size > 0) ?  sol.begin().operator->() : &dummy;
    _mumps_par.isol_loc = (sol_size > 0) ? perm.begin().operator->() : &idummy;
  } else { // sequential
    // 3.1b. inplace resolution: copy b in x and put x as rhs
    _mumps_par.icntl[21-1] = 0; // 0: sol goes on the proc=0, inplace in rhs
    _mumps_par.icntl[10-1] = 0; // nstep of iterative refinement (not in distributed version?)
    x = b;
    _mumps_par.nrhs = 1;
    _mumps_par.rhs = x.begin().operator->();
  }
  // 3.2. run mumps solver
  _mumps_par.icntl [9-1] = 1; // solve A*x=b : ctrl=1 ; otherwise solve A^t*x=b TODO
  _mumps_par.job = 3;
  dmumps_c(&_mumps_par);
  check_macro (_mumps_par.infog[1-1] >= 0,
        "solver failed: mumps error infog(1)="<<_mumps_par.infog[1-1]
                              << ", infog(2)="<<_mumps_par.infog[2-1]);
  if (_mumps_par.infog[1-1] > 0) {
        warning_macro ("mumps warning infog(1)="<<_mumps_par.infog[1-1]
                                << ", infog(2)="<<_mumps_par.infog[2-1]);
  }
  if (is_distributed<M>::value) {
    // 3.3. redistribute the solution as b.ownership
    x.resize (b.ownership());
    for (size_t i = 0; i < sol_size; i++) {
      size_type dis_i = perm[i] - 1;
      assert_macro (dis_i < x.dis_size(), "invalid index");
      x.dis_entry(dis_i) = sol[i];
    }
    x.dis_entry_assembly();
  }
  return x;
}
template<class T, class M>
vec<T,M>
solver_mumps_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
  error_macro ("not yet");
  return vec<T,M>();
}
template<class T, class M>
solver_mumps_rep<T,M>::~solver_mumps_rep ()
{
  if (!_has_mumps_instance) return;
  _has_mumps_instance = false;
  // ----------------------------------
  //  step F : close the mumps instance
  // ----------------------------------
  _mumps_par.job = -2;
  dmumps_c(&_mumps_par);
  check_macro (_mumps_par.infog[1-1] == 0, "mumps: an error has occurred");
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
// TODO: code is only valid here for T=double
 
template class solver_mumps_rep<double,sequential>;
 
#ifdef _RHEOLEF_HAVE_MPI
template class solver_mumps_rep<double,distributed>;
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef
#endif // MUMPS