diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h index 142592050..9cab67a51 100644 --- a/Eigen/src/Sparse/SparseLU.h +++ b/Eigen/src/Sparse/SparseLU.h @@ -25,6 +25,12 @@ #ifndef EIGEN_SPARSELU_H #define EIGEN_SPARSELU_H +enum { + SvNoTrans = 0, + SvTranspose = 1, + SvAdjoint = 2 +}; + /** \ingroup Sparse_Module * * \class SparseLU @@ -115,7 +121,8 @@ class SparseLU //inline const MatrixType& matrixU() const { return m_matrixU; } template - bool solve(const MatrixBase &b, MatrixBase* x) const; + bool solve(const MatrixBase &b, MatrixBase* x, + const int transposed = SvNoTrans) const; /** \returns true if the factorization succeeded */ inline bool succeeded(void) const { return m_succeeded; } @@ -136,10 +143,17 @@ void SparseLU::compute(const MatrixType& a) 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 template -bool SparseLU::solve(const MatrixBase &b, MatrixBase* x) const +bool SparseLU::solve(const MatrixBase &b, MatrixBase* x, const int transposed) const { ei_assert(false && "not implemented yet"); return false; diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h index 6df94e35b..2aeccc534 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -318,7 +318,7 @@ class SparseLU : public SparseLU Scalar determinant() const; template - bool solve(const MatrixBase &b, MatrixBase* x) const; + bool solve(const MatrixBase &b, MatrixBase* x, const int transposed = SvNoTrans) const; void compute(const MatrixType& matrix); @@ -413,12 +413,22 @@ void SparseLU::compute(const MatrixType& a) template template -bool SparseLU::solve(const MatrixBase &b, MatrixBase *x) const +bool SparseLU::solve(const MatrixBase &b, + MatrixBase *x, const int transposed) const { const int size = m_matrix.rows(); const int rhsCols = b.cols(); 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.IterRefine = NOREFINE; @@ -443,6 +453,8 @@ bool SparseLU::solve(const MatrixBase &b, MatrixBa &m_sluStat, &info, Scalar()); StatFree(&m_sluStat); + // reset to previous state + m_sluOptions.Trans = NOTRANS; return info==0; } diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index d1090dfed..e1ec1ef35 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -191,6 +191,14 @@ template void sparse_solvers(int rows, int cols) VERIFY(refX.isApprox(x,test_precision()) && "LU: SuperLU"); } // std::cerr << refDet << " == " << slu.determinant() << "\n"; + if (slu.solve(b, &x, SvTranspose)) { + VERIFY(b.isApprox(m2.transpose() * x, test_precision())); + } + + if (slu.solve(b, &x, SvAdjoint)) { +// VERIFY(b.isApprox(m2.adjoint() * x, test_precision())); + } + if (count==0) { VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex }