Sparse module:

 * enable complex support for the CHOLMOD LLT backend
   using CHOLMOD's triangular solver
 * quick fix for complex support in SparseLLT::solve

1 files changed, 47 insertions, 22 deletions

diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp
index da9a245a2..8ce4e2264 100644
--- a/test/sparse_solvers.cpp
+++ b/test/sparse_solvers.cpp
@@ -24,6 +24,27 @@
 
 #include "sparse.h"
 
+template<typename Scalar> void
+initSPD(double density,
+        Matrix<Scalar,Dynamic,Dynamic>& refMat,
+        SparseMatrix<Scalar>& sparseMat)
+{
+  Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols());
+  initSparse(density,refMat,sparseMat);
+  refMat = refMat * refMat.adjoint();
+  for (int k=0; k<2; ++k)
+  {
+    initSparse(density,aux,sparseMat,ForceNonZeroDiag);
+    refMat += aux * aux.adjoint();
+  }
+  sparseMat.startFill();
+  for (int j=0 ; j<sparseMat.cols(); ++j)
+    for (int i=j ; i<sparseMat.rows(); ++i)
+      if (refMat(i,j)!=Scalar(0))
+        sparseMat.fill(i,j) = refMat(i,j);
+  sparseMat.endFill();
+}
+
 template<typename Scalar> void sparse_solvers(int rows, int cols)
 {
   double density = std::max(8./(rows*cols), 0.01);
@@ -56,7 +77,6 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
   }
 
   // test LLT
-  if (!NumTraits<Scalar>::IsComplex)
   {
     // TODO fix the issue with complex (see SparseLLT::solveInPlace)
     SparseMatrix<Scalar> m2(rows, cols);
@@ -65,49 +85,54 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
     DenseVector b = DenseVector::Random(cols);
     DenseVector refX(cols), x(cols);
 
-    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
-    refMat2 += refMat2.adjoint();
-    refMat2.diagonal() *= 0.5;
+    initSPD(density, refMat2, m2);
 
     refMat2.llt().solve(b, &refX);
     typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
-    //VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
+    if (!NumTraits<Scalar>::IsComplex)
+    {
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
+    }
     #ifdef EIGEN_CHOLMOD_SUPPORT
     x = b;
     SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x);
     VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
     #endif
-    #ifdef EIGEN_TAUCS_SUPPORT
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
-    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
-    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
-    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
-    #endif
+    if (!NumTraits<Scalar>::IsComplex)
+    {
+      #ifdef EIGEN_TAUCS_SUPPORT
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
+      #endif
+    }
   }
 
   // test LDLT
   if (!NumTraits<Scalar>::IsComplex)
   {
-    // TODO fix the issue with complex (see SparseLDT::solveInPlace)
+    // TODO fix the issue with complex (see SparseLDLT::solveInPlace)
     SparseMatrix<Scalar> m2(rows, cols);
     DenseMatrix refMat2(rows, cols);
 
     DenseVector b = DenseVector::Random(cols);
     DenseVector refX(cols), x(cols);
 
-    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
+    //initSPD(density, refMat2, m2);
+    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
     refMat2 += refMat2.adjoint();
     refMat2.diagonal() *= 0.5;
 
     refMat2.ldlt().solve(b, &refX);
-    typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
+    typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
     x = b;
     SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
     if (ldlt.succeeded())
@@ -177,6 +202,6 @@ void test_sparse_solvers()
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST( sparse_solvers<double>(8, 8) );
     CALL_SUBTEST( sparse_solvers<std::complex<double> >(16, 16) );
-    CALL_SUBTEST( sparse_solvers<double>(33, 33) );
+    CALL_SUBTEST( sparse_solvers<double>(101, 101) );
   }
 }