devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

add a minimum degree ordering routine based on CSparse (LGPL) and a new built-in sparse cholesky decomposition

Browse Source
This commit is contained in:
Gael Guennebaud
2010-11-04 09:58:22 +01:00
parent 5a4f77716d
commit 15e8ad686c
7 changed files with 1177 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files

137
unsupported/test/sparse_ldlt.cpp
View File

@@ -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);
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);
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);
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));
// 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");
// 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
}