diff options
| -rw-r--r-- | CMakeLists.txt | 7 | ||||
| -rw-r--r-- | Eigen/src/PaStiXSupport/PaStiXSupport.h | 16 | ||||
| -rw-r--r-- | Eigen/src/UmfPackSupport/UmfPackSupport.h | 8 | ||||
| -rw-r--r-- | test/CMakeLists.txt | 9 | ||||
| -rw-r--r-- | test/bicgstab.cpp | 6 | ||||
| -rw-r--r-- | test/cholmod_support.cpp | 18 | ||||
| -rw-r--r-- | test/conjugate_gradient.cpp | 6 | ||||
| -rw-r--r-- | test/pardiso_support.cpp | 10 | ||||
| -rw-r--r-- | test/pastix_support.cpp | 10 | ||||
| -rw-r--r-- | test/simplicial_cholesky.cpp | 6 | ||||
| -rw-r--r-- | test/sparse.h | 1 | ||||
| -rw-r--r-- | test/sparse_solver.h | 123 | ||||
| -rw-r--r-- | test/superlu_support.cpp | 14 | ||||
| -rw-r--r-- | test/umfpack_support.cpp | 14 |
14 files changed, 177 insertions, 71 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index d39976aee..bc77bcb2e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -338,6 +338,13 @@ if(EIGEN_BUILD_BTL) add_subdirectory(bench/btl EXCLUDE_FROM_ALL) endif(EIGEN_BUILD_BTL) +if(TEST_REAL_CASES) + if(NOT WIN32) + add_subdirectory(bench/spbench EXCLUDE_FROM_ALL) + set(ENV(EIGEN_MATRIX_DIR) ${TEST_REAL_CASES}) + endif(NOT WIN32) +endif(TEST_REAL_CASES) + ei_testing_print_summary() message(STATUS "") diff --git a/Eigen/src/PaStiXSupport/PaStiXSupport.h b/Eigen/src/PaStiXSupport/PaStiXSupport.h index 57ce0f106..7989a92ca 100644 --- a/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -136,10 +136,10 @@ namespace internal // WARNING It is assumed here that successive calls to this routine are done // with matrices having the same pattern. template <typename MatrixType> - void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans) + void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans, bool& hasTranspose) { eigen_assert(In.cols()==In.rows() && " Can only symmetrize the graph of a square matrix"); - if (StrMatTrans.outerSize() == 0) + if (!hasTranspose) { //First call to this routine, need to compute the structural pattern of In^T StrMatTrans = In.transpose(); // Set the elements of the matrix to zero @@ -148,6 +148,7 @@ namespace internal for (typename MatrixType::InnerIterator it(StrMatTrans, i); it; ++it) it.valueRef() = 0.0; } + hasTranspose = true; } Out = (StrMatTrans + In).eval(); } @@ -172,6 +173,7 @@ class PastixBase PastixBase():m_initisOk(false),m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) { m_pastixdata = 0; + m_hasTranspose = false; PastixInit(); } @@ -314,7 +316,6 @@ class PastixBase m_iparm(IPARM_END_TASK) = API_TASK_CLEAN; internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(), m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data()); - } Derived& compute (MatrixType& mat); @@ -328,6 +329,7 @@ class PastixBase mutable SparseMatrix<Scalar, ColMajor> m_mat_null; // An input null matrix mutable Matrix<Scalar, Dynamic,1> m_vec_null; // An input null vector mutable SparseMatrix<Scalar, ColMajor> m_StrMatTrans; // The transpose pattern of the input matrix + mutable bool m_hasTranspose; // The transpose of the current matrix has already been computed mutable int m_comm; // The MPI communicator identifier mutable Matrix<Index,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters @@ -357,6 +359,7 @@ void PastixBase<Derived>::PastixInit() PastixDestroy(); m_pastixdata = 0; m_iparm(IPARM_MODIFY_PARAMETER) = API_YES; + m_hasTranspose = false; } m_iparm(IPARM_START_TASK) = API_TASK_INIT; @@ -537,7 +540,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > temp = matrix; else { - internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans); + internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans, m_hasTranspose); } m_iparm[IPARM_SYM] = API_SYM_NO; m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; @@ -559,7 +562,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > temp = matrix; else { - internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans); + internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose); } m_iparm(IPARM_SYM) = API_SYM_NO; @@ -581,7 +584,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > temp = matrix; else { - internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans); + internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose); } m_iparm(IPARM_SYM) = API_SYM_NO; m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; @@ -591,6 +594,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > using Base::m_iparm; using Base::m_dparm; using Base::m_StrMatTrans; + using Base::m_hasTranspose; }; /** \ingroup PaStiXSupport_Module diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index 5921a86b0..636bba87d 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -124,9 +124,11 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N * \brief A sparse LU factorization and solver based on UmfPack * * This class allows to solve for A.X = B sparse linear problems via a LU factorization - * using the UmfPack library. The sparse matrix A must be column-major, squared and full rank. + * using the UmfPack library. The sparse matrix A must be in a compressed column-major form, squared and full rank. * The vectors or matrices X and B can be either dense or sparse. * + * WARNING The Eigen column-major SparseMatrix is not always in compressed form. + * The user should call makeCompressed() to get a matrix in CSC suitable for UMFPACK * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * * \sa \ref TutorialSparseDirectSolvers @@ -198,7 +200,9 @@ class UmfPackLU return m_q; } - /** Computes the sparse Cholesky decomposition of \a matrix */ + /** Computes the sparse Cholesky decomposition of \a matrix + * Note that the matrix should be in compressed format. Please, use makeCompressed() to get it !! + */ void compute(const MatrixType& matrix) { analyzePattern(matrix); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9ed11f601..b42958b9b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -18,6 +18,14 @@ set(LAPACK_FOUND TRUE) set(BLAS_LIBRARIES eigen_blas) set(LAPACK_LIBRARIES eigen_lapack) +if(TEST_REAL_CASES) + if(NOT WIN32) + add_definitions( -DEIGEN_MATRIXDIR="${TEST_REAL_CASES}" ) + else(NOT WIN32) + message(STATUS, "REAL CASES CAN NOT BE CURRENTLY TESTED ON WIN32") + endif(NOT WIN32) +endif(TEST_REAL_CASES) + set(SPARSE_LIBS " ") find_package(Cholmod) @@ -192,6 +200,7 @@ ei_add_test(vectorwiseop) ei_add_test(simplicial_cholesky) ei_add_test(conjugate_gradient) ei_add_test(bicgstab) +ei_add_test(gmres) if(UMFPACK_FOUND) diff --git a/test/bicgstab.cpp b/test/bicgstab.cpp index bb9da898f..2b6403583 100644 --- a/test/bicgstab.cpp +++ b/test/bicgstab.cpp @@ -40,8 +40,6 @@ template<typename T> void test_bicgstab_T() void test_bicgstab() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1(test_bicgstab_T<double>()); - CALL_SUBTEST_2(test_bicgstab_T<std::complex<double> >()); - } + CALL_SUBTEST_1(test_bicgstab_T<double>()); + CALL_SUBTEST_2(test_bicgstab_T<std::complex<double> >()); } diff --git a/test/cholmod_support.cpp b/test/cholmod_support.cpp index a5de86f60..9351558a9 100644 --- a/test/cholmod_support.cpp +++ b/test/cholmod_support.cpp @@ -42,18 +42,16 @@ template<typename T> void test_cholmod_T() check_sparse_spd_solving(ldlt_colmajor_lower); check_sparse_spd_solving(ldlt_colmajor_upper); -// check_sparse_spd_determinant(chol_colmajor_lower); -// check_sparse_spd_determinant(chol_colmajor_upper); -// check_sparse_spd_determinant(llt_colmajor_lower); -// check_sparse_spd_determinant(llt_colmajor_upper); -// check_sparse_spd_determinant(ldlt_colmajor_lower); -// check_sparse_spd_determinant(ldlt_colmajor_upper); + check_sparse_spd_determinant(chol_colmajor_lower); + check_sparse_spd_determinant(chol_colmajor_upper); + check_sparse_spd_determinant(llt_colmajor_lower); + check_sparse_spd_determinant(llt_colmajor_upper); + check_sparse_spd_determinant(ldlt_colmajor_lower); + check_sparse_spd_determinant(ldlt_colmajor_upper); } void test_cholmod_support() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1(test_cholmod_T<double>()); - CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >()); - } + CALL_SUBTEST_1(test_cholmod_T<double>()); + CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >()); } diff --git a/test/conjugate_gradient.cpp b/test/conjugate_gradient.cpp index e76327bac..f24f35817 100644 --- a/test/conjugate_gradient.cpp +++ b/test/conjugate_gradient.cpp @@ -40,8 +40,6 @@ template<typename T> void test_conjugate_gradient_T() void test_conjugate_gradient() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1(test_conjugate_gradient_T<double>()); - CALL_SUBTEST_2(test_conjugate_gradient_T<std::complex<double> >()); - } + CALL_SUBTEST_1(test_conjugate_gradient_T<double>()); + CALL_SUBTEST_2(test_conjugate_gradient_T<std::complex<double> >()); } diff --git a/test/pardiso_support.cpp b/test/pardiso_support.cpp index 11cb98e10..67efad6d8 100644 --- a/test/pardiso_support.cpp +++ b/test/pardiso_support.cpp @@ -22,10 +22,8 @@ template<typename T> void test_pardiso_T() void test_pardiso_support() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1(test_pardiso_T<float>()); - CALL_SUBTEST_2(test_pardiso_T<double>()); - CALL_SUBTEST_3(test_pardiso_T< std::complex<float> >()); - CALL_SUBTEST_4(test_pardiso_T< std::complex<double> >()); - } + CALL_SUBTEST_1(test_pardiso_T<float>()); + CALL_SUBTEST_2(test_pardiso_T<double>()); + CALL_SUBTEST_3(test_pardiso_T< std::complex<float> >()); + CALL_SUBTEST_4(test_pardiso_T< std::complex<double> >()); } diff --git a/test/pastix_support.cpp b/test/pastix_support.cpp index 80c5bbf27..dbce30d1c 100644 --- a/test/pastix_support.cpp +++ b/test/pastix_support.cpp @@ -52,10 +52,8 @@ template<typename T> void test_pastix_T_LU() void test_pastix_support() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1(test_pastix_T<float>()); - CALL_SUBTEST_2(test_pastix_T<double>()); - CALL_SUBTEST_3( (test_pastix_T_LU<std::complex<float> >()) ); - CALL_SUBTEST_4(test_pastix_T_LU<std::complex<double> >()); - } + CALL_SUBTEST_1(test_pastix_T<float>()); + CALL_SUBTEST_2(test_pastix_T<double>()); + CALL_SUBTEST_3( (test_pastix_T_LU<std::complex<float> >()) ); + CALL_SUBTEST_4(test_pastix_T_LU<std::complex<double> >()); }
\ No newline at end of file diff --git a/test/simplicial_cholesky.cpp b/test/simplicial_cholesky.cpp index 8ebc91a02..f1af0e467 100644 --- a/test/simplicial_cholesky.cpp +++ b/test/simplicial_cholesky.cpp @@ -50,8 +50,6 @@ template<typename T> void test_simplicial_cholesky_T() void test_simplicial_cholesky() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1(test_simplicial_cholesky_T<double>()); - CALL_SUBTEST_2(test_simplicial_cholesky_T<std::complex<double> >()); - } + CALL_SUBTEST_1(test_simplicial_cholesky_T<double>()); + CALL_SUBTEST_2(test_simplicial_cholesky_T<std::complex<double> >()); } diff --git a/test/sparse.h b/test/sparse.h index 1ce1d2ac4..84cd695f2 100644 --- a/test/sparse.h +++ b/test/sparse.h @@ -193,4 +193,5 @@ initSparse(double density, } } +#include <unsupported/Eigen/SparseExtra> #endif // EIGEN_TESTSPARSE_H diff --git a/test/sparse_solver.h b/test/sparse_solver.h index d6c111362..62c0e9495 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -74,6 +74,56 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, VERIFY(x.isApprox(refX,test_precision<Scalar>())); } +template<typename Scalar> +inline std::string get_matrixfolder() +{ + std::string mat_folder = EIGEN_MATRIXDIR; + if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value ) + mat_folder = mat_folder + static_cast<string>("/complex/"); + else + mat_folder = mat_folder + static_cast<string>("/real/"); + return mat_folder; +} + +template<typename Solver, typename Rhs> +void check_sparse_solving_real_cases(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const Rhs& refX) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef typename Mat::RealScalar RealScalar; + + Rhs x(b.rows(), b.cols()); + + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: factorization failed (check_sparse_solving_real_cases)\n"; + exit(0); + return; + } + x = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: solving failed\n"; + return; + } + + RealScalar res_error; + // Compute the norm of the relative error + if(refX.size() != 0) + res_error = (refX - x).norm()/refX.norm(); + else + { + // Compute the relative residual norm + res_error = (b - A * x).norm()/b.norm(); + } + if (res_error > test_precision<Scalar>() ){ + std::cerr << "Test " << g_test_stack.back() << " failed in "EI_PP_MAKE_STRING(__FILE__) + << " (" << EI_PP_MAKE_STRING(__LINE__) << ")" << std::endl << std::endl; + abort(); + } + +} template<typename Solver, typename DenseMat> void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) { @@ -121,6 +171,7 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; + typedef typename Mat::Index Index; typedef SparseMatrix<Scalar,ColMajor> SpMat; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; @@ -137,13 +188,37 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver) DenseVector b = DenseVector::Random(size); DenseMatrix dB(size,rhsCols); initSparse<Scalar>(density, dB, B); + + for (int i = 0; i < g_repeat; i++) { + check_sparse_solving(solver, A, b, dA, b); + check_sparse_solving(solver, halfA, b, dA, b); + check_sparse_solving(solver, A, dB, dA, dB); + check_sparse_solving(solver, halfA, dB, dA, dB); + check_sparse_solving(solver, A, B, dA, dB); + check_sparse_solving(solver, halfA, B, dA, dB); + } - check_sparse_solving(solver, A, b, dA, b); - check_sparse_solving(solver, halfA, b, dA, b); - check_sparse_solving(solver, A, dB, dA, dB); - check_sparse_solving(solver, halfA, dB, dA, dB); - check_sparse_solving(solver, A, B, dA, dB); - check_sparse_solving(solver, halfA, B, dA, dB); + // First, get the folder +#ifdef EIGEN_MATRIXDIR + if (internal::is_same<Scalar, float>::value + || internal::is_same<Scalar, std::complex<float> >::value) + return ; + + std::string mat_folder = get_matrixfolder<Scalar>(); + MatrixMarketIterator<Scalar> it(mat_folder); + for (; it; ++it) + { + if (it.sym() == SPD){ + Mat halfA; + PermutationMatrix<Dynamic, Dynamic, Index> pnull; + halfA.template selfadjointView<Solver::UpLo>() = it.matrix().template triangularView<Eigen::Lower>().twistedBy(pnull); + + std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; + check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); + check_sparse_solving_real_cases(solver, halfA, it.rhs(), it.refX()); + } + } +#endif } template<typename Solver> void check_sparse_spd_determinant(Solver& solver) @@ -156,9 +231,11 @@ template<typename Solver> void check_sparse_spd_determinant(Solver& solver) Mat A, halfA; DenseMatrix dA; generate_sparse_spd_problem(solver, A, halfA, dA, 30); - - check_sparse_determinant(solver, A, dA); - check_sparse_determinant(solver, halfA, dA ); + + for (int i = 0; i < g_repeat; i++) { + check_sparse_determinant(solver, A, dA); + check_sparse_determinant(solver, halfA, dA ); + } } template<typename Solver, typename DenseMat> @@ -194,9 +271,27 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver) DenseVector b = DenseVector::Random(size); DenseMatrix dB = DenseMatrix::Random(size,rhsCols); + A.makeCompressed(); + for (int i = 0; i < g_repeat; i++) { + check_sparse_solving(solver, A, b, dA, b); + check_sparse_solving(solver, A, dB, dA, dB); + } + + // First, get the folder +#ifdef EIGEN_MATRIXDIR + if (internal::is_same<Scalar, float>::value + || internal::is_same<Scalar, std::complex<float> >::value) + return ; + + std::string mat_folder = get_matrixfolder<Scalar>(); + MatrixMarketIterator<Scalar> it(mat_folder); + for (; it; ++it) + { + std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; + check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); + } +#endif - check_sparse_solving(solver, A, b, dA, b); - check_sparse_solving(solver, A, dB, dA, dB); } template<typename Solver> void check_sparse_square_determinant(Solver& solver) @@ -209,6 +304,8 @@ template<typename Solver> void check_sparse_square_determinant(Solver& solver) Mat A; DenseMatrix dA; generate_sparse_square_problem(solver, A, dA, 30); - - check_sparse_determinant(solver, A, dA); + A.makeCompressed(); + for (int i = 0; i < g_repeat; i++) { + check_sparse_determinant(solver, A, dA); + } } diff --git a/test/superlu_support.cpp b/test/superlu_support.cpp index dbabaf2e1..ad435943b 100644 --- a/test/superlu_support.cpp +++ b/test/superlu_support.cpp @@ -28,12 +28,10 @@ void test_superlu_support() { - for(int i = 0; i < g_repeat; i++) { - SuperLU<SparseMatrix<double> > superlu_double_colmajor; - SuperLU<SparseMatrix<std::complex<double> > > superlu_cplxdouble_colmajor; - CALL_SUBTEST_1( check_sparse_square_solving(superlu_double_colmajor) ); - CALL_SUBTEST_2( check_sparse_square_solving(superlu_cplxdouble_colmajor) ); - CALL_SUBTEST_1( check_sparse_square_determinant(superlu_double_colmajor) ); - CALL_SUBTEST_2( check_sparse_square_determinant(superlu_cplxdouble_colmajor) ); - } + SuperLU<SparseMatrix<double> > superlu_double_colmajor; + SuperLU<SparseMatrix<std::complex<double> > > superlu_cplxdouble_colmajor; + CALL_SUBTEST_1( check_sparse_square_solving(superlu_double_colmajor) ); + CALL_SUBTEST_2( check_sparse_square_solving(superlu_cplxdouble_colmajor) ); + CALL_SUBTEST_1( check_sparse_square_determinant(superlu_double_colmajor) ); + CALL_SUBTEST_2( check_sparse_square_determinant(superlu_cplxdouble_colmajor) ); } diff --git a/test/umfpack_support.cpp b/test/umfpack_support.cpp index 6e459a922..f900e92c6 100644 --- a/test/umfpack_support.cpp +++ b/test/umfpack_support.cpp @@ -28,12 +28,10 @@ void test_umfpack_support() { - for(int i = 0; i < g_repeat; i++) { - UmfPackLU<SparseMatrix<double> > umfpack_double_colmajor; - UmfPackLU<SparseMatrix<std::complex<double> > > umfpack_cplxdouble_colmajor; - CALL_SUBTEST_1(check_sparse_square_solving(umfpack_double_colmajor)); - CALL_SUBTEST_2(check_sparse_square_solving(umfpack_cplxdouble_colmajor)); - CALL_SUBTEST_1(check_sparse_square_determinant(umfpack_double_colmajor)); - CALL_SUBTEST_2(check_sparse_square_determinant(umfpack_cplxdouble_colmajor)); - } + UmfPackLU<SparseMatrix<double, ColMajor> > umfpack_double_colmajor; + UmfPackLU<SparseMatrix<std::complex<double> > > umfpack_cplxdouble_colmajor; + CALL_SUBTEST_1(check_sparse_square_solving(umfpack_double_colmajor)); + CALL_SUBTEST_2(check_sparse_square_solving(umfpack_cplxdouble_colmajor)); + CALL_SUBTEST_1(check_sparse_square_determinant(umfpack_double_colmajor)); + CALL_SUBTEST_2(check_sparse_square_determinant(umfpack_cplxdouble_colmajor)); } |