mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-05-31 10:34:02 +08:00
Add SimplicialNonHermitianLLT and SimplicialNonHermitianLDLT
This commit is contained in:
Maarten Baert
Charles Schlosser
parent
4dccaa587e
commit
35bf6c8edc
@ -58,6 +58,7 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived> {
|
||||
enum { UpLo = internal::traits<Derived>::UpLo };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef typename internal::traits<Derived>::DiagonalScalar DiagonalScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef CholMatrixType const* ConstCholMatrixPtr;
|
||||
@ -114,7 +115,7 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived> {
|
||||
*
|
||||
* \returns a reference to \c *this.
|
||||
*/
|
||||
Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) {
|
||||
Derived& setShift(const DiagonalScalar& offset, const DiagonalScalar& scale = 1) {
|
||||
m_shiftOffset = offset;
|
||||
m_shiftScale = scale;
|
||||
return derived();
|
||||
@ -178,18 +179,18 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived> {
|
||||
|
||||
protected:
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
template <bool DoLDLT>
|
||||
template <bool DoLDLT, bool NonHermitian>
|
||||
void compute(const MatrixType& matrix) {
|
||||
eigen_assert(matrix.rows() == matrix.cols());
|
||||
Index size = matrix.cols();
|
||||
CholMatrixType tmp(size, size);
|
||||
ConstCholMatrixPtr pmat;
|
||||
ordering(matrix, pmat, tmp);
|
||||
ordering<NonHermitian>(matrix, pmat, tmp);
|
||||
analyzePattern_preordered(*pmat, DoLDLT);
|
||||
factorize_preordered<DoLDLT>(*pmat);
|
||||
factorize_preordered<DoLDLT, NonHermitian>(*pmat);
|
||||
}
|
||||
|
||||
template <bool DoLDLT>
|
||||
template <bool DoLDLT, bool NonHermitian>
|
||||
void factorize(const MatrixType& a) {
|
||||
eigen_assert(a.rows() == a.cols());
|
||||
Index size = a.cols();
|
||||
@ -200,28 +201,33 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived> {
|
||||
// If there is no ordering, try to directly use the input matrix without any copy
|
||||
internal::simplicial_cholesky_grab_input<CholMatrixType, MatrixType>::run(a, pmat, tmp);
|
||||
} else {
|
||||
tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
|
||||
internal::permute_symm_to_symm<UpLo, Upper, NonHermitian>(a, tmp, m_P.indices().data());
|
||||
pmat = &tmp;
|
||||
}
|
||||
|
||||
factorize_preordered<DoLDLT>(*pmat);
|
||||
factorize_preordered<DoLDLT, NonHermitian>(*pmat);
|
||||
}
|
||||
|
||||
template <bool DoLDLT>
|
||||
template <bool DoLDLT, bool NonHermitian>
|
||||
void factorize_preordered(const CholMatrixType& a);
|
||||
|
||||
void analyzePattern(const MatrixType& a, bool doLDLT) {
|
||||
template <bool DoLDLT, bool NonHermitian>
|
||||
void analyzePattern(const MatrixType& a) {
|
||||
eigen_assert(a.rows() == a.cols());
|
||||
Index size = a.cols();
|
||||
CholMatrixType tmp(size, size);
|
||||
ConstCholMatrixPtr pmat;
|
||||
ordering(a, pmat, tmp);
|
||||
analyzePattern_preordered(*pmat, doLDLT);
|
||||
ordering<NonHermitian>(a, pmat, tmp);
|
||||
analyzePattern_preordered(*pmat, DoLDLT);
|
||||
}
|
||||
void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
|
||||
|
||||
template <bool NonHermitian>
|
||||
void ordering(const MatrixType& a, ConstCholMatrixPtr& pmat, CholMatrixType& ap);
|
||||
|
||||
inline DiagonalScalar getDiag(Scalar x) { return internal::traits<Derived>::getDiag(x); }
|
||||
inline Scalar getSymm(Scalar x) { return internal::traits<Derived>::getSymm(x); }
|
||||
|
||||
/** keeps off-diagonal entries; drops diagonal entries */
|
||||
struct keep_diag {
|
||||
inline bool operator()(const Index& row, const Index& col, const Scalar&) const { return row != col; }
|
||||
@ -238,8 +244,8 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived> {
|
||||
PermutationMatrix<Dynamic, Dynamic, StorageIndex> m_P; // the permutation
|
||||
PermutationMatrix<Dynamic, Dynamic, StorageIndex> m_Pinv; // the inverse permutation
|
||||
|
||||
RealScalar m_shiftOffset;
|
||||
RealScalar m_shiftScale;
|
||||
DiagonalScalar m_shiftOffset;
|
||||
DiagonalScalar m_shiftScale;
|
||||
};
|
||||
|
||||
template <typename MatrixType_, int UpLo_ = Lower,
|
||||
@ -248,6 +254,12 @@ class SimplicialLLT;
|
||||
template <typename MatrixType_, int UpLo_ = Lower,
|
||||
typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> >
|
||||
class SimplicialLDLT;
|
||||
template <typename MatrixType_, int UpLo_ = Lower,
|
||||
typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> >
|
||||
class SimplicialNonHermitianLLT;
|
||||
template <typename MatrixType_, int UpLo_ = Lower,
|
||||
typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> >
|
||||
class SimplicialNonHermitianLDLT;
|
||||
template <typename MatrixType_, int UpLo_ = Lower,
|
||||
typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> >
|
||||
class SimplicialCholesky;
|
||||
@ -260,12 +272,15 @@ struct traits<SimplicialLLT<MatrixType_, UpLo_, Ordering_> > {
|
||||
typedef Ordering_ OrderingType;
|
||||
enum { UpLo = UpLo_ };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar DiagonalScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
|
||||
typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
|
||||
static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
|
||||
static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
|
||||
static inline DiagonalScalar getDiag(Scalar x) { return numext::real(x); }
|
||||
static inline Scalar getSymm(Scalar x) { return numext::conj(x); }
|
||||
};
|
||||
|
||||
template <typename MatrixType_, int UpLo_, typename Ordering_>
|
||||
@ -274,12 +289,49 @@ struct traits<SimplicialLDLT<MatrixType_, UpLo_, Ordering_> > {
|
||||
typedef Ordering_ OrderingType;
|
||||
enum { UpLo = UpLo_ };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar DiagonalScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
|
||||
typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
|
||||
static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
|
||||
static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
|
||||
static inline DiagonalScalar getDiag(Scalar x) { return numext::real(x); }
|
||||
static inline Scalar getSymm(Scalar x) { return numext::conj(x); }
|
||||
};
|
||||
|
||||
template <typename MatrixType_, int UpLo_, typename Ordering_>
|
||||
struct traits<SimplicialNonHermitianLLT<MatrixType_, UpLo_, Ordering_> > {
|
||||
typedef MatrixType_ MatrixType;
|
||||
typedef Ordering_ OrderingType;
|
||||
enum { UpLo = UpLo_ };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::Scalar DiagonalScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
|
||||
typedef TriangularView<const typename CholMatrixType::ConstTransposeReturnType, Eigen::Upper> MatrixU;
|
||||
static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
|
||||
static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.transpose()); }
|
||||
static inline DiagonalScalar getDiag(Scalar x) { return x; }
|
||||
static inline Scalar getSymm(Scalar x) { return x; }
|
||||
};
|
||||
|
||||
template <typename MatrixType_, int UpLo_, typename Ordering_>
|
||||
struct traits<SimplicialNonHermitianLDLT<MatrixType_, UpLo_, Ordering_> > {
|
||||
typedef MatrixType_ MatrixType;
|
||||
typedef Ordering_ OrderingType;
|
||||
enum { UpLo = UpLo_ };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::Scalar DiagonalScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
|
||||
typedef TriangularView<const typename CholMatrixType::ConstTransposeReturnType, Eigen::UnitUpper> MatrixU;
|
||||
static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
|
||||
static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.transpose()); }
|
||||
static inline DiagonalScalar getDiag(Scalar x) { return x; }
|
||||
static inline Scalar getSymm(Scalar x) { return x; }
|
||||
};
|
||||
|
||||
template <typename MatrixType_, int UpLo_, typename Ordering_>
|
||||
@ -287,6 +339,10 @@ struct traits<SimplicialCholesky<MatrixType_, UpLo_, Ordering_> > {
|
||||
typedef MatrixType_ MatrixType;
|
||||
typedef Ordering_ OrderingType;
|
||||
enum { UpLo = UpLo_ };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar DiagonalScalar;
|
||||
static inline DiagonalScalar getDiag(Scalar x) { return numext::real(x); }
|
||||
static inline Scalar getSymm(Scalar x) { return numext::conj(x); }
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
@ -346,7 +402,7 @@ class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<MatrixType_, U
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialLLT& compute(const MatrixType& matrix) {
|
||||
Base::template compute<false>(matrix);
|
||||
Base::template compute<false, false>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
@ -356,7 +412,7 @@ class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<MatrixType_, U
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a) { Base::analyzePattern(a, false); }
|
||||
void analyzePattern(const MatrixType& a) { Base::template analyzePattern<false, false>(a); }
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
@ -364,7 +420,7 @@ class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<MatrixType_, U
|
||||
*
|
||||
* \sa analyzePattern()
|
||||
*/
|
||||
void factorize(const MatrixType& a) { Base::template factorize<false>(a); }
|
||||
void factorize(const MatrixType& a) { Base::template factorize<false, false>(a); }
|
||||
|
||||
/** \returns the determinant of the underlying matrix from the current factorization */
|
||||
Scalar determinant() const {
|
||||
@ -434,7 +490,7 @@ class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<MatrixType_,
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialLDLT& compute(const MatrixType& matrix) {
|
||||
Base::template compute<true>(matrix);
|
||||
Base::template compute<true, false>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
@ -444,7 +500,7 @@ class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<MatrixType_,
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a) { Base::analyzePattern(a, true); }
|
||||
void analyzePattern(const MatrixType& a) { Base::template analyzePattern<true, false>(a); }
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
@ -452,7 +508,177 @@ class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<MatrixType_,
|
||||
*
|
||||
* \sa analyzePattern()
|
||||
*/
|
||||
void factorize(const MatrixType& a) { Base::template factorize<true>(a); }
|
||||
void factorize(const MatrixType& a) { Base::template factorize<true, false>(a); }
|
||||
|
||||
/** \returns the determinant of the underlying matrix from the current factorization */
|
||||
Scalar determinant() const { return Base::m_diag.prod(); }
|
||||
};
|
||||
|
||||
/** \ingroup SparseCholesky_Module
|
||||
* \class SimplicialNonHermitianLLT
|
||||
* \brief A direct sparse LLT Cholesky factorizations, for symmetric non-hermitian matrices.
|
||||
*
|
||||
* This class provides a LL^T Cholesky factorizations of sparse matrices that are
|
||||
* symmetric but not hermitian. For real matrices, this is equivalent to the regular LLT factorization.
|
||||
* The factorization allows for solving A.X = B where X and B can be either dense or sparse.
|
||||
*
|
||||
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
||||
* such that the factorized matrix is P A P^-1.
|
||||
*
|
||||
* \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
|
||||
* \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower
|
||||
* or Upper. Default is Lower.
|
||||
* \tparam Ordering_ The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
|
||||
*
|
||||
* \implsparsesolverconcept
|
||||
*
|
||||
* \sa class SimplicialNonHermitianLDLT, SimplicialLLT, class AMDOrdering, class NaturalOrdering
|
||||
*/
|
||||
template <typename MatrixType_, int UpLo_, typename Ordering_>
|
||||
class SimplicialNonHermitianLLT
|
||||
: public SimplicialCholeskyBase<SimplicialNonHermitianLLT<MatrixType_, UpLo_, Ordering_> > {
|
||||
public:
|
||||
typedef MatrixType_ MatrixType;
|
||||
enum { UpLo = UpLo_ };
|
||||
typedef SimplicialCholeskyBase<SimplicialNonHermitianLLT> Base;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef Matrix<Scalar, Dynamic, 1> VectorType;
|
||||
typedef internal::traits<SimplicialNonHermitianLLT> Traits;
|
||||
typedef typename Traits::MatrixL MatrixL;
|
||||
typedef typename Traits::MatrixU MatrixU;
|
||||
|
||||
public:
|
||||
/** Default constructor */
|
||||
SimplicialNonHermitianLLT() : Base() {}
|
||||
|
||||
/** Constructs and performs the LLT factorization of \a matrix */
|
||||
explicit SimplicialNonHermitianLLT(const MatrixType& matrix) : Base(matrix) {}
|
||||
|
||||
/** \returns an expression of the factor L */
|
||||
inline const MatrixL matrixL() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
|
||||
return Traits::getL(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** \returns an expression of the factor U (= L^*) */
|
||||
inline const MatrixU matrixU() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
|
||||
return Traits::getU(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialNonHermitianLLT& compute(const MatrixType& matrix) {
|
||||
Base::template compute<false, true>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
/** Performs a symbolic decomposition on the sparcity of \a matrix.
|
||||
*
|
||||
* This function is particularly useful when solving for several problems having the same structure.
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a) { Base::template analyzePattern<false, true>(a); }
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
* The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
||||
*
|
||||
* \sa analyzePattern()
|
||||
*/
|
||||
void factorize(const MatrixType& a) { Base::template factorize<false, true>(a); }
|
||||
|
||||
/** \returns the determinant of the underlying matrix from the current factorization */
|
||||
Scalar determinant() const {
|
||||
Scalar detL = Base::m_matrix.diagonal().prod();
|
||||
return detL * detL;
|
||||
}
|
||||
};
|
||||
|
||||
/** \ingroup SparseCholesky_Module
|
||||
* \class SimplicialNonHermitianLDLT
|
||||
* \brief A direct sparse LDLT Cholesky factorizations without square root, for symmetric non-hermitian matrices.
|
||||
*
|
||||
* This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
|
||||
* symmetric but not hermitian. For real matrices, this is equivalent to the regular LDLT factorization.
|
||||
* The factorization allows for solving A.X = B where X and B can be either dense or sparse.
|
||||
*
|
||||
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
||||
* such that the factorized matrix is P A P^-1.
|
||||
*
|
||||
* \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
|
||||
* \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower
|
||||
* or Upper. Default is Lower.
|
||||
* \tparam Ordering_ The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
|
||||
*
|
||||
* \implsparsesolverconcept
|
||||
*
|
||||
* \sa class SimplicialNonHermitianLLT, SimplicialLDLT, class AMDOrdering, class NaturalOrdering
|
||||
*/
|
||||
template <typename MatrixType_, int UpLo_, typename Ordering_>
|
||||
class SimplicialNonHermitianLDLT
|
||||
: public SimplicialCholeskyBase<SimplicialNonHermitianLDLT<MatrixType_, UpLo_, Ordering_> > {
|
||||
public:
|
||||
typedef MatrixType_ MatrixType;
|
||||
enum { UpLo = UpLo_ };
|
||||
typedef SimplicialCholeskyBase<SimplicialNonHermitianLDLT> Base;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef Matrix<Scalar, Dynamic, 1> VectorType;
|
||||
typedef internal::traits<SimplicialNonHermitianLDLT> Traits;
|
||||
typedef typename Traits::MatrixL MatrixL;
|
||||
typedef typename Traits::MatrixU MatrixU;
|
||||
|
||||
public:
|
||||
/** Default constructor */
|
||||
SimplicialNonHermitianLDLT() : Base() {}
|
||||
|
||||
/** Constructs and performs the LLT factorization of \a matrix */
|
||||
explicit SimplicialNonHermitianLDLT(const MatrixType& matrix) : Base(matrix) {}
|
||||
|
||||
/** \returns a vector expression of the diagonal D */
|
||||
inline const VectorType vectorD() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
||||
return Base::m_diag;
|
||||
}
|
||||
/** \returns an expression of the factor L */
|
||||
inline const MatrixL matrixL() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
||||
return Traits::getL(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** \returns an expression of the factor U (= L^*) */
|
||||
inline const MatrixU matrixU() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
||||
return Traits::getU(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialNonHermitianLDLT& compute(const MatrixType& matrix) {
|
||||
Base::template compute<true, true>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
/** Performs a symbolic decomposition on the sparcity of \a matrix.
|
||||
*
|
||||
* This function is particularly useful when solving for several problems having the same structure.
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a) { Base::template analyzePattern<true, true>(a); }
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
* The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
||||
*
|
||||
* \sa analyzePattern()
|
||||
*/
|
||||
void factorize(const MatrixType& a) { Base::template factorize<true, true>(a); }
|
||||
|
||||
/** \returns the determinant of the underlying matrix from the current factorization */
|
||||
Scalar determinant() const { return Base::m_diag.prod(); }
|
||||
@ -475,7 +701,6 @@ class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<Matr
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef Matrix<Scalar, Dynamic, 1> VectorType;
|
||||
typedef internal::traits<SimplicialCholesky> Traits;
|
||||
typedef internal::traits<SimplicialLDLT<MatrixType, UpLo> > LDLTTraits;
|
||||
typedef internal::traits<SimplicialLLT<MatrixType, UpLo> > LLTTraits;
|
||||
|
||||
@ -511,9 +736,9 @@ class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<Matr
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialCholesky& compute(const MatrixType& matrix) {
|
||||
if (m_LDLT)
|
||||
Base::template compute<true>(matrix);
|
||||
Base::template compute<true, false>(matrix);
|
||||
else
|
||||
Base::template compute<false>(matrix);
|
||||
Base::template compute<false, false>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
@ -523,7 +748,12 @@ class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<Matr
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a) { Base::analyzePattern(a, m_LDLT); }
|
||||
void analyzePattern(const MatrixType& a) {
|
||||
if (m_LDLT)
|
||||
Base::template analyzePattern<true, false>(a);
|
||||
else
|
||||
Base::template analyzePattern<false, false>(a);
|
||||
}
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
@ -533,9 +763,9 @@ class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<Matr
|
||||
*/
|
||||
void factorize(const MatrixType& a) {
|
||||
if (m_LDLT)
|
||||
Base::template factorize<true>(a);
|
||||
Base::template factorize<true, false>(a);
|
||||
else
|
||||
Base::template factorize<false>(a);
|
||||
Base::template factorize<false, false>(a);
|
||||
}
|
||||
|
||||
/** \internal */
|
||||
@ -594,6 +824,7 @@ class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<Matr
|
||||
};
|
||||
|
||||
template <typename Derived>
|
||||
template <bool NonHermitian>
|
||||
void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr& pmat, CholMatrixType& ap) {
|
||||
eigen_assert(a.rows() == a.cols());
|
||||
const Index size = a.rows();
|
||||
@ -602,7 +833,7 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMat
|
||||
if (!internal::is_same<OrderingType, NaturalOrdering<Index> >::value) {
|
||||
{
|
||||
CholMatrixType C;
|
||||
C = a.template selfadjointView<UpLo>();
|
||||
internal::permute_symm_to_fullsymm<UpLo, NonHermitian>(a, C, NULL);
|
||||
|
||||
OrderingType ordering;
|
||||
ordering(C, m_Pinv);
|
||||
@ -614,14 +845,14 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMat
|
||||
m_P.resize(0);
|
||||
|
||||
ap.resize(size, size);
|
||||
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
|
||||
internal::permute_symm_to_symm<UpLo, Upper, NonHermitian>(a, ap, m_P.indices().data());
|
||||
} else {
|
||||
m_Pinv.resize(0);
|
||||
m_P.resize(0);
|
||||
if (int(UpLo) == int(Lower) || MatrixType::IsRowMajor) {
|
||||
// we have to transpose the lower part to to the upper one
|
||||
ap.resize(size, size);
|
||||
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
|
||||
internal::permute_symm_to_symm<UpLo, Upper, NonHermitian>(a, ap, NULL);
|
||||
} else
|
||||
internal::simplicial_cholesky_grab_input<CholMatrixType, MatrixType>::run(a, pmat, ap);
|
||||
}
|
||||
|
||||
@ -67,7 +67,7 @@ void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrix
|
||||
}
|
||||
|
||||
template <typename Derived>
|
||||
template <bool DoLDLT>
|
||||
template <bool DoLDLT, bool NonHermitian>
|
||||
void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap) {
|
||||
using std::sqrt;
|
||||
|
||||
@ -97,7 +97,7 @@ void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType&
|
||||
for (typename CholMatrixType::InnerIterator it(ap, k); it; ++it) {
|
||||
StorageIndex i = it.index();
|
||||
if (i <= k) {
|
||||
y[i] += numext::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */
|
||||
y[i] += getSymm(it.value()); /* scatter A(i,k) into Y (sum duplicates) */
|
||||
Index len;
|
||||
for (len = 0; tags[i] != k; i = m_parent[i]) {
|
||||
pattern[len++] = i; /* L(k,i) is nonzero */
|
||||
@ -109,8 +109,8 @@ void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType&
|
||||
|
||||
/* compute numerical values kth row of L (a sparse triangular solve) */
|
||||
|
||||
RealScalar d =
|
||||
numext::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
|
||||
DiagonalScalar d =
|
||||
getDiag(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
|
||||
y[k] = Scalar(0);
|
||||
for (; top < size; ++top) {
|
||||
Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
|
||||
@ -120,14 +120,14 @@ void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType&
|
||||
/* the nonzero entry L(k,i) */
|
||||
Scalar l_ki;
|
||||
if (DoLDLT)
|
||||
l_ki = yi / numext::real(m_diag[i]);
|
||||
l_ki = yi / getDiag(m_diag[i]);
|
||||
else
|
||||
yi = l_ki = yi / Lx[Lp[i]];
|
||||
|
||||
Index p2 = Lp[i] + m_nonZerosPerCol[i];
|
||||
Index p;
|
||||
for (p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) y[Li[p]] -= numext::conj(Lx[p]) * yi;
|
||||
d -= numext::real(l_ki * numext::conj(yi));
|
||||
for (p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) y[Li[p]] -= getSymm(Lx[p]) * yi;
|
||||
d -= getDiag(l_ki * getSymm(yi));
|
||||
Li[p] = k; /* store L(k,i) in column form of L */
|
||||
Lx[p] = l_ki;
|
||||
++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
|
||||
@ -141,7 +141,7 @@ void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType&
|
||||
} else {
|
||||
Index p = Lp[k] + m_nonZerosPerCol[k]++;
|
||||
Li[p] = k; /* store L(k,k) = sqrt (d) in column k */
|
||||
if (d <= RealScalar(0)) {
|
||||
if (NonHermitian ? d == RealScalar(0) : numext::real(d) <= RealScalar(0)) {
|
||||
ok = false; /* failure, matrix is not positive definite */
|
||||
break;
|
||||
}
|
||||
|
||||
@ -34,13 +34,13 @@ namespace internal {
|
||||
template <typename MatrixType, unsigned int Mode>
|
||||
struct traits<SparseSelfAdjointView<MatrixType, Mode> > : traits<MatrixType> {};
|
||||
|
||||
template <int SrcMode, int DstMode, typename MatrixType, int DestOrder>
|
||||
template <int SrcMode, int DstMode, bool NonHermitian, typename MatrixType, int DestOrder>
|
||||
void permute_symm_to_symm(
|
||||
const MatrixType& mat,
|
||||
SparseMatrix<typename MatrixType::Scalar, DestOrder, typename MatrixType::StorageIndex>& _dest,
|
||||
const typename MatrixType::StorageIndex* perm = 0);
|
||||
|
||||
template <int Mode, typename MatrixType, int DestOrder>
|
||||
template <int Mode, bool NonHermitian, typename MatrixType, int DestOrder>
|
||||
void permute_symm_to_fullsymm(
|
||||
const MatrixType& mat,
|
||||
SparseMatrix<typename MatrixType::Scalar, DestOrder, typename MatrixType::StorageIndex>& _dest,
|
||||
@ -234,7 +234,7 @@ struct Assignment<DstXprType, SrcXprType, Functor, SparseSelfAdjoint2Sparse> {
|
||||
template <typename DestScalar, int StorageOrder>
|
||||
static void run(SparseMatrix<DestScalar, StorageOrder, StorageIndex>& dst, const SrcXprType& src,
|
||||
const AssignOpType& /*func*/) {
|
||||
internal::permute_symm_to_fullsymm<SrcXprType::Mode>(src.matrix(), dst);
|
||||
internal::permute_symm_to_fullsymm<SrcXprType::Mode, false>(src.matrix(), dst);
|
||||
}
|
||||
|
||||
// FIXME: the handling of += and -= in sparse matrices should be cleanup so that next two overloads could be reduced
|
||||
@ -405,7 +405,7 @@ struct product_evaluator<Product<Lhs, RhsView, DefaultProduct>, ProductTag, Spar
|
||||
***************************************************************************/
|
||||
namespace internal {
|
||||
|
||||
template <int Mode, typename MatrixType, int DestOrder>
|
||||
template <int Mode, bool NonHermitian, typename MatrixType, int DestOrder>
|
||||
void permute_symm_to_fullsymm(
|
||||
const MatrixType& mat,
|
||||
SparseMatrix<typename MatrixType::Scalar, DestOrder, typename MatrixType::StorageIndex>& _dest,
|
||||
@ -476,13 +476,13 @@ void permute_symm_to_fullsymm(
|
||||
dest.valuePtr()[k] = it.value();
|
||||
k = count[ip]++;
|
||||
dest.innerIndexPtr()[k] = jp;
|
||||
dest.valuePtr()[k] = numext::conj(it.value());
|
||||
dest.valuePtr()[k] = (NonHermitian ? it.value() : numext::conj(it.value()));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <int SrcMode_, int DstMode_, typename MatrixType, int DstOrder>
|
||||
template <int SrcMode_, int DstMode_, bool NonHermitian, typename MatrixType, int DstOrder>
|
||||
void permute_symm_to_symm(const MatrixType& mat,
|
||||
SparseMatrix<typename MatrixType::Scalar, DstOrder, typename MatrixType::StorageIndex>& _dest,
|
||||
const typename MatrixType::StorageIndex* perm) {
|
||||
@ -534,7 +534,7 @@ void permute_symm_to_symm(const MatrixType& mat,
|
||||
|
||||
if (!StorageOrderMatch) std::swap(ip, jp);
|
||||
if (((int(DstMode) == int(Lower) && ip < jp) || (int(DstMode) == int(Upper) && ip > jp)))
|
||||
dest.valuePtr()[k] = numext::conj(it.value());
|
||||
dest.valuePtr()[k] = (NonHermitian ? it.value() : numext::conj(it.value()));
|
||||
else
|
||||
dest.valuePtr()[k] = it.value();
|
||||
}
|
||||
@ -595,14 +595,14 @@ struct Assignment<DstXprType, SparseSymmetricPermutationProduct<MatrixType, Mode
|
||||
const internal::assign_op<Scalar, typename MatrixType::Scalar>&) {
|
||||
// internal::permute_symm_to_fullsymm<Mode>(m_matrix,_dest,m_perm.indices().data());
|
||||
SparseMatrix<Scalar, (Options & RowMajor) == RowMajor ? ColMajor : RowMajor, DstIndex> tmp;
|
||||
internal::permute_symm_to_fullsymm<Mode>(src.matrix(), tmp, src.perm().indices().data());
|
||||
internal::permute_symm_to_fullsymm<Mode, false>(src.matrix(), tmp, src.perm().indices().data());
|
||||
dst = tmp;
|
||||
}
|
||||
|
||||
template <typename DestType, unsigned int DestMode>
|
||||
static void run(SparseSelfAdjointView<DestType, DestMode>& dst, const SrcXprType& src,
|
||||
const internal::assign_op<Scalar, typename MatrixType::Scalar>&) {
|
||||
internal::permute_symm_to_symm<Mode, DestMode>(src.matrix(), dst.matrix(), src.perm().indices().data());
|
||||
internal::permute_symm_to_symm<Mode, DestMode, false>(src.matrix(), dst.matrix(), src.perm().indices().data());
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@ -20,6 +20,12 @@ void test_simplicial_cholesky_T() {
|
||||
SimplicialLDLT<SparseMatrixType, Upper> ldlt_colmajor_upper_amd;
|
||||
SimplicialLDLT<SparseMatrixType, Lower, NaturalOrdering<I_> > ldlt_colmajor_lower_nat;
|
||||
SimplicialLDLT<SparseMatrixType, Upper, NaturalOrdering<I_> > ldlt_colmajor_upper_nat;
|
||||
SimplicialNonHermitianLLT<SparseMatrixType, Lower> nhllt_colmajor_lower_amd;
|
||||
SimplicialNonHermitianLLT<SparseMatrixType, Upper> nhllt_colmajor_upper_amd;
|
||||
SimplicialNonHermitianLDLT<SparseMatrixType, Lower> nhldlt_colmajor_lower_amd;
|
||||
SimplicialNonHermitianLDLT<SparseMatrixType, Upper> nhldlt_colmajor_upper_amd;
|
||||
SimplicialNonHermitianLDLT<SparseMatrixType, Lower, NaturalOrdering<I_> > nhldlt_colmajor_lower_nat;
|
||||
SimplicialNonHermitianLDLT<SparseMatrixType, Upper, NaturalOrdering<I_> > nhldlt_colmajor_upper_nat;
|
||||
|
||||
check_sparse_spd_solving(chol_colmajor_lower_amd);
|
||||
check_sparse_spd_solving(chol_colmajor_upper_amd);
|
||||
@ -27,6 +33,10 @@ void test_simplicial_cholesky_T() {
|
||||
check_sparse_spd_solving(llt_colmajor_upper_amd);
|
||||
check_sparse_spd_solving(ldlt_colmajor_lower_amd);
|
||||
check_sparse_spd_solving(ldlt_colmajor_upper_amd);
|
||||
check_sparse_nonhermitian_solving(nhllt_colmajor_lower_amd);
|
||||
check_sparse_nonhermitian_solving(nhllt_colmajor_upper_amd);
|
||||
check_sparse_nonhermitian_solving(nhldlt_colmajor_lower_amd);
|
||||
check_sparse_nonhermitian_solving(nhldlt_colmajor_upper_amd);
|
||||
|
||||
check_sparse_spd_determinant(chol_colmajor_lower_amd);
|
||||
check_sparse_spd_determinant(chol_colmajor_upper_amd);
|
||||
@ -34,9 +44,15 @@ void test_simplicial_cholesky_T() {
|
||||
check_sparse_spd_determinant(llt_colmajor_upper_amd);
|
||||
check_sparse_spd_determinant(ldlt_colmajor_lower_amd);
|
||||
check_sparse_spd_determinant(ldlt_colmajor_upper_amd);
|
||||
check_sparse_nonhermitian_determinant(nhllt_colmajor_lower_amd);
|
||||
check_sparse_nonhermitian_determinant(nhllt_colmajor_upper_amd);
|
||||
check_sparse_nonhermitian_determinant(nhldlt_colmajor_lower_amd);
|
||||
check_sparse_nonhermitian_determinant(nhldlt_colmajor_upper_amd);
|
||||
|
||||
check_sparse_spd_solving(ldlt_colmajor_lower_nat, (std::min)(300, EIGEN_TEST_MAX_SIZE), 1000);
|
||||
check_sparse_spd_solving(ldlt_colmajor_upper_nat, (std::min)(300, EIGEN_TEST_MAX_SIZE), 1000);
|
||||
check_sparse_nonhermitian_solving(nhldlt_colmajor_lower_nat, (std::min)(300, EIGEN_TEST_MAX_SIZE), 1000);
|
||||
check_sparse_nonhermitian_solving(nhldlt_colmajor_upper_nat, (std::min)(300, EIGEN_TEST_MAX_SIZE), 1000);
|
||||
}
|
||||
|
||||
EIGEN_DECLARE_TEST(simplicial_cholesky) {
|
||||
|
||||
@ -484,6 +484,96 @@ void check_sparse_spd_determinant(Solver& solver) {
|
||||
}
|
||||
}
|
||||
|
||||
template <typename Solver, typename DenseMat>
|
||||
int generate_sparse_nonhermitian_problem(Solver&, typename Solver::MatrixType& A, typename Solver::MatrixType& halfA,
|
||||
DenseMat& dA, int maxSize = 300) {
|
||||
typedef typename Solver::MatrixType Mat;
|
||||
typedef typename Mat::Scalar Scalar;
|
||||
typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
|
||||
|
||||
int size = internal::random<int>(1, maxSize);
|
||||
double density = (std::max)(8. / static_cast<double>(size * size), 0.01);
|
||||
|
||||
Mat M(size, size);
|
||||
DenseMatrix dM(size, size);
|
||||
|
||||
initSparse<Scalar>(density, dM, M, ForceNonZeroDiag);
|
||||
|
||||
A = M * M.transpose();
|
||||
dA = dM * dM.transpose();
|
||||
|
||||
halfA.resize(size, size);
|
||||
if (Solver::UpLo == (Lower | Upper))
|
||||
halfA = A;
|
||||
else
|
||||
halfA = A.template triangularView<Solver::UpLo>();
|
||||
|
||||
return size;
|
||||
}
|
||||
|
||||
template <typename Solver>
|
||||
void check_sparse_nonhermitian_solving(Solver& solver, int maxSize = (std::min)(300, EIGEN_TEST_MAX_SIZE),
|
||||
int maxRealWorldSize = 100000) {
|
||||
typedef typename Solver::MatrixType Mat;
|
||||
typedef typename Mat::Scalar Scalar;
|
||||
typedef typename Mat::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> SpMat;
|
||||
typedef SparseVector<Scalar, 0, StorageIndex> SpVec;
|
||||
typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
|
||||
typedef Matrix<Scalar, Dynamic, 1> DenseVector;
|
||||
|
||||
// generate the problem
|
||||
Mat A, halfA;
|
||||
DenseMatrix dA;
|
||||
for (int i = 0; i < g_repeat; i++) {
|
||||
int size = generate_sparse_nonhermitian_problem(solver, A, halfA, dA, maxSize);
|
||||
|
||||
// generate the right hand sides
|
||||
int rhsCols = internal::random<int>(1, 16);
|
||||
double density = (std::max)(8. / static_cast<double>(size * rhsCols), 0.1);
|
||||
SpMat B(size, rhsCols);
|
||||
DenseVector b = DenseVector::Random(size);
|
||||
DenseMatrix dB(size, rhsCols);
|
||||
initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
|
||||
SpVec c = B.col(0);
|
||||
DenseVector dc = dB.col(0);
|
||||
|
||||
CALL_SUBTEST(check_sparse_solving(solver, A, b, dA, b));
|
||||
CALL_SUBTEST(check_sparse_solving(solver, halfA, b, dA, b));
|
||||
CALL_SUBTEST(check_sparse_solving(solver, A, dB, dA, dB));
|
||||
CALL_SUBTEST(check_sparse_solving(solver, halfA, dB, dA, dB));
|
||||
CALL_SUBTEST(check_sparse_solving(solver, A, B, dA, dB));
|
||||
CALL_SUBTEST(check_sparse_solving(solver, halfA, B, dA, dB));
|
||||
CALL_SUBTEST(check_sparse_solving(solver, A, c, dA, dc));
|
||||
CALL_SUBTEST(check_sparse_solving(solver, halfA, c, dA, dc));
|
||||
|
||||
// check only once
|
||||
if (i == 0) {
|
||||
b = DenseVector::Zero(size);
|
||||
check_sparse_solving(solver, A, b, dA, b);
|
||||
}
|
||||
}
|
||||
|
||||
EIGEN_UNUSED_VARIABLE(maxRealWorldSize);
|
||||
}
|
||||
|
||||
template <typename Solver>
|
||||
void check_sparse_nonhermitian_determinant(Solver& solver) {
|
||||
typedef typename Solver::MatrixType Mat;
|
||||
typedef typename Mat::Scalar Scalar;
|
||||
typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
|
||||
|
||||
// generate the problem
|
||||
Mat A, halfA;
|
||||
DenseMatrix dA;
|
||||
generate_sparse_nonhermitian_problem(solver, A, halfA, dA, 30);
|
||||
|
||||
for (int i = 0; i < g_repeat; i++) {
|
||||
check_sparse_determinant(solver, A, dA);
|
||||
check_sparse_determinant(solver, halfA, dA);
|
||||
}
|
||||
}
|
||||
|
||||
template <typename Solver>
|
||||
void check_sparse_zero_matrix(Solver& solver) {
|
||||
typedef typename Solver::MatrixType Mat;
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user