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;
