diff options
-rw-r--r-- | Eigen/src/Sparse/CholmodSupport.h | 9 | ||||
-rw-r--r-- | Eigen/src/Sparse/TaucsSupport.h | 12 | ||||
-rw-r--r-- | test/sparse_solvers.cpp | 34 |
3 files changed, 34 insertions, 21 deletions
diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index a9ef2d3a4..cf407240f 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -31,22 +31,22 @@ void ei_cholmod_configure_matrix(CholmodType& mat) if (ei_is_same_type<Scalar,float>::ret) { mat.xtype = CHOLMOD_REAL; - mat.dtype = 1; + mat.dtype = CHOLMOD_SINGLE; } else if (ei_is_same_type<Scalar,double>::ret) { mat.xtype = CHOLMOD_REAL; - mat.dtype = 0; + mat.dtype = CHOLMOD_DOUBLE; } else if (ei_is_same_type<Scalar,std::complex<float> >::ret) { mat.xtype = CHOLMOD_COMPLEX; - mat.dtype = 1; + mat.dtype = CHOLMOD_SINGLE; } else if (ei_is_same_type<Scalar,std::complex<double> >::ret) { mat.xtype = CHOLMOD_COMPLEX; - mat.dtype = 0; + mat.dtype = CHOLMOD_DOUBLE; } else { @@ -74,6 +74,7 @@ cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix() ei_cholmod_configure_matrix<Scalar>(res); + if (Derived::Flags & SelfAdjoint) { if (Derived::Flags & Upper) diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h index 0caa8cbed..2a1963f5b 100644 --- a/Eigen/src/Sparse/TaucsSupport.h +++ b/Eigen/src/Sparse/TaucsSupport.h @@ -50,6 +50,7 @@ taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix() ei_assert(false && "Scalar type not supported by TAUCS"); } + // FIXME 1) shapes are not in the Flags and 2) it seems Taucs ignores these flags anyway and only accept lower symmetric matrices if (Flags & Upper) res.flags |= TAUCS_UPPER; if (Flags & Lower) @@ -86,6 +87,7 @@ class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType> using Base::m_flags; using Base::m_matrix; using Base::m_status; + using Base::m_succeeded; public: @@ -126,10 +128,16 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a) m_taucsSupernodalFactor = 0; } + taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix(); + if (m_flags & IncompleteFactorization) { - taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix(); taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0); + if(!taucsRes) + { + m_succeeded = false; + return; + } // the matrix returned by Taucs is not necessarily sorted, // so let's copy it in two steps DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsRes); @@ -141,7 +149,6 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a) } else { - taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix(); if ( (m_flags & SupernodalLeftLooking) || ((!(m_flags & SupernodalMultifrontal)) && (m_flags & MemoryEfficient)) ) { @@ -154,6 +161,7 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a) } m_status = (m_status & ~IncompleteFactorization) | CompleteFactorization | MatrixLIsDirty; } + m_succeeded = true; } template<typename MatrixType> diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index fab2ab56e..00df1bffd 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -131,16 +131,20 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) #endif #ifdef EIGEN_TAUCS_SUPPORT - x = b; - SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,IncompleteFactorization).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)"); // TODO fix TAUCS with complexes - x = b; - SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)"); - x = b; - SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)"); + if (!NumTraits<Scalar>::IsComplex) + { + x = b; +// SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,IncompleteFactorization).solveInPlace(x); +// VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)"); + + x = b; + SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)"); + x = b; + SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)"); + } #endif } @@ -154,9 +158,9 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) DenseVector b = DenseVector::Random(cols); DenseVector refX(cols), x(cols); - //initSPD(density, refMat2, m2); +// initSPD(density, refMat2, m2); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0); - refMat2 += refMat2.adjoint(); + refMat2 += (refMat2.adjoint()).eval(); refMat2.diagonal() *= 0.5; refX = refMat2.llt().solve(b); // FIXME use LLT to compute the reference because LDLT seems to fail with large matrices @@ -202,7 +206,7 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) } if (slu.solve(b, &x, SvAdjoint)) { -// VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>())); + VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>())); } if (count==0) { @@ -236,8 +240,8 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) void test_sparse_solvers() { for(int i = 0; i < g_repeat; i++) { -// CALL_SUBTEST(sparse_solvers<double>(8, 8) ); - CALL_SUBTEST(sparse_solvers<std::complex<double> >(16, 16) ); -// CALL_SUBTEST(sparse_solvers<double>(100, 100) ); + CALL_SUBTEST_1(sparse_solvers<double>(8, 8) ); + CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(16, 16) ); + CALL_SUBTEST_1(sparse_solvers<double>(100, 100) ); } } |