GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-22 01:29:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

sparse module: add preliminary support for direct sparse LU solver

Browse Source 
using SuperLU. Calling SuperLU was very painful, but it was worth it,
it seems to be damn fast !
This commit is contained in:
Gael Guennebaud 2008-10-19 15:26:28 +00:00
parent 6be0131774
commit ecc6c43dba
5 changed files with 562 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
Eigen/Sparse
View File

@ -29,6 +29,27 @@
#endif
#ifdef EIGEN_SUPERLU_SUPPORT
typedef int int_t;
#include "superlu/slu_Cnames.h"
#include "superlu/supermatrix.h"
#include "superlu/slu_util.h"
namespace SuperLU_S {
#include "superlu/slu_sdefs.h"
}
namespace SuperLU_D {
#include "superlu/slu_ddefs.h"
}
namespace SuperLU_C {
#include "superlu/slu_cdefs.h"
}
namespace SuperLU_Z {
#include "superlu/slu_zdefs.h"
}
namespace Eigen { struct SluMatrix; }
#endif
namespace Eigen {
#include "src/Sparse/SparseUtil.h"
@ -54,6 +75,10 @@ namespace Eigen {
# include "src/Sparse/TaucsSupport.h"
#endif
#ifdef EIGEN_SUPERLU_SUPPORT
# include "src/Sparse/SuperLUSupport.h"
#endif
} // namespace Eigen
#endif // EIGEN_SPARSE_MODULE_H

148
Eigen/src/Sparse/SparseLU.h Normal file
View File

@ -0,0 +1,148 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_SPARSELU_H
#define EIGEN_SPARSELU_H
/** \ingroup Sparse_Module
*
* \class SparseLU
*
* \brief LU decomposition of a sparse matrix and associated features
*
* \param MatrixType the type of the matrix of which we are computing the LU factorization
*
* \sa class LU, class SparseLLT
*/
template<typename MatrixType, int Backend = DefaultBackend>
class SparseLU
{
protected:
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef SparseMatrix<Scalar,Lower> LUMatrixType;
enum {
MatrixLUIsDirty = 0x10000
};
public:
/** Creates a dummy LU factorization object with flags \a flags. */
SparseLU(int flags = 0)
: m_flags(flags), m_status(0)
{
m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
}
/** Creates a LU object and compute the respective factorization of \a matrix using
* flags \a flags. */
SparseLU(const MatrixType& matrix, int flags = 0)
: /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0)
{
m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
compute(matrix);
}
/** Sets the relative threshold value used to prune zero coefficients during the decomposition.
*
* Setting a value greater than zero speeds up computation, and yields to an imcomplete
* factorization with fewer non zero coefficients. Such approximate factors are especially
* useful to initialize an iterative solver.
*
* Note that the exact meaning of this parameter might depends on the actual
* backend. Moreover, not all backends support this feature.
*
* \sa precision() */
void setPrecision(RealScalar v) { m_precision = v; }
/** \returns the current precision.
*
* \sa setPrecision() */
RealScalar precision() const { return m_precision; }
/** Sets the flags. Possible values are:
* - CompleteFactorization
* - IncompleteFactorization
* - MemoryEfficient
* - one of the ordering methods
* - etc...
*
* \sa flags() */
void setFlags(int f) { m_flags = f; }
/** \returns the current flags */
int flags() const { return m_flags; }
void setOrderingMethod(int m)
{
ei_assert(m&~OrderingMask == 0 && m!=0 && "invalid ordering method");
m_flags = m_flags&~OrderingMask | m&OrderingMask;
}
int orderingMethod() const
{
return m_flags&OrderingMask;
}
/** Computes/re-computes the LU factorization */
void compute(const MatrixType& matrix);
/** \returns the lower triangular matrix L */
//inline const MatrixType& matrixL() const { return m_matrixL; }
/** \returns the upper triangular matrix U */
//inline const MatrixType& matrixU() const { return m_matrixU; }
template<typename BDerived, typename XDerived>
bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
/** \returns true if the factorization succeeded */
inline bool succeeded(void) const { return m_succeeded; }
protected:
RealScalar m_precision;
int m_flags;
mutable int m_status;
bool m_succeeded;
};
/** Computes / recomputes the LU decomposition of matrix \a a
* using the default algorithm.
*/
template<typename MatrixType, int Backend>
void SparseLU<MatrixType,Backend>::compute(const MatrixType& a)
{
ei_assert(false && "not implemented yet");
}
/** Computes *x = U^-1 L^-1 b */
template<typename MatrixType, int Backend>
template<typename BDerived, typename XDerived>
bool SparseLU<MatrixType,Backend>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const
{
ei_assert(false && "not implemented yet");
return false;
}
#endif // EIGEN_SPARSELU_H

7
Eigen/src/Sparse/SparseMatrix.h
View File

@ -247,7 +247,7 @@ class SparseMatrix
typedef typename ei_cleantype<OtherCopy>::type _OtherCopy;
resize(other.rows(), other.cols());
Map<VectorXi>(m_outerIndex,outerSize()).setZero();
Eigen::Map<VectorXi>(m_outerIndex,outerSize()).setZero();
// pass 1
// FIXME the above copy could be merged with that pass
for (int j=0; j<otherCopy.outerSize(); ++j)
@ -317,6 +317,11 @@ class SparseMatrix
cholmod_sparse asCholmodMatrix();
#endif
#ifdef EIGEN_SUPERLU_SUPPORT
static SparseMatrix Map(SluMatrix& sluMatrix);
SluMatrix asSluMatrix();
#endif
/** Destructor */
inline ~SparseMatrix()
{

21
Eigen/src/Sparse/SparseUtil.h
View File

@ -41,12 +41,23 @@ enum SparseBackend {
// solver flags
enum {
CompleteFactorization = 0x0000, // the default
IncompleteFactorization = 0x0001,
MemoryEfficient = 0x0002,
CompleteFactorization = 0x0000, // the default
IncompleteFactorization = 0x0001,
MemoryEfficient = 0x0002,
// For LLT Cholesky:
SupernodalMultifrontal = 0x0010,
SupernodalLeftLooking = 0x0020
SupernodalMultifrontal = 0x0010,
SupernodalLeftLooking = 0x0020,
// Ordering methods:
NaturalOrdering = 0x0100, // the default
MinimumDegree_AT_PLUS_A = 0x0200,
MinimumDegree_ATA = 0x0300,
ColApproxMinimumDegree = 0x0400,
Metis = 0x0500,
Scotch = 0x0600,
Chaco = 0x0700,
OrderingMask = 0x0f00
};
template<typename Derived> class SparseMatrixBase;

367
Eigen/src/Sparse/SuperLUSupport.h Normal file
View File

@ -0,0 +1,367 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_SUPERLUSUPPORT_H
#define EIGEN_SUPERLUSUPPORT_H
// declaration of gssvx taken from GMM++
#define DECL_GSSVX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
int *perm_c, int *perm_r, int *etree, char *equed, \
FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
SuperMatrix *U, void *work, int lwork, \
SuperMatrix *B, SuperMatrix *X, \
FLOATTYPE *recip_pivot_growth, \
FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
SuperLUStat_t *stats, int *info, KEYTYPE) { \
NAMESPACE::mem_usage_t mem_usage; \
NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
U, work, lwork, B, X, recip_pivot_growth, rcond, \
ferr, berr, &mem_usage, stats, info); \
return mem_usage.for_lu; /* bytes used by the factor storage */ \
}
DECL_GSSVX(SuperLU_S,sgssvx,float,float)
DECL_GSSVX(SuperLU_C,cgssvx,float,std::complex<float>)
DECL_GSSVX(SuperLU_D,dgssvx,double,double)
DECL_GSSVX(SuperLU_Z,zgssvx,double,std::complex<double>)
template<typename MatrixType>
struct SluMatrixMapHelper;
/** \internal
*
* A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
* and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
*
* This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
*/
struct SluMatrix : SuperMatrix
{
SluMatrix() {}
SluMatrix(const SluMatrix& other)
: SuperMatrix(other)
{
Store = &nnz;
nnz = other.nnz;
values = other.values;
innerInd = other.innerInd;
outerInd = other.outerInd;
}
struct
{
union {int nnz;int lda;};
void *values;
int *innerInd;
int *outerInd;
};
void setStorageType(Stype_t t)
{
Stype = t;
if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
Store = &nnz;
else
{
ei_assert(false && "storage type not supported");
Store = 0;
}
}
template<typename Scalar>
void setScalarType()
{
if (ei_is_same_type<Scalar,float>::ret)
Dtype = SLU_S;
else if (ei_is_same_type<Scalar,double>::ret)
Dtype = SLU_D;
else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
Dtype = SLU_C;
else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
Dtype = SLU_Z;
else
{
ei_assert(false && "Scalar type not supported by SuperLU");
}
}
template<typename MatrixType>
static SluMatrix Map(MatrixType& mat)
{
SluMatrix res;
SluMatrixMapHelper<MatrixType>::run(mat, res);
return res;
}
};
template<typename Scalar, int Rows, int Cols, int StorageOrder, int MRows, int MCols>
struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,StorageOrder,MRows,MCols> >
{
typedef Matrix<Scalar,Rows,Cols,StorageOrder,MRows,MCols> MatrixType;
static void run(MatrixType& mat, SluMatrix& res)
{
assert(StorageOrder==0 && "row-major dense matrices is not supported by SuperLU");
res.setStorageType(SLU_DN);
res.setScalarType<Scalar>();
res.Mtype = SLU_GE;
res.nrow = mat.rows();
res.ncol = mat.cols();
res.lda = mat.stride();
res.values = mat.data();
}
};
template<typename Scalar, int Flags>
struct SluMatrixMapHelper<SparseMatrix<Scalar,Flags> >
{
typedef SparseMatrix<Scalar,Flags> MatrixType;
static void run(MatrixType& mat, SluMatrix& res)
{
if (Flags&RowMajorBit)
{
res.setStorageType(SLU_NR);
res.nrow = mat.cols();
res.ncol = mat.rows();
}
else
{
res.setStorageType(SLU_NC);
res.nrow = mat.rows();
res.ncol = mat.cols();
}
res.Mtype = SLU_GE;
res.nnz = mat.nonZeros();
res.values = mat._valuePtr();
res.innerInd = mat._innerIndexPtr();
res.outerInd = mat._outerIndexPtr();
res.setScalarType<Scalar>();
// FIXME the following is not very accurate
if (Flags & Upper)
res.Mtype = SLU_TRU;
if (Flags & Lower)
res.Mtype = SLU_TRL;
if (Flags & SelfAdjoint)
ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU");
}
};
template<typename Scalar, int Flags>
SluMatrix SparseMatrix<Scalar,Flags>::asSluMatrix()
{
return SluMatrix::Map(*this);
}
template<typename Scalar, int Flags>
SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(SluMatrix& sluMat)
{
SparseMatrix res;
if (Flags&RowMajorBit)
{
assert(sluMat.Stype == SLU_NR);
res.m_innerSize = sluMat.ncol;
res.m_outerSize = sluMat.nrow;
}
else
{
assert(sluMat.Stype == SLU_NC);
res.m_innerSize = sluMat.nrow;
res.m_outerSize = sluMat.ncol;
}
res.m_outerIndex = sluMat.outerInd;
SparseArray<Scalar> data = SparseArray<Scalar>::Map(
sluMat.innerInd,
reinterpret_cast<Scalar*>(sluMat.values),
sluMat.outerInd[res.m_outerSize]);
res.m_data.swap(data);
res.markAsRValue();
return res;
}
template<typename MatrixType>
class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
{
protected:
typedef SparseLU<MatrixType> Base;
typedef typename Base::Scalar Scalar;
typedef typename Base::RealScalar RealScalar;
typedef Matrix<Scalar,Dynamic,1> Vector;
using Base::m_flags;
using Base::m_status;
public:
SparseLU(int flags = NaturalOrdering)
: Base(flags)
{
}
SparseLU(const MatrixType& matrix, int flags = NaturalOrdering)
: Base(matrix, flags)
{
compute(matrix);
}
~SparseLU()
{
}
template<typename BDerived, typename XDerived>
bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
void compute(const MatrixType& matrix);
protected:
// cached data to reduce reallocation:
mutable SparseMatrix<Scalar> m_matrix;
mutable SluMatrix m_sluA;
mutable SuperMatrix m_sluL, m_sluU,;
mutable SluMatrix m_sluB, m_sluX;
mutable SuperLUStat_t m_sluStat;
mutable superlu_options_t m_sluOptions;
mutable std::vector<int> m_sluEtree, m_sluPermR, m_sluPermC;
mutable std::vector<RealScalar> m_sluRscale, m_sluCscale;
mutable std::vector<RealScalar> m_sluFerr, m_sluBerr;
mutable char m_sluEqued;
};
template<typename MatrixType>
void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
{
const int size = a.rows();
m_matrix = a;
set_default_options(&m_sluOptions);
m_sluOptions.ColPerm = NATURAL;
m_sluOptions.PrintStat = NO;
m_sluOptions.ConditionNumber = NO;
m_sluOptions.Trans = NOTRANS;
switch (Base::orderingMethod())
{
case NaturalOrdering : m_sluOptions.ColPerm = NATURAL; break;
case MinimumDegree_AT_PLUS_A : m_sluOptions.ColPerm = MMD_AT_PLUS_A; break;
case MinimumDegree_ATA : m_sluOptions.ColPerm = MMD_ATA; break;
case ColApproxMinimumDegree : m_sluOptions.ColPerm = COLAMD; break;
default:
std::cerr << "Eigen: ordering method \"" << Base::orderingMethod() << "\" not supported by the SuperLU backend\n";
m_sluOptions.ColPerm = NATURAL;
};
m_sluA = m_matrix.asSluMatrix();
memset(&m_sluL,0,sizeof m_sluL);
memset(&m_sluU,0,sizeof m_sluU);
m_sluEqued = 'B';
int info = 0;
m_sluPermR.resize(size);
m_sluPermC.resize(size);
m_sluRscale.resize(size);
m_sluCscale.resize(size);
m_sluEtree.resize(size);
RealScalar recip_pivot_gross, rcond;
RealScalar ferr, berr;
// set empty B and X
m_sluB.setStorageType(SLU_DN);
m_sluB.setScalarType<Scalar>();
m_sluB.Mtype = SLU_GE;
m_sluB.values = 0;
m_sluB.nrow = m_sluB.ncol = 0;
m_sluB.lda = size;
m_sluX = m_sluB;
StatInit(&m_sluStat);
SuperLU_gssvx(&m_sluOptions, &m_sluA, &m_sluPermC[0], &m_sluPermR[0], &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&ferr, &berr,
&m_sluStat, &info, Scalar());
StatFree(&m_sluStat);
// FIXME how to check for errors ???
m_succeeded = true;
}
// template<typename MatrixType>
// inline const MatrixType&
// SparseLU<MatrixType,SuperLU>::matrixL() const
// {
// ei_assert(false && "matrixL() is Not supported by the SuperLU backend");
// return m_matrix;
// }
//
// template<typename MatrixType>
// inline const MatrixType&
// SparseLU<MatrixType,SuperLU>::matrixU() const
// {
// ei_assert(false && "matrixU() is Not supported by the SuperLU backend");
// return m_matrix;
// }
template<typename MatrixType>
template<typename BDerived,typename XDerived>
bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const
{
const int size = m_matrix.rows();
const int rhsCols = b.cols();
ei_assert(size==b.rows());
m_sluOptions.Fact = FACTORED;
m_sluOptions.IterRefine = NOREFINE;
m_sluFerr.resize(rhsCols);
m_sluBerr.resize(rhsCols);
m_sluB = SluMatrix::Map(b.const_cast_derived());
m_sluX = SluMatrix::Map(x->derived());
StatInit(&m_sluStat);
int info = 0;
RealScalar recip_pivot_gross, rcond;
SuperLU_gssvx(
&m_sluOptions, &m_sluA,
&m_sluPermC[0], &m_sluPermR[0],
&m_sluEtree[0], &m_sluEqued,
&m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&m_sluFerr[0], &m_sluBerr[0],
&m_sluStat, &info, Scalar());
StatFree(&m_sluStat);
}
#endif // EIGEN_SUPERLUSUPPORT_H