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

patch from Moritz Lenz to allow solving transposed problem with superlu

Browse Source
This commit is contained in:
Gael Guennebaud 2009-04-10 19:54:43 +00:00
parent fb3078fb62
commit 804a239d30
3 changed files with 39 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
Eigen/src/Sparse/SparseLU.h
View File

@ -25,6 +25,12 @@
#ifndef EIGEN_SPARSELU_H #ifndef EIGEN_SPARSELU_H
#define EIGEN_SPARSELU_H #define EIGEN_SPARSELU_H
enum {
SvNoTrans = 0,
SvTranspose = 1,
SvAdjoint = 2
};
/** \ingroup Sparse_Module /** \ingroup Sparse_Module
* *
* \class SparseLU * \class SparseLU
@ -115,7 +121,8 @@ class SparseLU
//inline const MatrixType& matrixU() const { return m_matrixU; } //inline const MatrixType& matrixU() const { return m_matrixU; }
template<typename BDerived, typename XDerived> template<typename BDerived, typename XDerived>
bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const; bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
const int transposed = SvNoTrans) const;
/** \returns true if the factorization succeeded */ /** \returns true if the factorization succeeded */
inline bool succeeded(void) const { return m_succeeded; } inline bool succeeded(void) const { return m_succeeded; }
@ -136,10 +143,17 @@ void SparseLU<MatrixType,Backend>::compute(const MatrixType& a)
ei_assert(false && "not implemented yet"); ei_assert(false && "not implemented yet");
} }
/** Computes *x = U^-1 L^-1 b */ /** Computes *x = U^-1 L^-1 b
*
* If \a transpose is set to SvTranspose or SvAdjoint, the solution
* of the transposed/adjoint system is computed instead.
*
* Not all backends implement the solution of the transposed or
* adjoint system.
*/
template<typename MatrixType, int Backend> template<typename MatrixType, int Backend>
template<typename BDerived, typename XDerived> template<typename BDerived, typename XDerived>
bool SparseLU<MatrixType,Backend>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const bool SparseLU<MatrixType,Backend>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const
{ {
ei_assert(false && "not implemented yet"); ei_assert(false && "not implemented yet");
return false; return false;

16
Eigen/src/Sparse/SuperLUSupport.h
View File

@ -318,7 +318,7 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
Scalar determinant() const; Scalar determinant() const;
template<typename BDerived, typename XDerived> template<typename BDerived, typename XDerived>
bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const; bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed = SvNoTrans) const;
void compute(const MatrixType& matrix); void compute(const MatrixType& matrix);
@ -413,12 +413,22 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
template<typename MatrixType> template<typename MatrixType>
template<typename BDerived,typename XDerived> template<typename BDerived,typename XDerived>
bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b,
MatrixBase<XDerived> *x, const int transposed) const
{ {
const int size = m_matrix.rows(); const int size = m_matrix.rows();
const int rhsCols = b.cols(); const int rhsCols = b.cols();
ei_assert(size==b.rows()); ei_assert(size==b.rows());
switch (transposed) {
case SvNoTrans : m_sluOptions.Trans = NOTRANS; break;
case SvTranspose : m_sluOptions.Trans = TRANS; break;
case SvAdjoint : m_sluOptions.Trans = CONJ; break;
default:
std::cerr << "Eigen: tranpsiotion option \"" << transposed << "\" not supported by the SuperLU backend\n";
m_sluOptions.Trans = NOTRANS;
}
m_sluOptions.Fact = FACTORED; m_sluOptions.Fact = FACTORED;
m_sluOptions.IterRefine = NOREFINE; m_sluOptions.IterRefine = NOREFINE;
@ -443,6 +453,8 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBa
&m_sluStat, &info, Scalar()); &m_sluStat, &info, Scalar());
StatFree(&m_sluStat); StatFree(&m_sluStat);
// reset to previous state
m_sluOptions.Trans = NOTRANS;
return info==0; return info==0;
} }

8
test/sparse_solvers.cpp
View File

@ -191,6 +191,14 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU");
} }
// std::cerr << refDet << " == " << slu.determinant() << "\n"; // std::cerr << refDet << " == " << slu.determinant() << "\n";
if (slu.solve(b, &x, SvTranspose)) {
VERIFY(b.isApprox(m2.transpose() * x, test_precision<Scalar>()));
}
if (slu.solve(b, &x, SvAdjoint)) {
// VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>()));
}
if (count==0) { if (count==0) {
VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
} }