From 10c6bcdc2eac68d0b48ecfe4d391079610eb6124 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Mon, 19 Dec 2016 14:07:42 +0100 Subject: [PATCH] Add support for long indexes and for (real-valued) row-major matrices to CholmodSupport module --- Eigen/src/CholmodSupport/CholmodSupport.h | 71 ++++++++++++++++++----- test/cholmod_support.cpp | 42 +++++++++----- 2 files changed, 84 insertions(+), 29 deletions(-) diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h index 571972023..a07efb1d5 100644 --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h @@ -32,7 +32,7 @@ template<> struct cholmod_configure_matrix > { } }; -// Other scalar types are not yet suppotred by Cholmod +// Other scalar types are not yet supported by Cholmod // template<> struct cholmod_configure_matrix { // template // static void run(CholmodType& mat) { @@ -124,6 +124,9 @@ cholmod_sparse viewAsCholmod(const SparseSelfAdjointView::IsComplex == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); + if(_Options & RowMajorBit) res.stype *=-1; return res; } @@ -159,6 +162,44 @@ MappedSparseMatrix viewAsEigen(cholmod_sparse& cm) static_cast(cm.p), static_cast(cm.i),static_cast(cm.x) ); } +namespace internal { + +// template specializations for int and long that call the correct cholmod method + +#define EIGEN_CHOLMOD_SPECIALIZE0(ret, name) \ + template ret cm_ ## name (cholmod_common &Common) { return cholmod_ ## name (&Common); } \ + template<> ret cm_ ## name (cholmod_common &Common) { return cholmod_l_ ## name (&Common); } + +#define EIGEN_CHOLMOD_SPECIALIZE1(ret, name, t1, a1) \ + template ret cm_ ## name (t1& a1, cholmod_common &Common) { return cholmod_ ## name (&a1, &Common); } \ + template<> ret cm_ ## name (t1& a1, cholmod_common &Common) { return cholmod_l_ ## name (&a1, &Common); } + +EIGEN_CHOLMOD_SPECIALIZE0(int, start) +EIGEN_CHOLMOD_SPECIALIZE0(int, finish) + +EIGEN_CHOLMOD_SPECIALIZE1(int, free_factor, cholmod_factor*, L) +EIGEN_CHOLMOD_SPECIALIZE1(int, free_dense, cholmod_dense*, X) +EIGEN_CHOLMOD_SPECIALIZE1(int, free_sparse, cholmod_sparse*, A) + +EIGEN_CHOLMOD_SPECIALIZE1(cholmod_factor*, analyze, cholmod_sparse, A) + +template cholmod_dense* cm_solve (int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) { return cholmod_solve (sys, &L, &B, &Common); } +template<> cholmod_dense* cm_solve (int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) { return cholmod_l_solve (sys, &L, &B, &Common); } + +template cholmod_sparse* cm_spsolve (int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) { return cholmod_spsolve (sys, &L, &B, &Common); } +template<> cholmod_sparse* cm_spsolve (int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) { return cholmod_l_spsolve (sys, &L, &B, &Common); } + +template +int cm_factorize_p (cholmod_sparse* A, double beta[2], _StorageIndex* fset, size_t fsize, cholmod_factor* L, cholmod_common &Common) { return cholmod_factorize_p (A, beta, fset, fsize, L, &Common); } +template<> +int cm_factorize_p (cholmod_sparse* A, double beta[2], long* fset, size_t fsize, cholmod_factor* L, cholmod_common &Common) { return cholmod_l_factorize_p (A, beta, fset, fsize, L, &Common); } + +#undef EIGEN_CHOLMOD_SPECIALIZE0 +#undef EIGEN_CHOLMOD_SPECIALIZE1 + +} // namespace internal + + enum CholmodMode { CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt }; @@ -195,7 +236,7 @@ class CholmodBase : public SparseSolverBase { EIGEN_STATIC_ASSERT((internal::is_same::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); m_shiftOffset[0] = m_shiftOffset[1] = 0.0; - cholmod_start(&m_cholmod); + internal::cm_start(m_cholmod); } explicit CholmodBase(const MatrixType& matrix) @@ -203,15 +244,15 @@ class CholmodBase : public SparseSolverBase { EIGEN_STATIC_ASSERT((internal::is_same::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); m_shiftOffset[0] = m_shiftOffset[1] = 0.0; - cholmod_start(&m_cholmod); + internal::cm_start(m_cholmod); compute(matrix); } ~CholmodBase() { if(m_cholmodFactor) - cholmod_free_factor(&m_cholmodFactor, &m_cholmod); - cholmod_finish(&m_cholmod); + internal::cm_free_factor(m_cholmodFactor, m_cholmod); + internal::cm_finish(m_cholmod); } inline StorageIndex cols() const { return internal::convert_index(m_cholmodFactor->n); } @@ -219,7 +260,7 @@ class CholmodBase : public SparseSolverBase /** \brief Reports whether previous computation was successful. * - * \returns \c Success if computation was succesful, + * \returns \c Success if computation was successful, * \c NumericalIssue if the matrix.appears to be negative. */ ComputationInfo info() const @@ -246,11 +287,11 @@ class CholmodBase : public SparseSolverBase { if(m_cholmodFactor) { - cholmod_free_factor(&m_cholmodFactor, &m_cholmod); + internal::cm_free_factor(m_cholmodFactor, m_cholmod); m_cholmodFactor = 0; } cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView()); - m_cholmodFactor = cholmod_analyze(&A, &m_cholmod); + m_cholmodFactor = internal::cm_analyze(A, m_cholmod); this->m_isInitialized = true; this->m_info = Success; @@ -268,7 +309,7 @@ class CholmodBase : public SparseSolverBase { eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView()); - cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod); + internal::cm_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, m_cholmod); // If the factorization failed, minor is the column at which it did. On success minor == n. this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue); @@ -289,19 +330,20 @@ class CholmodBase : public SparseSolverBase EIGEN_UNUSED_VARIABLE(size); eigen_assert(size==b.rows()); - // Cholmod needs column-major stoarge without inner-stride, which corresponds to the default behavior of Ref. + // Cholmod needs column-major storage without inner-stride, which corresponds to the default behavior of Ref. Ref > b_ref(b.derived()); cholmod_dense b_cd = viewAsCholmod(b_ref); - cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod); + cholmod_dense* x_cd = internal::cm_solve(CHOLMOD_A, *m_cholmodFactor, b_cd, m_cholmod); if(!x_cd) { this->m_info = NumericalIssue; return; } // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) + // NOTE Actually, the copy can be avoided by calling cholmod_solve2 instead of cholmod_solve dest = Matrix::Map(reinterpret_cast(x_cd->x),b.rows(),b.cols()); - cholmod_free_dense(&x_cd, &m_cholmod); + internal::cm_free_dense(x_cd, m_cholmod); } /** \internal */ @@ -316,15 +358,16 @@ class CholmodBase : public SparseSolverBase // note: cs stands for Cholmod Sparse Ref > b_ref(b.const_cast_derived()); cholmod_sparse b_cs = viewAsCholmod(b_ref); - cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod); + cholmod_sparse* x_cs = internal::cm_spsolve(CHOLMOD_A, *m_cholmodFactor, b_cs, m_cholmod); if(!x_cs) { this->m_info = NumericalIssue; return; } // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) + // NOTE cholmod_spsolve in fact just calls the dense solver for blocks of 4 columns at a time (similar to Eigen's sparse solver) dest.derived() = viewAsEigen(*x_cs); - cholmod_free_sparse(&x_cs, &m_cholmod); + internal::cm_free_sparse(x_cs, m_cholmod); } #endif // EIGEN_PARSED_BY_DOXYGEN diff --git a/test/cholmod_support.cpp b/test/cholmod_support.cpp index a7eda28f7..931207334 100644 --- a/test/cholmod_support.cpp +++ b/test/cholmod_support.cpp @@ -12,21 +12,21 @@ #include -template void test_cholmod_T() +template void test_cholmod_ST() { - CholmodDecomposition, Lower> g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt); - CholmodDecomposition, Upper> g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt); - CholmodDecomposition, Lower> g_llt_colmajor_lower; g_llt_colmajor_lower.setMode(CholmodSimplicialLLt); - CholmodDecomposition, Upper> g_llt_colmajor_upper; g_llt_colmajor_upper.setMode(CholmodSimplicialLLt); - CholmodDecomposition, Lower> g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt); - CholmodDecomposition, Upper> g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt); + CholmodDecomposition g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt); + CholmodDecomposition g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt); + CholmodDecomposition g_llt_colmajor_lower; g_llt_colmajor_lower.setMode(CholmodSimplicialLLt); + CholmodDecomposition g_llt_colmajor_upper; g_llt_colmajor_upper.setMode(CholmodSimplicialLLt); + CholmodDecomposition g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt); + CholmodDecomposition g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt); - CholmodSupernodalLLT, Lower> chol_colmajor_lower; - CholmodSupernodalLLT, Upper> chol_colmajor_upper; - CholmodSimplicialLLT, Lower> llt_colmajor_lower; - CholmodSimplicialLLT, Upper> llt_colmajor_upper; - CholmodSimplicialLDLT, Lower> ldlt_colmajor_lower; - CholmodSimplicialLDLT, Upper> ldlt_colmajor_upper; + CholmodSupernodalLLT chol_colmajor_lower; + CholmodSupernodalLLT chol_colmajor_upper; + CholmodSimplicialLLT llt_colmajor_lower; + CholmodSimplicialLLT llt_colmajor_upper; + CholmodSimplicialLDLT ldlt_colmajor_lower; + CholmodSimplicialLDLT ldlt_colmajor_upper; check_sparse_spd_solving(g_chol_colmajor_lower); check_sparse_spd_solving(g_chol_colmajor_upper); @@ -50,8 +50,20 @@ template void test_cholmod_T() check_sparse_spd_determinant(ldlt_colmajor_upper); } +template void test_cholmod_T() +{ + test_cholmod_ST >(); +} + void test_cholmod_support() { - CALL_SUBTEST_1(test_cholmod_T()); - CALL_SUBTEST_2(test_cholmod_T >()); + CALL_SUBTEST_11( (test_cholmod_T()) ); + CALL_SUBTEST_12( (test_cholmod_T()) ); + CALL_SUBTEST_13( (test_cholmod_T()) ); + CALL_SUBTEST_14( (test_cholmod_T()) ); + CALL_SUBTEST_21( (test_cholmod_T, ColMajor, int >()) ); + CALL_SUBTEST_22( (test_cholmod_T, ColMajor, long>()) ); + // TODO complex row-major matrices do not work at the moment: + // CALL_SUBTEST_23( (test_cholmod_T, RowMajor, int >()) ); + // CALL_SUBTEST_24( (test_cholmod_T, RowMajor, long>()) ); }