devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

add specialization of check_sparse_solving() for SuperLU solver, in order to test adjoint and transpose solves

Browse Source
This commit is contained in:
Ralf Hannemann-Tamas
2021-02-08 22:00:31 +00:00
committed by Rasmus Munk Larsen
parent b578930657
commit 984d010b7b
4 changed files with 357 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -18,6 +18,63 @@ template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename
template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType;
template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType;
template <bool Conjugate,class SparseLUType>
class SparseLUTransposeView : public SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> >
{
protected:
typedef SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> > APIBase;
using APIBase::m_isInitialized;
public:
typedef typename SparseLUType::Scalar Scalar;
typedef typename SparseLUType::StorageIndex StorageIndex;
typedef typename SparseLUType::MatrixType MatrixType;
typedef typename SparseLUType::OrderingType OrderingType;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
SparseLUTransposeView() : m_sparseLU(nullptr) {}
SparseLUTransposeView(const SparseLUTransposeView& view) {
this->m_sparseLU = view.m_sparseLU;
}
void setIsInitialized(const bool isInitialized) {this->m_isInitialized = isInitialized;}
void setSparseLU(SparseLUType* sparseLU) {m_sparseLU = sparseLU;}
using APIBase::_solve_impl;
template<typename Rhs, typename Dest>
bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const
{
Dest& X(X_base.derived());
eigen_assert(m_sparseLU->info() == Success && "The matrix should be factorized first");
EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
// this ugly const_cast_derived() helps to detect aliasing when applying the permutations
for(Index j = 0; j < B.cols(); ++j){
X.col(j) = m_sparseLU->colsPermutation() * B.const_cast_derived().col(j);
}
//Forward substitution with transposed or adjoint of U
m_sparseLU->matrixU().template solveTransposedInPlace<Conjugate>(X);
//Backward substitution with transposed or adjoint of L
m_sparseLU->matrixL().template solveTransposedInPlace<Conjugate>(X);
// Permute back the solution
for (Index j = 0; j < B.cols(); ++j)
X.col(j) = m_sparseLU->rowsPermutation().transpose() * X.col(j);
return true;
}
inline Index rows() const { return m_sparseLU->rows(); }
inline Index cols() const { return m_sparseLU->cols(); }
private:
SparseLUType *m_sparseLU;
SparseLUTransposeView& operator=(const SparseLUTransposeView&);
};
/** \ingroup SparseLU_Module
* \class SparseLU
*
@@ -97,6 +154,7 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
};
public:
SparseLU():m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
{
initperfvalues();
@@ -128,6 +186,45 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
//Factorize
factorize(matrix);
}
/** \returns an expression of the transposed of the factored matrix.
*
* A typical usage is to solve for the transposed problem A^T x = b:
* \code
* solver.compute(A);
* x = solver.transpose().solve(b);
* \endcode
*
* \sa adjoint(), solve()
*/
const SparseLUTransposeView<false,SparseLU<_MatrixType,_OrderingType> > transpose()
{
SparseLUTransposeView<false, SparseLU<_MatrixType,_OrderingType> > transposeView;
transposeView.setSparseLU(this);
transposeView.setIsInitialized(this->m_isInitialized);
return transposeView;
}
/** \returns an expression of the adjoint of the factored matrix
*
* A typical usage is to solve for the adjoint problem A' x = b:
* \code
* solver.compute(A);
* x = solver.adjoint().solve(b);
* \endcode
*
* For real scalar types, this function is equivalent to transpose().
*
* \sa transpose(), solve()
*/
const SparseLUTransposeView<true, SparseLU<_MatrixType,_OrderingType> > adjoint()
{
SparseLUTransposeView<true, SparseLU<_MatrixType,_OrderingType> > adjointView;
adjointView.setSparseLU(this);
adjointView.setIsInitialized(this->m_isInitialized);
return adjointView;
}
inline Index rows() const { return m_mat.rows(); }
inline Index cols() const { return m_mat.cols(); }
@@ -394,7 +491,6 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
private:
// Disable copy constructor
SparseLU (const SparseLU& );
}; // End class SparseLU
@@ -712,6 +808,12 @@ struct SparseLUMatrixLReturnType : internal::no_assignment_operator
{
m_mapL.solveInPlace(X);
}
template<bool Conjugate, typename Dest>
void solveTransposedInPlace( MatrixBase<Dest> &X) const
{
m_mapL.template solveTransposedInPlace<Conjugate>(X);
}
const MappedSupernodalType& m_mapL;
};
@@ -766,6 +868,52 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator
}
} // End For U-solve
}
template<bool Conjugate, typename Dest> void solveTransposedInPlace(MatrixBase<Dest> &X) const
{
using numext::conj;
Index nrhs = X.cols();
Index n = X.rows();
// Forward solve with U
for (Index k = 0; k <= m_mapL.nsuper(); 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];
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(jcol, j) -= X(irow, j) * (Conjugate? conj(it.value()): it.value());
}
}
}
if (nsupc == 1)
{
for (Index j = 0; j < nrhs; j++)
{
X(fsupc, j) /= (Conjugate? conj(m_mapL.valuePtr()[luptr]) : m_mapL.valuePtr()[luptr]);
}
}
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) );
if(Conjugate)
U = A.adjoint().template triangularView<Lower>().solve(U);
else
U = A.transpose().template triangularView<Lower>().solve(U);
}
}// End For U-solve
}
const MatrixLType& m_mapL;
const MatrixUType& m_mapU;
};

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

@@ -156,6 +156,9 @@ class MappedSuperNodalMatrix
class InnerIterator;
template<typename Dest>
void solveInPlace( MatrixBase<Dest>&X) const;
template<bool Conjugate, typename Dest>
void solveTransposedInPlace( MatrixBase<Dest>&X) const;
@@ -294,6 +297,77 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) co
}
}
template<typename Scalar, typename Index_>
template<bool Conjugate, typename Dest>
void MappedSuperNodalMatrix<Scalar,Index_>::solveTransposedInPlace( MatrixBase<Dest>&X) const
{
using numext::conj;
Index n = int(X.rows());
Index nrhs = Index(X.cols());
const Scalar * Lval = valuePtr(); // Nonzero values
Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector
work.setZero();
for (Index k = nsuper(); k >= 0; k--)
{
Index fsupc = supToCol()[k]; // First column of the current supernode
Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column
Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode
Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode
Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode
Index irow; //Current index row
if (nsupc == 1 )
{
for (Index j = 0; j < nrhs; j++)
{
InnerIterator it(*this, fsupc);
++it; // Skip the diagonal element
for (; it; ++it)
{
irow = it.row();
X(fsupc,j) -= X(irow, j) * (Conjugate?conj(it.value()):it.value());
}
}
}
else
{
// The supernode has more than one column
Index luptr = colIndexPtr()[fsupc];
Index lda = colIndexPtr()[fsupc+1] - luptr;
//Begin Gather
for (Index j = 0; j < nrhs; j++)
{
Index iptr = istart + nsupc;
for (Index i = 0; i < nrow; i++)
{
irow = rowIndex()[iptr];
work.topRows(nrow)(i,j)= X(irow,j); // Gather operation
iptr++;
}
}
// 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) );
if(Conjugate)
U = U - A.adjoint() * work.topRows(nrow);
else
U = U - A.transpose() * work.topRows(nrow);
// Triangular solve (of transposed diagonal block)
new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
if(Conjugate)
U = A.adjoint().template triangularView<UnitUpper>().solve(U);
else
U = A.transpose().template triangularView<UnitUpper>().solve(U);
}
}
}
} // end namespace internal
} // end namespace Eigen