aboutsummaryrefslogtreecommitdiffhomepage
path: root/unsupported/test/sparse_ldlt.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'unsupported/test/sparse_ldlt.cpp')
-rw-r--r--unsupported/test/sparse_ldlt.cpp131
1 files changed, 109 insertions, 22 deletions
diff --git a/unsupported/test/sparse_ldlt.cpp b/unsupported/test/sparse_ldlt.cpp
index 275839670..035e21123 100644
--- a/unsupported/test/sparse_ldlt.cpp
+++ b/unsupported/test/sparse_ldlt.cpp
@@ -31,6 +31,8 @@
template<typename Scalar> void sparse_ldlt(int rows, int cols)
{
+ static bool odd = true;
+ odd = !odd;
double density = std::max(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
@@ -42,41 +44,126 @@ template<typename Scalar> void sparse_ldlt(int rows, int cols)
DenseVector refX(cols), x(cols);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
- for(int i=0; i<rows; ++i)
- m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
-
- refX = refMat2.template selfadjointView<Upper>().ldlt().solve(b);
+
+ SparseMatrix<Scalar> m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
+ DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
+
+ refX = refMat3.template selfadjointView<Upper>().ldlt().solve(b);
typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix;
x = b;
- SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
+ SparseLDLT<SparseSelfAdjointMatrix> ldlt(m3);
if (ldlt.succeeded())
ldlt.solveInPlace(x);
else
std::cerr << "warning LDLT failed\n";
- VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+// VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
+
+
+
+ // new Simplicial LLT
+
+
+ // new API
+ {
+ SparseMatrix<Scalar> m2(rows, cols);
+ DenseMatrix refMat2(rows, cols);
-#ifdef EIGEN_CHOLMOD_SUPPORT
- x = b;
- SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt2(m2);
- if (ldlt2.succeeded())
- ldlt2.solveInPlace(x);
- else
- std::cerr << "warning LDLT failed\n";
+ DenseVector b = DenseVector::Random(cols);
+ DenseVector ref_x(cols), x(cols);
+ DenseMatrix B = DenseMatrix::Random(rows,cols);
+ DenseMatrix ref_X(rows,cols), X(rows,cols);
- VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
- VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solveInPlace");
+ initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0);
+ for(int i=0; i<rows; ++i)
+ m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
+
+
+ SparseMatrix<Scalar> m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
+ DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
+
+ m3_lo.template selfadjointView<Lower>().rankUpdate(m2,0);
+ m3_up.template selfadjointView<Upper>().rankUpdate(m2,0);
+
+ // with a single vector as the rhs
+ ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
- SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt3(m2);
- if (ldlt3.succeeded())
- x = ldlt3.solve(b);
- else
- std::cerr << "warning LDLT failed\n";
+ 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, single 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, 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");
+
+
+ // 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);
+ 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");
+
+
+// // 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));
- VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
- VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solve");
+// spX = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3).solve(spB);
+// VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
+//
+// spX = SimplicialCholesky<SparseMatrix<Scalar>, Upper>(m3).solve(spB);
+// VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
+ }
+
+
+
+// for(int i=0; i<rows; ++i)
+// m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
+//
+// refX = refMat2.template selfadjointView<Upper>().ldlt().solve(b);
+// typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix;
+// x = b;
+// SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
+// if (ldlt.succeeded())
+// ldlt.solveInPlace(x);
+// else
+// std::cerr << "warning LDLT failed\n";
+//
+// VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+// VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
+
+#ifdef EIGEN_CHOLMOD_SUPPORT
+// x = b;
+// SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt2(m2);
+// if (ldlt2.succeeded())
+// ldlt2.solveInPlace(x);
+// else
+// std::cerr << "warning LDLT failed\n";
+//
+// VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+// VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solveInPlace");
+//
+//
+// SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt3(m2);
+// if (ldlt3.succeeded())
+// x = ldlt3.solve(b);
+// else
+// std::cerr << "warning LDLT failed\n";
+//
+// VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+// VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solve");
#endif
}
generated by cgit v1.2.3 (git 2.39.1) at 2025-01-10 08:56:23 +0000