clean a bit AMD and SimplicialCholesky and add support for partly stored selfadjoint matrices
diff --git a/unsupported/Eigen/SparseExtra b/unsupported/Eigen/SparseExtra
index 45b1543..3611d6f 100644
--- a/unsupported/Eigen/SparseExtra
+++ b/unsupported/Eigen/SparseExtra
@@ -51,7 +51,7 @@
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 a716305..f521832 100644
--- a/unsupported/Eigen/src/SparseExtra/Amd.h
+++ b/unsupported/Eigen/src/SparseExtra/Amd.h
@@ -96,18 +96,14 @@
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 @@
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 @@
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 @@
}
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 14f18d3..3e502c0 100644
--- a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
+++ b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
@@ -177,6 +177,7 @@
: m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
{
cholmod_start(&m_cholmod);
+ setMode(CholmodLDLt);
}
CholmodDecomposition(const MatrixType& matrix)
@@ -351,6 +352,10 @@
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 6b22d15..302be1c 100644
--- a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
+++ b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
@@ -71,7 +71,7 @@
/** \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 @@
*
* \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 @@
*
* \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 @@
}
*/
#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 @@
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 035e211..0766627 100644
--- a/unsupported/test/sparse_ldlt.cpp
+++ b/unsupported/test/sparse_ldlt.cpp
@@ -96,30 +96,30 @@
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");
//
