GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-15 01:21:48 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
(cherry picked from commit ab2b26fbc27cd03c1d75ea8c2cce22fdd2bcc45b)
This commit is contained in:
Antonio Sánchez 2022-11-29 19:37:03 +00:00 committed by Antonio Sanchez
parent 5547205092
commit 26adb0e5af
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

@ -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<typename Dest> void solveInPlace(MatrixBase<Dest> &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<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);
}
@ -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<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)
U = A.adjoint().template triangularView<Lower>().solve(U);
else

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

@ -274,9 +274,8 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) co
// Triangular solve
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);
// Matrix-vector product
new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
work.topRows(nrow).noalias() = A * U;
@ -349,7 +348,7 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveTransposedInPlace( MatrixBase<D
// Matrix-vector product with transposed submatrix
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)
U = U - A.adjoint() * work.topRows(nrow);
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(oldb.isApprox(bm) && "sparse solver testing: the rhs should not be modified!");
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: