Add SimplicialNonHermitianLLT and SimplicialNonHermitianLDLT
diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index 423287b..f3ce975 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h
@@ -58,6 +58,7 @@ 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 @@ * * \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 @@ 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 @@ // 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 @@ 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, @@ -250,6 +256,12 @@ 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; namespace internal { @@ -260,12 +272,15 @@ 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 @@ 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 @@ 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 @@ /** 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 @@ * * \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 @@ * * \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 @@ /** 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 @@ * * \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 @@ * * \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 @@ 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 @@ /** 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 @@ * * \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 @@ */ 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 @@ }; 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 @@ 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 @@ 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); }
diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h index abfbbe6..0b13c56 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
@@ -67,7 +67,7 @@ } 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 @@ 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 @@ /* 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 @@ /* 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 @@ } 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 2f087c9..3402bae 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h
@@ -34,13 +34,13 @@ 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 @@ 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 @@ ***************************************************************************/ 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 @@ 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 @@ 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 @@ 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()); } };
diff --git a/test/simplicial_cholesky.cpp b/test/simplicial_cholesky.cpp index ca67496..ed93218 100644 --- a/test/simplicial_cholesky.cpp +++ b/test/simplicial_cholesky.cpp
@@ -20,6 +20,12 @@ 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 @@ 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 @@ 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 033df83..50cb463 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h
@@ -484,6 +484,96 @@ } } +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;