diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index 423287b5b..f3ce975cf 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -58,6 +58,7 @@ class SimplicialCholeskyBase : public SparseSolverBase { enum { UpLo = internal::traits::UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; + typedef typename internal::traits::DiagonalScalar DiagonalScalar; typedef typename MatrixType::StorageIndex StorageIndex; typedef SparseMatrix CholMatrixType; typedef CholMatrixType const* ConstCholMatrixPtr; @@ -114,7 +115,7 @@ class SimplicialCholeskyBase : public SparseSolverBase { * * \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 { protected: /** Computes the sparse Cholesky decomposition of \a matrix */ - template + template 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(matrix, pmat, tmp); analyzePattern_preordered(*pmat, DoLDLT); - factorize_preordered(*pmat); + factorize_preordered(*pmat); } - template + template void factorize(const MatrixType& a) { eigen_assert(a.rows() == a.cols()); Index size = a.cols(); @@ -200,28 +201,33 @@ class SimplicialCholeskyBase : public SparseSolverBase { // If there is no ordering, try to directly use the input matrix without any copy internal::simplicial_cholesky_grab_input::run(a, pmat, tmp); } else { - tmp.template selfadjointView() = a.template selfadjointView().twistedBy(m_P); + internal::permute_symm_to_symm(a, tmp, m_P.indices().data()); pmat = &tmp; } - factorize_preordered(*pmat); + factorize_preordered(*pmat); } - template + template void factorize_preordered(const CholMatrixType& a); - void analyzePattern(const MatrixType& a, bool doLDLT) { + template + 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(a, pmat, tmp); + analyzePattern_preordered(*pmat, DoLDLT); } void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); + template void ordering(const MatrixType& a, ConstCholMatrixPtr& pmat, CholMatrixType& ap); + inline DiagonalScalar getDiag(Scalar x) { return internal::traits::getDiag(x); } + inline Scalar getSymm(Scalar x) { return internal::traits::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 { PermutationMatrix m_P; // the permutation PermutationMatrix m_Pinv; // the inverse permutation - RealScalar m_shiftOffset; - RealScalar m_shiftScale; + DiagonalScalar m_shiftOffset; + DiagonalScalar m_shiftScale; }; template > class SimplicialLDLT; +template > +class SimplicialNonHermitianLLT; +template > +class SimplicialNonHermitianLDLT; template > class SimplicialCholesky; @@ -260,12 +272,15 @@ struct traits > { typedef Ordering_ OrderingType; enum { UpLo = UpLo_ }; typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar DiagonalScalar; typedef typename MatrixType::StorageIndex StorageIndex; typedef SparseMatrix CholMatrixType; typedef TriangularView MatrixL; typedef TriangularView 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 @@ -274,12 +289,49 @@ struct traits > { typedef Ordering_ OrderingType; enum { UpLo = UpLo_ }; typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar DiagonalScalar; typedef typename MatrixType::StorageIndex StorageIndex; typedef SparseMatrix CholMatrixType; typedef TriangularView MatrixL; typedef TriangularView 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 +struct traits > { + 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 CholMatrixType; + typedef TriangularView MatrixL; + typedef TriangularView 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 +struct traits > { + 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 CholMatrixType; + typedef TriangularView MatrixL; + typedef TriangularView 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 @@ -287,6 +339,10 @@ struct traits > { 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(matrix); + Base::template compute(matrix); return *this; } @@ -356,7 +412,7 @@ class SimplicialLLT : public SimplicialCholeskyBase(a); } /** Performs a numeric decomposition of \a matrix * @@ -364,7 +420,7 @@ class SimplicialLLT : public SimplicialCholeskyBase(a); } + void factorize(const MatrixType& a) { Base::template factorize(a); } /** \returns the determinant of the underlying matrix from the current factorization */ Scalar determinant() const { @@ -434,7 +490,7 @@ class SimplicialLDLT : public SimplicialCholeskyBase(matrix); + Base::template compute(matrix); return *this; } @@ -444,7 +500,7 @@ class SimplicialLDLT : public SimplicialCholeskyBase(a); } /** Performs a numeric decomposition of \a matrix * @@ -452,7 +508,177 @@ class SimplicialLDLT : public SimplicialCholeskyBase(a); } + void factorize(const MatrixType& a) { Base::template factorize(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 +class SimplicialNonHermitianLLT + : public SimplicialCholeskyBase > { + public: + typedef MatrixType_ MatrixType; + enum { UpLo = UpLo_ }; + typedef SimplicialCholeskyBase Base; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::StorageIndex StorageIndex; + typedef SparseMatrix CholMatrixType; + typedef Matrix VectorType; + typedef internal::traits 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(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(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(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 +class SimplicialNonHermitianLDLT + : public SimplicialCholeskyBase > { + public: + typedef MatrixType_ MatrixType; + enum { UpLo = UpLo_ }; + typedef SimplicialCholeskyBase Base; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::StorageIndex StorageIndex; + typedef SparseMatrix CholMatrixType; + typedef Matrix VectorType; + typedef internal::traits 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(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(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(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 CholMatrixType; typedef Matrix VectorType; - typedef internal::traits Traits; typedef internal::traits > LDLTTraits; typedef internal::traits > LLTTraits; @@ -511,9 +736,9 @@ class SimplicialCholesky : public SimplicialCholeskyBase(matrix); + Base::template compute(matrix); else - Base::template compute(matrix); + Base::template compute(matrix); return *this; } @@ -523,7 +748,12 @@ class SimplicialCholesky : public SimplicialCholeskyBase(a); + else + Base::template analyzePattern(a); + } /** Performs a numeric decomposition of \a matrix * @@ -533,9 +763,9 @@ class SimplicialCholesky : public SimplicialCholeskyBase(a); + Base::template factorize(a); else - Base::template factorize(a); + Base::template factorize(a); } /** \internal */ @@ -594,6 +824,7 @@ class SimplicialCholesky : public SimplicialCholeskyBase +template void SimplicialCholeskyBase::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::ordering(const MatrixType& a, ConstCholMat if (!internal::is_same >::value) { { CholMatrixType C; - C = a.template selfadjointView(); + internal::permute_symm_to_fullsymm(a, C, NULL); OrderingType ordering; ordering(C, m_Pinv); @@ -614,14 +845,14 @@ void SimplicialCholeskyBase::ordering(const MatrixType& a, ConstCholMat m_P.resize(0); ap.resize(size, size); - ap.template selfadjointView() = a.template selfadjointView().twistedBy(m_P); + internal::permute_symm_to_symm(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() = a.template selfadjointView(); + internal::permute_symm_to_symm(a, ap, NULL); } else internal::simplicial_cholesky_grab_input::run(a, pmat, ap); } diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h index abfbbe6bc..0b13c56b5 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h @@ -67,7 +67,7 @@ void SimplicialCholeskyBase::analyzePattern_preordered(const CholMatrix } template -template +template void SimplicialCholeskyBase::factorize_preordered(const CholMatrixType& ap) { using std::sqrt; @@ -97,7 +97,7 @@ void SimplicialCholeskyBase::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::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::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::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; } diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 2f087c9c4..3402baeb6 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -34,13 +34,13 @@ namespace internal { template struct traits > : traits {}; -template +template void permute_symm_to_symm( const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::StorageIndex* perm = 0); -template +template void permute_symm_to_fullsymm( const MatrixType& mat, SparseMatrix& _dest, @@ -234,7 +234,7 @@ struct Assignment { template static void run(SparseMatrix& dst, const SrcXprType& src, const AssignOpType& /*func*/) { - internal::permute_symm_to_fullsymm(src.matrix(), dst); + internal::permute_symm_to_fullsymm(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, ProductTag, Spar ***************************************************************************/ namespace internal { -template +template void permute_symm_to_fullsymm( const MatrixType& mat, SparseMatrix& _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 +template void permute_symm_to_symm(const MatrixType& mat, SparseMatrix& _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&) { // internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); SparseMatrix tmp; - internal::permute_symm_to_fullsymm(src.matrix(), tmp, src.perm().indices().data()); + internal::permute_symm_to_fullsymm(src.matrix(), tmp, src.perm().indices().data()); dst = tmp; } template static void run(SparseSelfAdjointView& dst, const SrcXprType& src, const internal::assign_op&) { - internal::permute_symm_to_symm(src.matrix(), dst.matrix(), src.perm().indices().data()); + internal::permute_symm_to_symm(src.matrix(), dst.matrix(), src.perm().indices().data()); } }; diff --git a/test/simplicial_cholesky.cpp b/test/simplicial_cholesky.cpp index ca6749671..ed9321871 100644 --- a/test/simplicial_cholesky.cpp +++ b/test/simplicial_cholesky.cpp @@ -20,6 +20,12 @@ void test_simplicial_cholesky_T() { SimplicialLDLT ldlt_colmajor_upper_amd; SimplicialLDLT > ldlt_colmajor_lower_nat; SimplicialLDLT > ldlt_colmajor_upper_nat; + SimplicialNonHermitianLLT nhllt_colmajor_lower_amd; + SimplicialNonHermitianLLT nhllt_colmajor_upper_amd; + SimplicialNonHermitianLDLT nhldlt_colmajor_lower_amd; + SimplicialNonHermitianLDLT nhldlt_colmajor_upper_amd; + SimplicialNonHermitianLDLT > nhldlt_colmajor_lower_nat; + SimplicialNonHermitianLDLT > 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) { diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 033df835f..50cb463de 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -484,6 +484,96 @@ void check_sparse_spd_determinant(Solver& solver) { } } +template +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 DenseMatrix; + + int size = internal::random(1, maxSize); + double density = (std::max)(8. / static_cast(size * size), 0.01); + + Mat M(size, size); + DenseMatrix dM(size, size); + + initSparse(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(); + + return size; +} + +template +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 SpMat; + typedef SparseVector SpVec; + typedef Matrix DenseMatrix; + typedef Matrix 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(1, 16); + double density = (std::max)(8. / static_cast(size * rhsCols), 0.1); + SpMat B(size, rhsCols); + DenseVector b = DenseVector::Random(size); + DenseMatrix dB(size, rhsCols); + initSparse(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 +void check_sparse_nonhermitian_determinant(Solver& solver) { + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix 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 void check_sparse_zero_matrix(Solver& solver) { typedef typename Solver::MatrixType Mat;