clean a bit AMD and SimplicialCholesky and add support for partly stored selfadjoint matrices

5 files changed, 230 insertions, 236 deletions

diff --git a/unsupported/Eigen/SparseExtra b/unsupported/Eigen/SparseExtra
index 45b154339..3611d6f32 100644
--- a/unsupported/Eigen/SparseExtra
+++ b/unsupported/Eigen/SparseExtra
@@ -51,7 +51,7 @@ enum {
   OrderingMask                = 0x0f00
 };
 
-#include "src/misc/Solve.h"
+#include "../../Eigen/src/misc/Solve.h"
 
 #include "src/SparseExtra/RandomSetter.h"
 #include "src/SparseExtra/Solve.h"
diff --git a/unsupported/Eigen/src/SparseExtra/Amd.h b/unsupported/Eigen/src/SparseExtra/Amd.h
index a716305fd..f52183232 100644
--- a/unsupported/Eigen/src/SparseExtra/Amd.h
+++ b/unsupported/Eigen/src/SparseExtra/Amd.h
@@ -96,18 +96,14 @@ Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Ind
   return k;
 }
 
-/** keeps off-diagonal entries; drops diagonal entries */
-template<typename Index, typename Scalar>
-struct keep_diag {
-  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
-  {
-    return row!=col;
-  }
-};
 
-/** p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
-template<typename Scalar, int Options, typename Index>
-int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index>& A)  /* order 0:natural, 1:Chol, 2:LU, 3:QR */
+/** \internal
+  * Approximate minimum degree ordering algorithm.
+  * \returns the permutation P reducing the fill-in of the input matrix \a C
+  * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional.
+  * On exit the values of C are destroyed */
+template<typename Scalar, typename Index>
+void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic>& perm)
 {
   typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
   
@@ -115,49 +111,13 @@ int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index>
       k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
       ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
   unsigned int h;
-  /* --- Construct matrix C ----------------------------------------------- */
-  if(order <= 0 || order > 3) return (NULL); /* check */
   
-  Index m = A.rows();
-  Index n = A.cols();
+  Index n = C.cols();
   dense = std::max<Index> (16, 10 * sqrt ((double) n));   /* find dense threshold */
   dense = std::min<Index> (n-2, dense);
-  CCS C;
-  if(order == 1 && n == m) // Cholesky
-  {
-    C = A + SparseMatrix<Scalar,Options,Index>(A.adjoint());
-  }
-  else if(order == 2)  // LU
-  {
-    CCS AT = A.adjoint();
-    // drop dense columns from AT
-    Index* ATp = AT._outerIndexPtr();
-    Index* ATi = AT._innerIndexPtr();
-    Index p2;
-    Index j;
-    for(p2 = 0, j = 0; j < m; j++)
-    {
-      Index p = ATp[j];                         // column j of AT starts here
-      ATp[j] = p2;                              // new column j starts here
-      if(ATp[j+1] - p > dense) continue;        // skip dense col j
-      for(; p < ATp[j+1]; p++)
-        ATi[p2++] = ATi[p];
-    }
-    ATp[m] = p2;                                // finalize AT
-    // TODO this could be implemented using a sparse filter expression
-    // TODO do a cheap selfadjoint rank update
-    C = AT * AT.adjoint();                    // C=A'*A with no dense rows
-  }
-  else  // QR
-  {
-    C = A.adjoint() * A;
-  }
-    
-  C.prune(keep_diag<Index,Scalar>());
-  
   
-  Index cnz = A.nonZeros();
-  Index* P = new Index[n+1];     /* allocate result */
+  Index cnz = C.nonZeros();
+  perm.resize(n+1);
   t = cnz + cnz/5 + 2*n;                 /* add elbow room to C */
   C.resizeNonZeros(t);
   
@@ -170,7 +130,7 @@ int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index>
   Index* degree  = W + 5*(n+1);
   Index* w       = W + 6*(n+1);
   Index* hhead   = W + 7*(n+1);
-  Index* last    = P;                              /* use P as workspace for last */
+  Index* last    = perm.indices().data();                              /* use P as workspace for last */
   
   /* --- Initialize quotient graph ---------------------------------------- */
   Index* Cp = C._outerIndexPtr();
@@ -475,11 +435,12 @@ int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index>
   }
   for(k = 0, i = 0; i <= n; i++)       /* postorder the assembly tree */
   {
-    if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, P, w);
+    if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, perm.indices().data(), w);
   }
+  
+  perm.indices().conservativeResize(n);
 
   delete[] W;
-  return P;
 }
 
 } // namespace internal
diff --git a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
index 14f18d30f..3e502c0aa 100644
--- a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
+++ b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
@@ -177,6 +177,7 @@ class CholmodDecomposition
       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
     {
       cholmod_start(&m_cholmod);
+      setMode(CholmodLDLt);
     }
 
     CholmodDecomposition(const MatrixType& matrix)
@@ -351,6 +352,10 @@ class CholmodDecomposition
       cholmod_free_sparse(&x_cs, &m_cholmod);
     }
     #endif // EIGEN_PARSED_BY_DOXYGEN
+    
+    template<typename Stream>
+    void dumpMemory(Stream& s)
+    {}
 
   protected:
     mutable cholmod_common m_cholmod;
diff --git a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
index 6b22d1502..302be1cce 100644
--- a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
+++ b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
@@ -71,7 +71,7 @@ enum SimplicialCholeskyMode {
 /** \brief A direct sparse Cholesky factorization
   *
   * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization.
-  * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices
+  * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
   * X and B can be either dense or sparse.
   *
   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
@@ -180,73 +180,8 @@ class SimplicialCholesky
       * 
       * \sa factorize()
       */
-    void analyzePattern(const MatrixType& a)
-    {
-      eigen_assert(a.rows()==a.cols());
-      const Index size = a.rows();
-      m_matrix.resize(size, size);
-      m_parent.resize(size);
-      m_nonZerosPerCol.resize(size);
-      
-      Index* tags = ei_aligned_stack_new(Index, size);
-      
-      // TODO allows to configure the permutation
-      const Index* P = internal::minimum_degree_ordering(1, a);
-      const Index* Pinv = 0;
-      if(P)
-      {
-        m_P.indices() = VectorXi::Map(P,size);
-        m_Pinv        = m_P.inverse();
-        Pinv          = m_Pinv.indices().data();
-      }
-      else
-      {
-        m_P.resize(0);
-        m_Pinv.resize(0);
-      }
-
-      for(Index k = 0; k < size; ++k)
-      {
-        /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
-        m_parent[k] = -1;             /* parent of k is not yet known */
-        tags[k] = k;                  /* mark node k as visited */
-        m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
-        Index kk = P ? P[k] : k;      /* kth original, or permuted, column */
-        for(typename MatrixType::InnerIterator it(a,kk); it; ++it)
-        {
-          /* A (i,k) is nonzero (original or permuted A) */
-          Index i = Pinv ? Pinv[it.index()] : it.index();
-          if(i < k)
-          {
-            /* follow path from i to root of etree, stop at flagged node */
-            for(; tags[i] != k; i = m_parent[i])
-            {
-              /* find parent of i if not yet determined */
-              if (m_parent[i] == -1)
-                m_parent[i] = k;
-              ++m_nonZerosPerCol[i];        /* L (k,i) is nonzero */
-              tags[i] = k;                  /* mark i as visited */
-            }
-          }
-        }
-      }
-      
-      // release worspace
-      ei_aligned_stack_delete(Index, tags, size);
-      
-      /* construct Lp index array from m_nonZerosPerCol column counts */
-      Index* Lp = m_matrix._outerIndexPtr();
-      Lp[0] = 0;
-      for(Index k = 0; k < size; ++k)
-        Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (m_LDLt ? 0 : 1);
-
-      m_matrix.resizeNonZeros(Lp[size]);
-      
-      m_isInitialized     = true;
-      m_info              = Success;
-      m_analysisIsOk      = true;
-      m_factorizationIsOk = false;
-    }
+    void analyzePattern(const MatrixType& a);
+    
     
     /** Performs a numeric decomposition of \a matrix
       *
@@ -254,111 +189,17 @@ class SimplicialCholesky
       *
       * \sa analyzePattern()
       */
-    void factorize(const MatrixType& a)
-    {
-      eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
-      eigen_assert(a.rows()==a.cols());
-      const Index size = a.rows();
-      eigen_assert(m_parent.size()==size);
-      eigen_assert(m_nonZerosPerCol.size()==size);
-
-      const Index* Lp = m_matrix._outerIndexPtr();
-      Index* Li = m_matrix._innerIndexPtr();
-      Scalar* Lx = m_matrix._valuePtr();
-
-      Scalar* y       = ei_aligned_stack_new(Scalar, size);
-      Index* pattern  = ei_aligned_stack_new(Index, size);
-      Index* tags     = ei_aligned_stack_new(Index, size);
-      
-      Index* P    = 0;
-      Index* Pinv = 0;
-      
-      if(m_P.size()==size)
-      {
-        P     = m_P.indices().data();
-        Pinv  = m_Pinv.indices().data();
-      }
-      
-      bool ok = true;
-      
-      m_diag.resize(m_LDLt ? size : 0);
-      
-      for(Index k = 0; k < size; ++k)
-      {
-        /* compute nonzero pattern of kth row of L, in topological order */
-        y[k] = 0.0;                     /* Y(0:k) is now all zero */
-        Index top = size;               /* stack for pattern is empty */
-        tags[k] = k;                    /* mark node k as visited */
-        m_nonZerosPerCol[k] = 0;        /* count of nonzeros in column k of L */
-        Index kk = (P) ? (P[k]) : (k);  /* kth original, or permuted, column */
-        for(typename MatrixType::InnerIterator it(a,kk); it; ++it)
-        {
-          Index i = Pinv ? Pinv[it.index()] : it.index();
-          if(i <= k)
-          {
-            y[i] += internal::conj(it.value());            /* scatter A(i,k) into Y (sum duplicates) */
-            Index len;
-            for(len = 0; tags[i] != k; i = m_parent[i])
-            {
-              pattern[len++] = i;     /* L(k,i) is nonzero */
-              tags[i] = k;            /* mark i as visited */
-            }
-            while(len > 0)
-              pattern[--top] = pattern[--len];
-          }
-        }
-
-        /* compute numerical values kth row of L (a sparse triangular solve) */
-        Scalar d = y[k];                  // get D(k,k) and clear Y(k)
-        y[k] = 0.0;
-        for(; top < size; ++top)
-        {
-          if(1)
-          {
-            Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
-            Scalar yi = y[i];             /* get and clear Y(i) */
-            y[i] = 0.0;
-            
-            /* the nonzero entry L(k,i) */
-            Scalar l_ki;
-            if(m_LDLt)
-              l_ki = yi / m_diag[i];       
-            else
-              yi = l_ki = yi / Lx[Lp[i]];
-            
-            Index p2 = Lp[i] + m_nonZerosPerCol[i];
-            Index p;
-            for(p = Lp[i] + (m_LDLt ? 0 : 1); p < p2; ++p)
-              y[Li[p]] -= internal::conj(Lx[p]) * yi;
-            d -= l_ki * internal::conj(yi);
-            Li[p] = k;                          /* store L(k,i) in column form of L */
-            Lx[p] = l_ki;
-            ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
-          }
-        }
-        if(m_LDLt)
-          m_diag[k] = d;
-        else
-        {
-          Index p = Lp[k]+m_nonZerosPerCol[k]++;
-          Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
-          Lx[p] = internal::sqrt(d) ;
-        }
-        if(d == Scalar(0))
-        {
-          ok = false;                         /* failure, D(k,k) is zero */
-          break;
-        }
-      }
-
-      // release workspace
-      ei_aligned_stack_delete(Scalar, y, size);
-      ei_aligned_stack_delete(Index, pattern, size);
-      ei_aligned_stack_delete(Index, tags, size);
-      
-      m_info = ok ? Success : NumericalIssue;
-      m_factorizationIsOk = true;
-    }
+    void factorize(const MatrixType& a);
+    
+    /** \returns the permutation P
+      * \sa permutationPinv() */
+    const PermutationMatrix<Dynamic>& permutationP() const
+    { return m_P; }
+    
+    /** \returns the inverse P^-1 of the permutation P
+      * \sa permutationP() */
+    const PermutationMatrix<Dynamic>& permutationPinv() const
+    { return m_Pinv; }
     
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** \internal */
@@ -408,8 +249,29 @@ class SimplicialCholesky
     }
     */
     #endif // EIGEN_PARSED_BY_DOXYGEN
+    
+    template<typename Stream>
+    void dumpMemory(Stream& s)
+    {
+      int total = 0;
+      s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
+      s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
+      s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
+      s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
+      s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
+      s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
+      s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
+    }
 
   protected:
+    /** keeps off-diagonal entries; drops diagonal entries */
+    struct keep_diag {
+      inline bool operator() (const Index& row, const Index& col, const Scalar&) const
+      {
+        return row!=col;
+      }
+    };
+
     mutable ComputationInfo m_info;
     bool m_isInitialized;
     bool m_factorizationIsOk;
@@ -424,6 +286,172 @@ class SimplicialCholesky
     PermutationMatrix<Dynamic> m_Pinv;  // the inverse permutation
 };
 
+template<typename _MatrixType, int _UpLo>
+void SimplicialCholesky<_MatrixType,_UpLo>::analyzePattern(const MatrixType& a)
+{
+  eigen_assert(a.rows()==a.cols());
+  const Index size = a.rows();
+  m_matrix.resize(size, size);
+  m_parent.resize(size);
+  m_nonZerosPerCol.resize(size);
+  
+  Index* tags = ei_aligned_stack_new(Index, size);
+  
+  // TODO allows to configure the permutation
+  {
+    CholMatrixType C;
+    C = a.template selfadjointView<UpLo>();
+    // remove diagonal entries:
+    C.prune(keep_diag());
+    internal::minimum_degree_ordering(C, m_P);
+  }
+  
+  if(m_P.size()>0)
+    m_Pinv  = m_P.inverse();
+  else
+    m_Pinv.resize(0);
+  
+  SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
+  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
+  
+  for(Index k = 0; k < size; ++k)
+  {
+    /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
+    m_parent[k] = -1;             /* parent of k is not yet known */
+    tags[k] = k;                  /* mark node k as visited */
+    m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
+    for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
+    {
+      Index i = it.index();
+      if(i < k)
+      {
+        /* follow path from i to root of etree, stop at flagged node */
+        for(; tags[i] != k; i = m_parent[i])
+        {
+          /* find parent of i if not yet determined */
+          if (m_parent[i] == -1)
+            m_parent[i] = k;
+          m_nonZerosPerCol[i]++;        /* L (k,i) is nonzero */
+          tags[i] = k;                  /* mark i as visited */
+        }
+      }
+    }
+  }
+  
+  // release workspace
+  ei_aligned_stack_delete(Index, tags, size);
+  
+  /* construct Lp index array from m_nonZerosPerCol column counts */
+  Index* Lp = m_matrix._outerIndexPtr();
+  Lp[0] = 0;
+  for(Index k = 0; k < size; ++k)
+    Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (m_LDLt ? 0 : 1);
+
+  m_matrix.resizeNonZeros(Lp[size]);
+  
+  m_isInitialized     = true;
+  m_info              = Success;
+  m_analysisIsOk      = true;
+  m_factorizationIsOk = false;
+}
+
+
+template<typename _MatrixType, int _UpLo>
+void SimplicialCholesky<_MatrixType,_UpLo>::factorize(const MatrixType& a)
+{
+  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+  eigen_assert(a.rows()==a.cols());
+  const Index size = a.rows();
+  eigen_assert(m_parent.size()==size);
+  eigen_assert(m_nonZerosPerCol.size()==size);
+
+  const Index* Lp = m_matrix._outerIndexPtr();
+  Index* Li = m_matrix._innerIndexPtr();
+  Scalar* Lx = m_matrix._valuePtr();
+
+  Scalar* y       = ei_aligned_stack_new(Scalar, size);
+  Index* pattern  = ei_aligned_stack_new(Index, size);
+  Index* tags     = ei_aligned_stack_new(Index, size);
+
+  SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
+  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
+  
+  bool ok = true;
+  m_diag.resize(m_LDLt ? size : 0);
+  
+  for(Index k = 0; k < size; ++k)
+  {
+    // compute nonzero pattern of kth row of L, in topological order
+    y[k] = 0.0;                     // Y(0:k) is now all zero
+    Index top = size;               // stack for pattern is empty
+    tags[k] = k;                    // mark node k as visited
+    m_nonZerosPerCol[k] = 0;        // count of nonzeros in column k of L
+    for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
+    {
+      Index i = it.index();
+      if(i <= k)
+      {
+        y[i] += internal::conj(it.value());            /* scatter A(i,k) into Y (sum duplicates) */
+        Index len;
+        for(len = 0; tags[i] != k; i = m_parent[i])
+        {
+          pattern[len++] = i;     /* L(k,i) is nonzero */
+          tags[i] = k;            /* mark i as visited */
+        }
+        while(len > 0)
+          pattern[--top] = pattern[--len];
+      }
+    }
+
+    /* compute numerical values kth row of L (a sparse triangular solve) */
+    Scalar d = y[k];                  // get D(k,k) and clear Y(k)
+    y[k] = 0.0;
+    for(; top < size; ++top)
+    {
+      Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
+      Scalar yi = y[i];             /* get and clear Y(i) */
+      y[i] = 0.0;
+      
+      /* the nonzero entry L(k,i) */
+      Scalar l_ki;
+      if(m_LDLt)
+        l_ki = yi / m_diag[i];       
+      else
+        yi = l_ki = yi / Lx[Lp[i]];
+      
+      Index p2 = Lp[i] + m_nonZerosPerCol[i];
+      Index p;
+      for(p = Lp[i] + (m_LDLt ? 0 : 1); p < p2; ++p)
+        y[Li[p]] -= internal::conj(Lx[p]) * yi;
+      d -= l_ki * internal::conj(yi);
+      Li[p] = k;                          /* store L(k,i) in column form of L */
+      Lx[p] = l_ki;
+      ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
+    }
+    if(m_LDLt)
+      m_diag[k] = d;
+    else
+    {
+      Index p = Lp[k]+m_nonZerosPerCol[k]++;
+      Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
+      Lx[p] = internal::sqrt(d) ;
+    }
+    if(d == Scalar(0))
+    {
+      ok = false;                         /* failure, D(k,k) is zero */
+      break;
+    }
+  }
+
+  // release workspace
+  ei_aligned_stack_delete(Scalar, y, size);
+  ei_aligned_stack_delete(Index, pattern, size);
+  ei_aligned_stack_delete(Index, tags, size);
+  
+  m_info = ok ? Success : NumericalIssue;
+  m_factorizationIsOk = true;
+}
+
 namespace internal {
   
 template<typename _MatrixType, int _UpLo, typename Rhs>
diff --git a/unsupported/test/sparse_ldlt.cpp b/unsupported/test/sparse_ldlt.cpp
index 035e21123..07666273c 100644
--- a/unsupported/test/sparse_ldlt.cpp
+++ b/unsupported/test/sparse_ldlt.cpp
@@ -96,30 +96,30 @@ template<typename Scalar> void sparse_ldlt(int rows, int cols)
     x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(b);
     VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, single dense rhs");
     
-//     x = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3_lo).solve(b);
-//     VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, lower only, single dense rhs");
+    x = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3_lo).solve(b);
+    VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, lower only, single dense rhs");
     
-//     x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>(m3_up).solve(b);
-//     VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, upper only, single dense rhs");
+    x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>(m3_up).solve(b);
+    VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, upper only, single dense rhs");
     
     
     // with multiple rhs
     ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B);
 
-    X = SimplicialCholesky<SparseMatrix<Scalar>, Lower>()/*.setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt)*/.compute(m3).solve(B);
+    X = SimplicialCholesky<SparseMatrix<Scalar>, Lower>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
     VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, lower, multiple dense rhs");
     
-//     X = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
-//     VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, multiple dense rhs");
+    X = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
+    VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, multiple dense rhs");
     
     
-//     // with a sparse rhs
+    // with a sparse rhs
 //     SparseMatrix<Scalar> 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<SparseMatrix<Scalar>, Lower>(m3).solve(spB);
 //     VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
 //