| Re: [eigen] Sample code for MUMPS backend |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen <eigen@xxxxxxxxxxxxxxxxxxx>
- Subject: Re: [eigen] Sample code for MUMPS backend
- From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
- Date: Wed, 20 Jan 2010 09:54:06 +0100
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=+S8sTITrENlRn4O7MHlDiAAC0UMJ11yuuO7Ah17vmus=; b=BhWaXV6aoDNCBIQXPeOGMmvDJLjBIxLgLkOo9fXFDEvcrGW35ujTuBXkbuPPB4Btke n+QNlRzENnwlQawLJMLfOfrzpll8HnAqOWwylPlb/4XdHnr1l1Hak/BXu4L3UvD35xge yiQWC6c2tqSZQbRLVYepEwTJ9+rbRI1JM0lWY=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=j6s7ecQojj/Xp4Km5vsLYObt4/AWDj+ypgX8cTxvNN0AuNu30QlYR9C84bt1jt4pkU Y+2ibLARqGBOL3orByGL9QtvA6NXbzsypZ0xQRlgtCbsVJ7njxt/HYseNANXfNzOUALX eyCtM0dGObGs80hOuaVpASJTdrUtI450AjfnQ=
Hi,
thanks a lot for this piece of code. I will sure integrate it in Eigen
as MUMPS is definitely one of the best libs for solving linear
systems. I don't when though!
Also, I'm thinking that the current design of the sparse solvers is
not necessarily very good. So if you have some additional ideas, feel
free. (I'll start a new topic on that).
gael.
On Mon, Jan 18, 2010 at 11:31 AM, Sebastian Wolff
<sw@xxxxxxxxxxxxxxxxxxxx> wrote:
> 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);
> }
>
>
>
>
>
>
>
>