modify the unit tests of sparse linear solvers to enable tests on real matrices, from MatrixMarket for instance

11 files changed, 154 insertions, 63 deletions

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));
 }