From 26adb0e5aff4de53f048a63b382d5e1a1ad2e96b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Antonio=20S=C3=A1nchez?= Date: Tue, 29 Nov 2022 19:37:03 +0000 Subject: [PATCH] Fix sparseLU solver when destination has a non-unit stride. (cherry picked from commit ab2b26fbc27cd03c1d75ea8c2cce22fdd2bcc45b) --- Eigen/src/SparseLU/SparseLU.h | 10 ++++------ Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h | 7 +++---- test/sparse_solver.h | 7 +++++++ 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index f9143bc02..761b95c98 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -35,8 +35,8 @@ public: MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; - SparseLUTransposeView() : m_sparseLU(NULL) {} - SparseLUTransposeView(const SparseLUTransposeView& view) { + SparseLUTransposeView() : APIBase(), m_sparseLU(NULL) {} + SparseLUTransposeView(const SparseLUTransposeView& view) : APIBase() { this->m_sparseLU = view.m_sparseLU; } void setIsInitialized(const bool isInitialized) {this->m_isInitialized = isInitialized;} @@ -833,7 +833,6 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator 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--) { @@ -853,7 +852,7 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator { // FIXME: the following lines should use Block expressions and not Map! Map, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); - Map< Matrix, 0, OuterStride<> > U (&(X.coeffRef(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); + typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc); U = A.template triangularView().solve(U); } @@ -876,7 +875,6 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator { using numext::conj; Index nrhs = X.cols(); - Index n = X.rows(); // Forward solve with U for (Index k = 0; k <= m_mapL.nsuper(); k++) { @@ -907,7 +905,7 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator 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) ); + typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc); if(Conjugate) U = A.adjoint().template triangularView().solve(U); else diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h index 0be293d17..fd5e9fa51 100644 --- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h +++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h @@ -274,9 +274,8 @@ void MappedSuperNodalMatrix::solveInPlace( MatrixBase&X) co // Triangular solve Map, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); - Map< Matrix, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); - U = A.template triangularView().solve(U); - + typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc); + U = A.template triangularView().solve(U); // Matrix-vector product new (&A) Map, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); work.topRows(nrow).noalias() = A * U; @@ -349,7 +348,7 @@ void MappedSuperNodalMatrix::solveTransposedInPlace( MatrixBase, 0, OuterStride<> > A( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); - Map< Matrix, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); + typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc); if(Conjugate) U = U - A.adjoint() * work.topRows(nrow); else diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 58927944b..6f95e2fa7 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -99,6 +99,13 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, VERIFY(solver.info() == Success && "solving failed when using Map"); VERIFY(oldb.isApprox(bm) && "sparse solver testing: the rhs should not be modified!"); VERIFY(xm.isApprox(refX,test_precision())); + + // Test with a Map and non-unit stride. + Eigen::Matrix out(2*xm.rows(), 2*xm.cols()); + out.setZero(); + Eigen::Map > outm(out.data(), xm.rows(), xm.cols(), Stride(2 * xm.rows(), 2)); + outm = solver.solve(bm); + VERIFY(outm.isApprox(refX,test_precision())); } // if not too large, do some extra check: