diff --git a/Eigen/Sparse b/Eigen/Sparse index df0310f89..0e5da9c3a 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -16,9 +16,16 @@ #ifdef EIGEN_TAUCS_SUPPORT + // taucs.h declares a lot of mess + #define isnan + #define finite + #define isinf extern "C" { #include "taucs.h" } + #undef isnan + #undef finite + #undef isinf #ifdef min #undef min diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index d1d10158a..b6c13350c 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -195,11 +195,12 @@ template template void SparseLLT::solveInPlace(MatrixBase &b) const { + if (m_status & MatrixLIsDirty) + matrixL(); + const int size = m_matrix.rows(); ei_assert(size==b.rows()); - if (m_status & MatrixLIsDirty) - matrixL(); Base::solveInPlace(b); } diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h index 7578a12c7..b7d4f5bbd 100644 --- a/Eigen/src/Sparse/SparseLLT.h +++ b/Eigen/src/Sparse/SparseLLT.h @@ -139,7 +139,6 @@ void SparseLLT::compute(const MatrixType& a) for (int j = 0; j < size; ++j) { Scalar x = ei_real(a.coeff(j,j)); - int endSize = size-j-1; // TODO better estimate of the density ! tempVector.init(density>0.001? IsDense : IsSparse); @@ -191,7 +190,8 @@ bool SparseLLT::solveInPlace(MatrixBase &b) const ei_assert(size==b.rows()); m_matrix.solveTriangularInPlace(b); - m_matrix.adjoint().solveTriangularInPlace(b); + // FIXME should be .adjoint() but it fails to compile... + m_matrix.transpose().solveTriangularInPlace(b); return true; } diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h index 5edf0fc27..22fef1ceb 100644 --- a/Eigen/src/Sparse/TaucsSupport.h +++ b/Eigen/src/Sparse/TaucsSupport.h @@ -151,7 +151,7 @@ void SparseLLT::compute(const MatrixType& a) else { // use the faster Multifrontal routine - m_taucsSupernodalFactor = taucs_ccs_factor_llt_ll(&taucsMatA); + m_taucsSupernodalFactor = taucs_ccs_factor_llt_mf(&taucsMatA); } m_status = (m_status & ~IncompleteFactorization) | CompleteFactorization | MatrixLIsDirty; } @@ -177,16 +177,19 @@ template template void SparseLLT::solveInPlace(MatrixBase &b) const { - const int size = m_matrix.rows(); - ei_assert(size==b.rows()); - if (m_status & MatrixLIsDirty) { -// ei_assert(!(m_status & SupernodalFactorIsDirty)); -// taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b); - //matrixL(); + // TODO use taucs's supernodal solver, in particular check types, storage order, etc. + // VectorXb x(b.rows()); + // for (int j=0; j" # On other platform use ctest as usual # -MACRO(EI_ADD_TEST testname) +macro(ei_add_test testname) - SET(targetname test_${testname}) + set(targetname test_${testname}) - SET(filename ${testname}.cpp) - ADD_EXECUTABLE(${targetname} ${filename}) + set(filename ${testname}.cpp) + add_executable(${targetname} ${filename}) - IF(NOT EIGEN_NO_ASSERTION_CHECKING) + if(NOT EIGEN_NO_ASSERTION_CHECKING) - SET_TARGET_PROPERTIES(${targetname} PROPERTIES COMPILE_FLAGS "-fexceptions") - OPTION(EIGEN_DEBUG_ASSERTS "Enable debuging of assertions" OFF) - IF(EIGEN_DEBUG_ASSERTS) - SET_TARGET_PROPERTIES(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_DEBUG_ASSERTS=1") - ENDIF(EIGEN_DEBUG_ASSERTS) + set_target_properties(${targetname} PROPERTIES COMPILE_FLAGS "-fexceptions") + option(EIGEN_DEBUG_ASSERTS "Enable debuging of assertions" OFF) + if(EIGEN_DEBUG_ASSERTS) + set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_DEBUG_ASSERTS=1") + endif(EIGEN_DEBUG_ASSERTS) - ELSE(NOT EIGEN_NO_ASSERTION_CHECKING) + else(NOT EIGEN_NO_ASSERTION_CHECKING) - SET_TARGET_PROPERTIES(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_NO_ASSERTION_CHECKING=1") + set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_NO_ASSERTION_CHECKING=1") - ENDIF(NOT EIGEN_NO_ASSERTION_CHECKING) + endif(NOT EIGEN_NO_ASSERTION_CHECKING) - IF(${ARGC} GREATER 1) - EI_ADD_TARGET_PROPERTY(${targetname} COMPILE_FLAGS "${ARGV1}") - ENDIF(${ARGC} GREATER 1) + if(${ARGC} GREATER 1) + ei_add_target_property(${targetname} COMPILE_FLAGS "${ARGV1}") + endif(${ARGC} GREATER 1) - EI_ADD_TARGET_PROPERTY(${targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}") + ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}") - IF(TEST_LIB) + if(TEST_LIB) target_link_libraries(${targetname} Eigen2) - ENDIF(TEST_LIB) + endif(TEST_LIB) - if(GSL_FOUND) - target_link_libraries(${targetname} ${GSL_LIBRARIES}) - endif(GSL_FOUND) + target_link_libraries(${targetname} ${EXTERNAL_LIBS}) - IF(WIN32) - ADD_TEST(${testname} "${targetname}") - ELSE(WIN32) - ADD_TEST(${testname} "${CMAKE_CURRENT_SOURCE_DIR}/runtest.sh" "${testname}") - ENDIF(WIN32) + if(WIN32) + add_test(${testname} "${targetname}") + else(WIN32) + add_test(${testname} "${CMAKE_CURRENT_SOURCE_DIR}/runtest.sh" "${testname}") + endif(WIN32) -ENDMACRO(EI_ADD_TEST) +endmacro(ei_add_test) -ENABLE_TESTING() +enable_testing() -IF(TEST_LIB) - ADD_DEFINITIONS("-DEIGEN_EXTERN_INSTANTIATIONS=1") -ENDIF(TEST_LIB) +if(TEST_LIB) + add_definitions("-DEIGEN_EXTERN_INSTANTIATIONS=1") +endif(TEST_LIB) -EI_ADD_TEST(meta) -EI_ADD_TEST(sizeof) -EI_ADD_TEST(dynalloc) -EI_ADD_TEST(nomalloc) -EI_ADD_TEST(packetmath) -EI_ADD_TEST(basicstuff) -EI_ADD_TEST(linearstructure) -EI_ADD_TEST(cwiseop) -EI_ADD_TEST(sum) -EI_ADD_TEST(product_small) -EI_ADD_TEST(product_large ${EI_OFLAG}) -EI_ADD_TEST(adjoint) -EI_ADD_TEST(submatrices) -EI_ADD_TEST(miscmatrices) -EI_ADD_TEST(commainitializer) -EI_ADD_TEST(smallvectors) -EI_ADD_TEST(map) -EI_ADD_TEST(array) -EI_ADD_TEST(triangular) -EI_ADD_TEST(cholesky) -EI_ADD_TEST(lu ${EI_OFLAG}) -EI_ADD_TEST(determinant) -EI_ADD_TEST(inverse) -EI_ADD_TEST(qr) -EI_ADD_TEST(eigensolver) -EI_ADD_TEST(svd) -EI_ADD_TEST(geometry) -EI_ADD_TEST(hyperplane) -EI_ADD_TEST(parametrizedline) -EI_ADD_TEST(regression) -EI_ADD_TEST(sparse ${EI_OFLAG}) +ei_add_test(meta) +ei_add_test(sizeof) +ei_add_test(dynalloc) +ei_add_test(nomalloc) +ei_add_test(packetmath) +ei_add_test(basicstuff) +ei_add_test(linearstructure) +ei_add_test(cwiseop) +ei_add_test(sum) +ei_add_test(product_small) +ei_add_test(product_large ${EI_OFLAG}) +ei_add_test(adjoint) +ei_add_test(submatrices) +ei_add_test(miscmatrices) +ei_add_test(commainitializer) +ei_add_test(smallvectors) +ei_add_test(map) +ei_add_test(array) +ei_add_test(triangular) +ei_add_test(cholesky) +ei_add_test(lu ${EI_OFLAG}) +ei_add_test(determinant) +ei_add_test(inverse) +ei_add_test(qr) +ei_add_test(eigensolver) +ei_add_test(svd) +ei_add_test(geometry) +ei_add_test(hyperplane) +ei_add_test(parametrizedline) +ei_add_test(regression) +ei_add_test(sparse ) -ENDIF(BUILD_TESTS) +endif(BUILD_TESTS) diff --git a/test/sparse.cpp b/test/sparse.cpp index f54835972..ca80e6362 100644 --- a/test/sparse.cpp +++ b/test/sparse.cpp @@ -23,6 +23,8 @@ // Eigen. If not, see . #include "main.h" +#include +#include #include enum { @@ -46,8 +48,7 @@ initSparse(double density, { Scalar v = (ei_random(0,1) < density) ? ei_random() : 0; if ((flags&ForceNonZeroDiag) && (i==j)) - while (ei_abs(v)<1e-2) - v = ei_random(); + v = ei_random(Scalar(5.),Scalar(20.)); if ((flags & MakeLowerTriangular) && j>i) v = 0; else if ((flags & MakeUpperTriangular) && j void sparse(int rows, int cols) refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); VERIFY_IS_APPROX(m, refMat); - + #if 0 // test InnerIterators and Block expressions for(int j=0; j void sparse(int rows, int cols) // TODO test row major } + #endif // test LLT { + SparseMatrix m2(rows, cols); + DenseMatrix refMat2(rows, cols); + + DenseVector b = DenseVector::Random(cols); + DenseVector refX(cols), x(cols); + + initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); + refMat2 += refMat2.adjoint(); + refMat2.diagonal() *= 0.5; + + refMat2.llt().solve(b, &refX); +// std::cerr << refMat2 << "\n\n" << refMat2.llt().matrixL() << "\n\n"; +// std::cerr << m2 << "\n\n"; + typedef SparseMatrix SparseSelfAdjointMatrix; + x = b; + SparseLLT (m2).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision()) && "LLT: default"); + #ifdef EIGEN_CHOLMOD_SUPPORT + x = b; + SparseLLT(m2).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision()) && "LLT: cholmod"); + #endif + #ifdef EIGEN_TAUCS_SUPPORT + x = b; + SparseLLT(m2,IncompleteFactorization).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision()) && "LLT: taucs (IncompleteFactorization)"); + x = b; + SparseLLT(m2,SupernodalMultifrontal).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision()) && "LLT: taucs (SupernodalMultifrontal)"); + x = b; + SparseLLT(m2,SupernodalLeftLooking).solveInPlace(x); + VERIFY(refX.isApprox(x,test_precision()) && "LLT: taucs (SupernodalLeftLooking)"); + #endif } }