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

Fix sparseLU solver when destination has a non-unit stride.

Browse Source
This commit is contained in:
Antonio Sánchez 2022-11-29 19:37:03 +00:00 committed by Rasmus Munk Larsen
parent 551eebc8ca
commit ab2b26fbc2
3 changed files with 14 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Eigen/src/SparseLU/SparseLU.h
View File

@ -37,8 +37,8 @@ public:
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
}; };
SparseLUTransposeView() : m_sparseLU(NULL) {} SparseLUTransposeView() : APIBase(), m_sparseLU(NULL) {}
SparseLUTransposeView(const SparseLUTransposeView& view) { SparseLUTransposeView(const SparseLUTransposeView& view) : APIBase() {
this->m_sparseLU = view.m_sparseLU; this->m_sparseLU = view.m_sparseLU;
} }
void setIsInitialized(const bool isInitialized) {this->m_isInitialized = isInitialized;} void setIsInitialized(const bool isInitialized) {this->m_isInitialized = isInitialized;}
@ -860,7 +860,6 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator
template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const
{ {
Index nrhs = X.cols(); Index nrhs = X.cols();
Index n = X.rows();
// Backward solve with U // Backward solve with U
for (Index k = m_mapL.nsuper(); k >= 0; k--) for (Index k = m_mapL.nsuper(); k >= 0; k--)
{ {
@ -880,7 +879,7 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator
{ {
// FIXME: the following lines should use Block expressions and not Map! // FIXME: the following lines should use Block expressions and not Map!
Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X.coeffRef(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc);
U = A.template triangularView<Upper>().solve(U); U = A.template triangularView<Upper>().solve(U);
} }
@ -903,7 +902,6 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator
{ {
using numext::conj; using numext::conj;
Index nrhs = X.cols(); Index nrhs = X.cols();
Index n = X.rows();
// Forward solve with U // Forward solve with U
for (Index k = 0; k <= m_mapL.nsuper(); k++) for (Index k = 0; k <= m_mapL.nsuper(); k++)
{ {
@ -934,7 +932,7 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator
else else
{ {
Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc);
if(Conjugate) if(Conjugate)
U = A.adjoint().template triangularView<Lower>().solve(U); U = A.adjoint().template triangularView<Lower>().solve(U);
else else

5
Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
View File

@ -276,9 +276,8 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) co
// Triangular solve // Triangular solve
Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc);
U = A.template triangularView<UnitLower>().solve(U); U = A.template triangularView<UnitLower>().solve(U);
// Matrix-vector product // Matrix-vector product
new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
work.topRows(nrow).noalias() = A * U; work.topRows(nrow).noalias() = A * U;
@ -351,7 +350,7 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveTransposedInPlace( MatrixBase<D
// Matrix-vector product with transposed submatrix // Matrix-vector product with transposed submatrix
Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); typename Dest::RowsBlockXpr U = X.derived().middleRows(fsupc, nsupc);
if(Conjugate) if(Conjugate)
U = U - A.adjoint() * work.topRows(nrow); U = U - A.adjoint() * work.topRows(nrow);
else else

7
test/sparse_solver.h
View File

@ -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(solver.info() == Success && "solving failed when using Map");
VERIFY(oldb.isApprox(bm) && "sparse solver testing: the rhs should not be modified!"); VERIFY(oldb.isApprox(bm) && "sparse solver testing: the rhs should not be modified!");
VERIFY(xm.isApprox(refX,test_precision<Scalar>())); VERIFY(xm.isApprox(refX,test_precision<Scalar>()));
// Test with a Map and non-unit stride.
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> out(2*xm.rows(), 2*xm.cols());
out.setZero();
Eigen::Map<DenseRhs, 0, Stride<Eigen::Dynamic, 2>> outm(out.data(), xm.rows(), xm.cols(), Stride<Eigen::Dynamic, 2>(2 * xm.rows(), 2));
outm = solver.solve(bm);
VERIFY(outm.isApprox(refX,test_precision<Scalar>()));
} }
// if not too large, do some extra check: // if not too large, do some extra check: