| [eigen] Sample code for MUMPS backend |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
Hi,
Benoit replied to my issue tracker entry I should send a piece of code
to this list, what I do now.
I want to provide you with some example code for the MUMPS sparse
solver backend. I am not familiar with the development of Eigen, but I
hope I can assist one one you such that the adjustments will be quite
easy. The code is excerpted from an open source project which will be
released under BSD later this year. So there will be no problem if you
redistribute it under LGPLv3 or something else.
The idea of a MUMPS backend is the following:
- MUMPS supports single and double precision, real and complex numbers.
For each, another interface routine is provided for reverse
communication. These routines can be encapsulated in a unified template
class named "TMumpsInterface".
- In Eigen 2.0.x, the solvers work in two stages: 1st factorize and
store the factorized form of the input matrix. The input matrix can be
descarded afterwards. 2nd using the factorized matrix object, one can
solve dense vectors or matrices (either SolveInPlace or Solve). MUMPS
works in SolveInPlace mode, btw.
- To me it was also important that the same solver object can be used
several times, that means multiple times for solving equation systems
with equal system matrix and different rhs; and multiple times for
solving equation systems with different system matrix. Therefore, the
factorized (MUMPS internal) structure must be cleaned before
factorizing a new matrix (NOTE: This is not successfully implemented in
the other solver backends at the moment!)
- The sample code is provided for LU factorization; MUMPS also supports
LDL and LL factorization (then, only the triangular parts must be
stored). The flags to be changed are commented.
- In my opinion, errors of the backend should be made accessible to the
user, also by string representations. I think, so far this is not
implemented in the other backends.
- The sample code is meant as pseudo C++ code. I used some typedefs for
matrices, sparse matrices and scalars. TScalar is double, Matrix is a
dynamically sized column-major double precision Eigen mantrix,
SparseMatrix is a double precision column-major Eigen SparseMatrix.
Best regards
Sebastian
Sample code:
Header:
#!C++
#include <Eigen/Core>
#include <Eigen/Sparse>
extern "C"
{
#include <dmumps_c.h>
#include <smumps_c.h>
}
template <typename T> struct TMumpsInterface {};
template <> struct TMumpsInterface<float> {
typedef SMUMPS_STRUC_C MUMPS_STRUC_C;
static void mumps_c(MUMPS_STRUC_C &info) { smumps_c(&info); }
};
template <> struct TMumpsInterface<double> {
typedef DMUMPS_STRUC_C MUMPS_STRUC_C;
static void mumps_c(MUMPS_STRUC_C &info) { dmumps_c(&info); }
};
/**
\brief LU solver for sparse matrices using MUMPS.
*/
class SparseLUMUMPS : public SparseLU
{
public:
// constructors ------------------------------------------------------
//! destructor
virtual ~SparseLUMUMPS();
//! default constructor
SparseLUMUMPS();
virtual void SetFlags(int flag);
virtual void SetPrecision(double prec);
virtual bool Solve(const MatrixBase &b, MatrixBase &x); // no const
function !!!
virtual bool SolveInPlace(MatrixBase &b); // no const function !!!
virtual bool Compute(const SparseMatrixBase &);
private:
typedef TMumpsInterface<TScalar> TSolverInterface;
typedef TSolverInterface::MUMPS_STRUC_C TSolverData;
TSolverData m_solver;
bool m_solvercontainsdata;
void init();
void clear();
const std::string mumpserror(int id) const;
};
Source:
#!C++
SparseLUMUMPS::SparseLUMUMPS() : SparseLU(), m_solvercontainsdata(false)
{
init();
}
SparseLUMUMPS::~SparseLUMUMPS()
{
clear();
}
void SparseLUMUMPS::init()
{
m_solver.comm_fortran = -987654; // for MPI (also required if MPI=disabled)
m_solver.par = 1; // will be used for factorization AND solution
m_solver.sym = 0; // matrix is not symmetric //m_solver.sym = 1; //
matrix is symmetric and positive definite
//m_solver.sym = 2; // matrix is symmetric and not positive definite
m_solver.job = -1; // initialize solver
m_solver.irn = 0; // rows
m_solver.jcn = 0; // columns
m_solver.a = 0; // values
m_solvercontainsdata = false;
succeeded_factorization() = false;
}
void SparseLUMUMPS::clear()
{
if (m_solvercontainsdata)
{
m_solver.job = -2;
TSolverInterface::mumps_c(m_solver);
if (m_solver.info[0] < 0)
{
std::stringstream msg;
msg << "MUMPS error at termination: " <<
mumpserror(m_solver.info[0]) << ".";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
m_solvercontainsdata = false;
succeeded_factorization() = false;
// also delete allocated data:
if (m_solver.irn != 0) delete[] m_solver.irn;
if (m_solver.jcn != 0) delete[] m_solver.jcn;
if (m_solver.a != 0) delete[] m_solver.a;
m_solver.irn = 0; // rows
m_solver.jcn = 0; // columns
m_solver.a = 0; // values
}
}
const std::string SparseLUMUMPS::mumpserror(int id) const
{
switch (id) {
case -6: case -10:
return "matrix is singular";
case -13 :
return "not enough memory";
case -9:
{
std::stringstream msg;
msg << "error " << id << ", increase ICNTL(14)";
return msg.str();
}
default:
{
std::stringstream msg;
msg << "error " << id;
return msg.str();
}
}
}
bool SparseLUMUMPS::Compute(const SparseMatrixBase &A)
{
#define ICNTL(I) icntl[(I)-1] if ( A.rows() != A.cols() )
{
std::stringstream msg;
msg << "Can not factorize a nonsquare matrix ("<<A.rows()<<" x
"<<A.cols()<<")";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
m_dimensions = std::pair<int,int>(A.rows(),A.cols()); // save
dimensions of matrix for later checks
clear(); // delete allocated data from previous run
init(); // reinit MUMPS engine parameters
TSolverInterface::mumps_c(m_solver); // let MUMPS initialize itself.
// create data vectors int * rows(0);
int * cols(0);
TScalar * values(0);
if (A.nonZeros() != 0)
{
try
{
rows = new int[A.nonZeros()];
cols = new int[A.nonZeros()];
values = new TScalar[A.nonZeros()];
}
GLOBAL_CATCH_THROW (logger_tng_tngmath, "Out of memory (could
not allocate work space).")
}
unsigned long int counter_value(0);
for (int k=0; k<A.outerSize(); ++k)
for (TSparseMatrix::InnerIterator it(A,k); it; ++it)
{
// if (it.row() > it.col()) continue; // upper triangular + diagonal
values[counter_value] = it.value();
rows [counter_value] = it.row()+1;
cols [counter_value] = it.col()+1;
++counter_value;
}
m_solvercontainsdata = true;
// prepare m_solver.n = A.rows();
m_solver.nz = counter_value;//A.nonZeros();
m_solver.irn = rows; // row indices
m_solver.jcn = cols; // column indices
m_solver.a = values; // values
m_solver.rhs = 0; // right hand side will be set later
// solver params
m_solver.icntl[0] = 0; // output errors
m_solver.icntl[1] = 0; // output detailed messages
m_solver.icntl[2] = 0; // output info messages
m_solver.icntl[3] = 0; // verbose level
m_solver.icntl[5] = 7; // automatically find out if to scale or
permute m_solver.icntl[6] = 7; // automatic pivot order for the
factorization m_solver.icntl[7] = 7; // simultaneous row and colum
iterative scaling (for not symmetric matrices)
//m_solver.icntl[7] = 1; // diagonal scaling (for symmetric matrices)
// m_solver.icntl[0] = 6; // output errors
// m_solver.icntl[1] = 6; // output detailed messages
// m_solver.icntl[2] = 6; // output info messages
// m_solver.icntl[3] = 0; // verbose level
// analysis phase
m_solver.job = 1;
TSolverInterface::mumps_c(m_solver); // call MUMPS
if (m_solver.info[0] < 0)
{
std::stringstream msg;
msg << "MUMPS error during analysis and reordering: " <<
mumpserror(m_solver.info[0]) << ".";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
// factorization
m_solver.job = 2;
TSolverInterface::mumps_c(m_solver); // call MUMPS
if (m_solver.info[0] < 0)
{
std::stringstream msg;
msg << "MUMPS error during factorization: " <<
mumpserror(m_solver.info[0]) << ".";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
/* possible strategy: what may happen if ICNTL(14) was not correctly set?
m_solver.ICNTL(14) /= 2;
int i_iter = 0;
do
{
++i_iter;
m_solver.ICNTL(14) *= 2;
TSolverInterface::mumps_c(m_solver);
}
while(m_solver.info[0] == -9 && i_iter < 3);
*/
// delete temporarily allocated data:
if (m_solver.irn != 0) delete[] m_solver.irn;
if (m_solver.jcn != 0) delete[] m_solver.jcn;
if (m_solver.a != 0) delete[] m_solver.a;
m_solver.irn = 0; // rows
m_solver.jcn = 0; // columns
m_solver.a = 0; // values
m_solver.job = 3; // next job will be "solution"
succeeded_factorization() = true; return succeeded_factorization();
}
bool SparseLUMUMPS::Solve(const MatrixBase &b, MatrixBase &x) {
if (b.rows() != m_dimensions.second)
{
std::stringstream msg;
msg << "Incompatible dimensions.";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
if (!succeeded_factorization() || !m_solvercontainsdata) return false;
x = b;
// MUMPS solution
m_solver.job = 3;
for (int column_rhs = 0; column_rhs != x.cols(); ++column_rhs)
{
m_solver.rhs = (TScalar*) (&(x.data()[column_rhs*x.rows()]));
TSolverInterface::mumps_c(m_solver); // call MUMPS
if (m_solver.info[0] < 0)
{
std::stringstream msg;
msg << "MUMPS error during solution: " <<
mumpserror(m_solver.info[0]) << ".";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
}
m_solver.rhs = 0;
return true;
}
bool SparseLUMUMPS::SolveInPlace(MatrixBase &b) {
if (b.rows() != m_dimensions.second)
{
std::stringstream msg;
msg << "Incompatible dimensions.";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
if (!succeeded_factorization() || !m_solvercontainsdata) return false;
// MUMPS solution
m_solver.job = 3;
for (int column_rhs = 0; column_rhs != b.cols(); ++column_rhs)
{
m_solver.rhs =
const_cast<double*>(&(b.data()[((size_t)column_rhs)*b.rows()]));
TSolverInterface::mumps_c(m_solver); // call MUMPS
if (m_solver.info[0] < 0)
{
std::stringstream msg;
msg << "MUMPS error during solution: " <<
mumpserror(m_solver.info[0]) << ".";
GLOBAL_THROW (logger_tng_tngmath, msg.str())
}
}
m_solver.rhs = 0;
return true;
}
void SparseLUMUMPS::SetFlags(int flag)
{
// solver().setFlags(flag);
}
void SparseLUMUMPS::SetPrecision(double prec)
{
// solver().setPrecision(prec);
}