From f804a319c81cb1629abb9bdc97dd74a2d2dec3d7 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Thu, 29 Mar 2012 14:32:54 +0200 Subject: modify the unit tests of sparse linear solvers to enable tests on real matrices, from MatrixMarket for instance --- test/sparse_solver.h | 123 +++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 110 insertions(+), 13 deletions(-) (limited to 'test/sparse_solver.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())); } +template +inline std::string get_matrixfolder() +{ + std::string mat_folder = EIGEN_MATRIXDIR; + if( internal::is_same >::value || internal::is_same >::value ) + mat_folder = mat_folder + static_cast("/complex/"); + else + mat_folder = mat_folder + static_cast("/real/"); + return mat_folder; +} + +template +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() ){ + 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 void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) { @@ -121,6 +171,7 @@ template void check_sparse_spd_solving(Solver& solver) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; + typedef typename Mat::Index Index; typedef SparseMatrix SpMat; typedef Matrix DenseMatrix; typedef Matrix DenseVector; @@ -137,13 +188,37 @@ template void check_sparse_spd_solving(Solver& solver) DenseVector b = DenseVector::Random(size); DenseMatrix dB(size,rhsCols); initSparse(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::value + || internal::is_same >::value) + return ; + + std::string mat_folder = get_matrixfolder(); + MatrixMarketIterator it(mat_folder); + for (; it; ++it) + { + if (it.sym() == SPD){ + Mat halfA; + PermutationMatrix pnull; + halfA.template selfadjointView() = it.matrix().template triangularView().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 void check_sparse_spd_determinant(Solver& solver) @@ -156,9 +231,11 @@ template 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 @@ -194,9 +271,27 @@ template 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::value + || internal::is_same >::value) + return ; + + std::string mat_folder = get_matrixfolder(); + MatrixMarketIterator 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 void check_sparse_square_determinant(Solver& solver) @@ -209,6 +304,8 @@ template 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); + } } -- cgit v1.2.3