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

UmfPack support: fix redundant evaluation/copies when calling compute() and support generic expressions as input

Browse Source
This commit is contained in:
Gael Guennebaud 2014-12-02 17:30:57 +01:00
parent d66b5a1d91
commit 97812ad0d3
1 changed files with 62 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
Eigen/src/UmfPackSupport/UmfPackSupport.h
View File

@ -107,6 +107,16 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N
return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
} }
namespace internal {
template<typename T> struct umfpack_helper_is_sparse_plain : false_type {};
template<typename Scalar, int Options, typename StorageIndex>
struct umfpack_helper_is_sparse_plain<SparseMatrix<Scalar,Options,StorageIndex> >
: true_type {};
template<typename Scalar, int Options, typename StorageIndex>
struct umfpack_helper_is_sparse_plain<MappedSparseMatrix<Scalar,Options,StorageIndex> >
: true_type {};
}
/** \ingroup UmfPackSupport_Module /** \ingroup UmfPackSupport_Module
* \brief A sparse LU factorization and solver based on UmfPack * \brief A sparse LU factorization and solver based on UmfPack
* *
@ -195,8 +205,11 @@ class UmfPackLU : internal::noncopyable
template<typename InputMatrixType> template<typename InputMatrixType>
void compute(const InputMatrixType& matrix) void compute(const InputMatrixType& matrix)
{ {
analyzePattern(matrix); if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
factorize(matrix); if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
grapInput(matrix.derived());
analyzePattern_impl();
factorize_impl();
} }
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
@ -234,22 +247,12 @@ class UmfPackLU : internal::noncopyable
template<typename InputMatrixType> template<typename InputMatrixType>
void analyzePattern(const InputMatrixType& matrix) void analyzePattern(const InputMatrixType& matrix)
{ {
if(m_symbolic) if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
umfpack_free_symbolic(&m_symbolic,Scalar()); if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
if(m_numeric)
umfpack_free_numeric(&m_numeric,Scalar());
grapInput(matrix); grapInput(matrix.derived());
int errorCode = 0; analyzePattern_impl();
errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
&m_symbolic, 0, 0);
m_isInitialized = true;
m_info = errorCode ? InvalidInput : Success;
m_analysisIsOk = true;
m_factorizationIsOk = false;
m_extractedDataAreDirty = true;
} }
/** Performs a numeric decomposition of \a matrix /** Performs a numeric decomposition of \a matrix
@ -265,15 +268,9 @@ class UmfPackLU : internal::noncopyable
if(m_numeric) if(m_numeric)
umfpack_free_numeric(&m_numeric,Scalar()); umfpack_free_numeric(&m_numeric,Scalar());
grapInput(matrix); grapInput(matrix.derived());
int errorCode; factorize_impl();
errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
m_symbolic, &m_numeric, 0, 0);
m_info = errorCode ? NumericalIssue : Success;
m_factorizationIsOk = true;
m_extractedDataAreDirty = true;
} }
#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_PARSED_BY_DOXYGEN
@ -288,7 +285,6 @@ class UmfPackLU : internal::noncopyable
protected: protected:
void init() void init()
{ {
m_info = InvalidInput; m_info = InvalidInput;
@ -302,7 +298,7 @@ class UmfPackLU : internal::noncopyable
} }
template<typename InputMatrixType> template<typename InputMatrixType>
void grapInput(const InputMatrixType& mat) void grapInput_impl(const InputMatrixType& mat, internal::true_type)
{ {
m_copyMatrix.resize(mat.rows(), mat.cols()); m_copyMatrix.resize(mat.rows(), mat.cols());
if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() ) if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() )
@ -320,6 +316,45 @@ class UmfPackLU : internal::noncopyable
m_valuePtr = mat.valuePtr(); m_valuePtr = mat.valuePtr();
} }
} }
template<typename InputMatrixType>
void grapInput_impl(const InputMatrixType& mat, internal::false_type)
{
m_copyMatrix = mat;
m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
m_valuePtr = m_copyMatrix.valuePtr();
}
template<typename InputMatrixType>
void grapInput(const InputMatrixType& mat)
{
grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain<InputMatrixType>());
}
void analyzePattern_impl()
{
int errorCode = 0;
errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
&m_symbolic, 0, 0);
m_isInitialized = true;
m_info = errorCode ? InvalidInput : Success;
m_analysisIsOk = true;
m_factorizationIsOk = false;
m_extractedDataAreDirty = true;
}
void factorize_impl()
{
int errorCode;
errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
m_symbolic, &m_numeric, 0, 0);
m_info = errorCode ? NumericalIssue : Success;
m_factorizationIsOk = true;
m_extractedDataAreDirty = true;
}
// cached data to reduce reallocation, etc. // cached data to reduce reallocation, etc.
mutable LUMatrixType m_l; mutable LUMatrixType m_l;