From f7bdbf69e1974cde4bbd57bdfc4691d42cb373c3 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Tue, 14 May 2013 17:15:23 +0200 Subject: [PATCH] Add support in SparseLU to solve with L and U factors independently --- Eigen/src/SparseLU/SparseLU.h | 144 +++++++++++------- .../src/SparseLU/SparseLU_SupernodalMatrix.h | 15 +- 2 files changed, 102 insertions(+), 57 deletions(-) diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index e78250084..5475e17f4 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -15,7 +15,8 @@ namespace Eigen { template class SparseLU; -template struct SparseLUMatrixLReturnType; +template struct SparseLUMatrixLReturnType; +template struct SparseLUMatrixUReturnType; /** \ingroup SparseLU_Module * \class SparseLU * @@ -123,8 +124,8 @@ class SparseLU : public internal::SparseLUImpl(m_Lstore); } + /** \returns an expression of the matrix U, + * The only operation available with this expression is the triangular solve + * \code + * y = b; matrixU().solveInPlace(y); + * \endcode + */ + SparseLUMatrixUReturnType > matrixU() const + { + return SparseLUMatrixUReturnType >(m_Lstore, m_Ustore); + } + + /** + * \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$ + * \sa colsPermutation() + */ + inline const PermutationType& rowsPermutation() const + { + return m_perm_r; + } + /** + * \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$ + * \sa rowsPermutation() + */ + inline const PermutationType& colsPermutation() const + { + return m_perm_c; + } /** Set the threshold used for a diagonal entry to be an acceptable pivot. */ void setPivotThreshold(const RealScalar& thresh) { @@ -195,59 +223,19 @@ class SparseLU : public internal::SparseLUImplmatrixL().solveInPlace(X); - - // Backward solve with U - for (Index k = m_Lstore.nsuper(); k >= 0; k--) - { - Index fsupc = m_Lstore.supToCol()[k]; - Index lda = m_Lstore.colIndexPtr()[fsupc+1] - m_Lstore.colIndexPtr()[fsupc]; // leading dimension - Index nsupc = m_Lstore.supToCol()[k+1] - fsupc; - Index luptr = m_Lstore.colIndexPtr()[fsupc]; - - if (nsupc == 1) - { - for (Index j = 0; j < nrhs; j++) - { - X(fsupc, j) /= m_Lstore.valuePtr()[luptr]; - } - } - else - { - Map, 0, OuterStride<> > A( &(m_Lstore.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); - Map< Matrix, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); - U = A.template triangularView().solve(U); - } - - for (Index j = 0; j < nrhs; ++j) - { - for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) - { - typename MappedSparseMatrix::InnerIterator it(m_Ustore, jcol); - for ( ; it; ++it) - { - Index irow = it.index(); - X(irow, j) -= X(jcol, j) * it.value(); - } - } - } - } // End For U-solve + //Forward substitution with L + this->matrixL().solveInPlace(X); + this->matrixU().solveInPlace(X); // Permute back the solution - for (Index j = 0; j < nrhs; ++j) - X.col(j) = m_perm_c.inverse() * X.col(j); + for (Index j = 0; j < B.cols(); ++j) + X.col(j) = colsPermutation().inverse() * X.col(j); return true; } @@ -584,6 +572,60 @@ struct SparseLUMatrixLReturnType const MappedSupernodalType& m_mapL; }; +template +struct SparseLUMatrixUReturnType +{ + typedef typename MatrixLType::Index Index; + typedef typename MatrixLType::Scalar Scalar; + SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU) + : m_mapL(mapL),m_mapU(mapU) + { } + Index rows() { return m_mapL.rows(); } + Index cols() { return m_mapL.cols(); } + + template void solveInPlace(MatrixBase &X) const + { + Index nrhs = X.cols(); + Index n = X.rows(); + // Backward solve with U + for (Index k = m_mapL.nsuper(); k >= 0; k--) + { + Index fsupc = m_mapL.supToCol()[k]; + Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension + Index nsupc = m_mapL.supToCol()[k+1] - fsupc; + Index luptr = m_mapL.colIndexPtr()[fsupc]; + + if (nsupc == 1) + { + for (Index j = 0; j < nrhs; j++) + { + X(fsupc, j) /= m_mapL.valuePtr()[luptr]; + } + } + else + { + Map, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); + Map< Matrix, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); + U = A.template triangularView().solve(U); + } + + for (Index j = 0; j < nrhs; ++j) + { + for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) + { + typename MatrixUType::InnerIterator it(m_mapU, jcol); + for ( ; it; ++it) + { + Index irow = it.index(); + X(irow, j) -= X(jcol, j) * it.value(); + } + } + } + } // End For U-solve + } + const MatrixLType& m_mapL; + const MatrixUType& m_mapU; +}; namespace internal { template diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h index 3eae95479..3836d1096 100644 --- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h +++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h @@ -189,13 +189,14 @@ class MappedSuperNodalMatrix::InnerIterator m_idval(mat.colIndexPtr()[outer]), m_startidval(m_idval), m_endidval(mat.colIndexPtr()[outer+1]), - m_idrow(mat.rowIndexPtr()[outer]) + m_idrow(mat.rowIndexPtr()[outer]), + m_endidrow(mat.rowIndexPtr()[outer+1]) {} inline InnerIterator& operator++() { m_idval++; m_idrow++; - return *this; + return *this; } inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; } @@ -209,7 +210,8 @@ class MappedSuperNodalMatrix::InnerIterator inline operator bool() const { - return ( (m_idval < m_endidval) && (m_idval >= m_startidval) ); + return ( (m_idval < m_endidval) && (m_idval >= m_startidval) + && (m_idrow < m_endidrow) ); } protected: @@ -220,6 +222,7 @@ class MappedSuperNodalMatrix::InnerIterator const Index m_startidval; // Start of the column value const Index m_endidval; // End of the column value Index m_idrow; //Index to browse the row indices + Index m_endidrow; // End index of row indices of the current column }; /** @@ -248,16 +251,16 @@ void MappedSuperNodalMatrix::solveInPlace( MatrixBase&X) con { for (Index j = 0; j < nrhs; j++) { - InnerIterator it(*this, fsupc); + InnerIterator it(*this, fsupc); ++it; // Skip the diagonal element for (; it; ++it) { irow = it.row(); - X(irow, j) -= X(fsupc, j) * it.value(); + X(irow, j) -= X(fsupc, j) * it.value(); } } } - else + else { // The supernode has more than one column Index luptr = colIndexPtr()[fsupc];