refactor sparse solving unit tests

8 files changed, 232 insertions, 210 deletions

diff --git a/unsupported/test/conjugate_gradient.cpp b/unsupported/test/conjugate_gradient.cpp
index 47cafaac3..91a205a77 100644
--- a/unsupported/test/conjugate_gradient.cpp
+++ b/unsupported/test/conjugate_gradient.cpp
@@ -22,7 +22,7 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#include "sparse_llt.h"
+#include "sparse_solver.h"
 #include <Eigen/IterativeSolvers>
 
 template<typename T> void test_conjugate_gradient_T()
@@ -32,10 +32,10 @@ template<typename T> void test_conjugate_gradient_T()
   ConjugateGradient<SparseMatrix<T>, Lower, IdentityPreconditioner> cg_colmajor_lower_I;
   ConjugateGradient<SparseMatrix<T>, Upper, IdentityPreconditioner> cg_colmajor_upper_I;
 
-  sparse_llt(cg_colmajor_lower_diag);
-  sparse_llt(cg_colmajor_upper_diag);
-  sparse_llt(cg_colmajor_lower_I);
-  sparse_llt(cg_colmajor_upper_I);
+  CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_diag)  );
+  CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_diag)  );
+  CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_I)     );
+  CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_I)     );
 }
 
 void test_conjugate_gradient()
diff --git a/unsupported/test/simplicial_cholesky.cpp b/unsupported/test/simplicial_cholesky.cpp
index 00a5e1096..b9a55ef36 100644
--- a/unsupported/test/simplicial_cholesky.cpp
+++ b/unsupported/test/simplicial_cholesky.cpp
@@ -22,7 +22,7 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#include "sparse_llt.h"
+#include "sparse_solver.h"
 
 template<typename T> void test_simplicial_cholesky_T()
 {
@@ -33,12 +33,19 @@ template<typename T> void test_simplicial_cholesky_T()
   SimplicialLDLt<SparseMatrix<T>, Lower> ldlt_colmajor_lower;
   SimplicialLDLt<SparseMatrix<T>, Upper> ldlt_colmajor_upper;
 
-  sparse_llt(chol_colmajor_lower);
-  sparse_llt(chol_colmajor_upper);
-  sparse_llt(llt_colmajor_lower);
-  sparse_llt(llt_colmajor_upper);
-  sparse_llt(ldlt_colmajor_lower);
-  sparse_llt(ldlt_colmajor_upper);
+  check_sparse_spd_solving(chol_colmajor_lower);
+  check_sparse_spd_solving(chol_colmajor_upper);
+  check_sparse_spd_solving(llt_colmajor_lower);
+  check_sparse_spd_solving(llt_colmajor_upper);
+  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);
 }
 
 void test_simplicial_cholesky()
diff --git a/unsupported/test/sparse_ldlt.cpp b/unsupported/test/sparse_ldlt.cpp
index c387f78ac..4cb61bf6c 100644
--- a/unsupported/test/sparse_ldlt.cpp
+++ b/unsupported/test/sparse_ldlt.cpp
@@ -132,17 +132,17 @@ template<typename Scalar,typename Index> void sparse_ldlt(int rows, int cols)
     
     
     // with a sparse rhs
-//     SparseMatrixType spB(rows,cols), spX(rows,cols);
-//     B.diagonal().array() += 1;
-//     spB = B.sparseView(0.5,1);
-//     
-//     ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB));
-// 
-//     spX = SimplicialCholesky<SparseMatrixType, Lower>(m3).solve(spB);
-//     VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
+    SparseMatrixType spB(rows,cols), spX(rows,cols);
+    B.diagonal().array() += 1;
+    spB = B.sparseView(0.5,1);
+    
+    ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB));
+
+    spX = SimplicialCholesky<SparseMatrixType, Lower>(m3).solve(spB);
+    VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: SimplicialCholesky solve, multiple sparse rhs");
 //     
-//     spX = SimplicialCholesky<SparseMatrixType, Upper>(m3).solve(spB);
-//     VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
+    spX = SimplicialCholesky<SparseMatrixType, Upper>(m3).solve(spB);
+    VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: SimplicialCholesky solve, multiple sparse rhs");
   }
   
   
diff --git a/unsupported/test/sparse_llt.h b/unsupported/test/sparse_llt.h
deleted file mode 100644
index 0c3cbe1cd..000000000
--- a/unsupported/test/sparse_llt.h
+++ /dev/null
@@ -1,92 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
-//
-// Eigen is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 3 of the License, or (at your option) any later version.
-//
-// Alternatively, you can redistribute it and/or
-// modify it under the terms of the GNU General Public License as
-// published by the Free Software Foundation; either version 2 of
-// the License, or (at your option) any later version.
-//
-// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
-// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License and a copy of the GNU General Public License along with
-// Eigen. If not, see <http://www.gnu.org/licenses/>.
-
-#include "sparse.h"
-#include <Eigen/SparseExtra>
-
-template<typename LLtSolver, typename Rhs, typename DenseMat, typename DenseRhs>
-void check_sllt(LLtSolver& llt, const typename LLtSolver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db)
-{
-  typedef typename LLtSolver::MatrixType Mat;
-  typedef typename Mat::Scalar Scalar;
-
-  DenseRhs refX = dA.ldlt().solve(db);
-  //Scalar refDet = dA.determinant();
-
-  Rhs x(b.rows(), b.cols());
-  Rhs oldb = b;
-
-  llt.compute(A);
-  if (llt.info() != Success)
-  {
-    std::cerr << "sparse LLt: factorization failed\n";
-    return;
-  }
-  x = llt.solve(b);
-  if (llt.info() != Success)
-  {
-    std::cerr << "sparse LLt: solving failed\n";
-    return;
-  }
-  VERIFY(oldb.isApprox(b) && "sparse LLt: the rhs should not be modified!");
-
-  VERIFY(refX.isApprox(x,test_precision<Scalar>()));
-
-  if(A.cols()<30)
-  {
-    //std::cout << refDet << " == " << lu.determinant() << "\n";
-    //VERIFY_IS_APPROX(refDet,llt.determinant());
-  }
-}
-
-template<typename LLtSolver> void sparse_llt(LLtSolver& llt)
-{
-  typedef typename LLtSolver::MatrixType Mat;
-  typedef typename Mat::Scalar Scalar;
-  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
-  typedef Matrix<Scalar,Dynamic,1> DenseVector;
-
-  std::vector<Vector2i> zeroCoords;
-  std::vector<Vector2i> nonzeroCoords;
-
-  int size = internal::random<int>(1,300);
-  double density = (std::max)(8./(size*size), 0.01);
-  //int rhsSize = internal::random<int>(1,10);
-
-  Mat m2(size, size);
-  DenseMatrix refMat2(size, size);
-
-  DenseVector b = DenseVector::Random(size);
-  DenseVector refX(size), x(size);
-
-  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
-
-  Mat m3 = m2 * m2.adjoint(), m3_half(size,size);
-  DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
-
-  m3_half.template selfadjointView<LLtSolver::UpLo>().rankUpdate(m2,0);
-
-  check_sllt(llt, m3, b, refMat3, b);
-  check_sllt(llt, m3_half, b, refMat3, b);
-}
diff --git a/unsupported/test/sparse_lu.h b/unsupported/test/sparse_lu.h
deleted file mode 100644
index 579d7de66..000000000
--- a/unsupported/test/sparse_lu.h
+++ /dev/null
@@ -1,90 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
-//
-// Eigen is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 3 of the License, or (at your option) any later version.
-//
-// Alternatively, you can redistribute it and/or
-// modify it under the terms of the GNU General Public License as
-// published by the Free Software Foundation; either version 2 of
-// the License, or (at your option) any later version.
-//
-// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
-// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License and a copy of the GNU General Public License along with
-// Eigen. If not, see <http://www.gnu.org/licenses/>.
-
-#include "sparse.h"
-#include <Eigen/SparseExtra>
-
-template<typename LUSolver, typename Rhs, typename DenseMat, typename DenseRhs>
-void check_slu(LUSolver& lu, const typename LUSolver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db)
-{
-  typedef typename LUSolver::MatrixType Mat;
-  typedef typename Mat::Scalar Scalar;
-
-  DenseRhs refX = dA.lu().solve(db);
-  Scalar refDet = dA.determinant();
-
-  Rhs x(b.rows(), b.cols());
-  Rhs oldb = b;
-
-  lu.compute(A);
-  if (lu.info() != Success)
-  {
-    std::cerr << "sparse LU: factorization failed\n";
-    return;
-  }
-  x = lu.solve(b);
-  if (lu.info() != Success)
-  {
-    std::cerr << "sparse LU: solving failed\n";
-    return;
-  }
-  VERIFY(oldb.isApprox(b) && "sparse LU: the rhs should not be modified!");
-
-  VERIFY(refX.isApprox(x,test_precision<Scalar>()));
-
-  if(A.cols()<30)
-  {
-    //std::cout << refDet << " == " << lu.determinant() << "\n";
-    VERIFY_IS_APPROX(refDet,lu.determinant());
-  }
-}
-
-template<typename LUSolver> void sparse_lu(LUSolver& lu)
-{
-  typedef typename LUSolver::MatrixType Mat;
-  typedef typename Mat::Scalar Scalar;
-  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
-  typedef Matrix<Scalar,Dynamic,1> DenseVector;
-
-  std::vector<Vector2i> zeroCoords;
-  std::vector<Vector2i> nonzeroCoords;
-
-  int size = internal::random<int>(1,300);
-  double density = (std::max)(8./(size*size), 0.01);
-  //int rhsSize = internal::random<int>(1,10);
-
-  Mat m2(size, size);
-  DenseMatrix refMat2(size, size);
-
-  DenseVector b = DenseVector::Random(size);
-  DenseVector refX(size), x(size);
-
-  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
-  check_slu(lu, m2, b, refMat2, b);
-
-  refMat2.setZero();
-  m2.setZero();
-  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
-  check_slu(lu, m2, b, refMat2, b);
-}
diff --git a/unsupported/test/sparse_solver.h b/unsupported/test/sparse_solver.h
new file mode 100644
index 000000000..d803ef539
--- /dev/null
+++ b/unsupported/test/sparse_solver.h
@@ -0,0 +1,193 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "sparse.h"
+#include <Eigen/SparseExtra>
+
+template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs>
+void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+
+  DenseRhs refX = dA.lu().solve(db);
+
+  Rhs x(b.rows(), b.cols());
+  Rhs oldb = b;
+
+  solver.compute(A);
+  if (solver.info() != Success)
+  {
+    std::cerr << "sparse SPD: factorization failed (check_sparse_solving)\n";
+    exit(0);
+    return;
+  }
+  x = solver.solve(b);
+  if (solver.info() != Success)
+  {
+    std::cerr << "sparse SPD: solving failed\n";
+    return;
+  }
+  VERIFY(oldb.isApprox(b) && "sparse SPD: the rhs should not be modified!");
+
+  VERIFY(x.isApprox(refX,test_precision<Scalar>()));
+}
+
+template<typename Solver, typename DenseMat>
+void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef typename Mat::RealScalar RealScalar;
+  
+  solver.compute(A);
+  if (solver.info() != Success)
+  {
+    std::cerr << "sparse SPD: factorization failed (check_sparse_determinant)\n";
+    return;
+  }
+
+  Scalar refDet = dA.determinant();
+  VERIFY_IS_APPROX(refDet,solver.determinant());
+}
+
+
+template<typename Solver, typename DenseMat>
+int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+
+  int size = internal::random<int>(1,maxSize);
+  double density = (std::max)(8./(size*size), 0.01);
+
+  Mat M(size, size);
+  DenseMatrix dM(size, size);
+
+  initSparse<Scalar>(density, dM, M, ForceNonZeroDiag);
+
+  A = M * M.adjoint();
+  dA = dM * dM.adjoint();
+  
+  halfA.resize(size,size);
+  halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M);
+  
+  return size;
+}
+
+template<typename Solver> void check_sparse_spd_solving(Solver& solver)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+
+  // generate the problem
+  Mat A, halfA;
+  DenseMatrix dA;
+  int size = generate_sparse_spd_problem(solver, A, halfA, dA);
+
+  // generate the right hand sides
+  int rhsCols = internal::random<int>(1,16);
+  double density = (std::max)(8./(size*rhsCols), 0.1);
+  Mat B(size,rhsCols);
+  DenseVector b = DenseVector::Random(size);
+  DenseMatrix dB(size,rhsCols);
+  initSparse<Scalar>(density, dB, B);
+
+  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);
+}
+
+template<typename Solver> void check_sparse_spd_determinant(Solver& solver)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+
+  // generate the problem
+  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 );
+}
+
+template<typename Solver, typename DenseMat>
+int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+
+  int size = internal::random<int>(1,maxSize);
+  double density = (std::max)(8./(size*size), 0.01);
+  
+  A.resize(size,size);
+  dA.resize(size,size);
+
+  initSparse<Scalar>(density, dA, A, ForceNonZeroDiag);
+  
+  return size;
+}
+
+template<typename Solver> void check_sparse_square_solving(Solver& solver)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+
+  int rhsCols = internal::random<int>(1,16);
+
+  Mat A;
+  DenseMatrix dA;
+  int size = generate_sparse_square_problem(solver, A, dA);
+
+  DenseVector b = DenseVector::Random(size);
+  DenseMatrix dB = DenseMatrix::Random(size,rhsCols);
+
+  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)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+
+  // generate the problem
+  Mat A;
+  DenseMatrix dA;
+  generate_sparse_square_problem(solver, A, dA, 30);
+
+  check_sparse_determinant(solver, A, dA);
+}
diff --git a/unsupported/test/superlu_support.cpp b/unsupported/test/superlu_support.cpp
index 31d623614..e54d72c1c 100644
--- a/unsupported/test/superlu_support.cpp
+++ b/unsupported/test/superlu_support.cpp
@@ -22,7 +22,7 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#include "sparse_lu.h"
+#include "sparse_solver.h"
 
 #ifdef EIGEN_SUPERLU_SUPPORT
 #include <Eigen/SuperLUSupport>
@@ -33,7 +33,9 @@ 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(sparse_lu(superlu_double_colmajor));
-    CALL_SUBTEST_1(sparse_lu(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/unsupported/test/umfpack_support.cpp b/unsupported/test/umfpack_support.cpp
index a2dde2196..16f688302 100644
--- a/unsupported/test/umfpack_support.cpp
+++ b/unsupported/test/umfpack_support.cpp
@@ -22,7 +22,7 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#include "sparse_lu.h"
+#include "sparse_solver.h"
 
 #ifdef EIGEN_UMFPACK_SUPPORT
 #include <Eigen/UmfPackSupport>
@@ -33,7 +33,9 @@ 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(sparse_lu(umfpack_double_colmajor));
-    CALL_SUBTEST_1(sparse_lu(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));
   }
 }