diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 117725c04..9adfe4758 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -92,6 +92,9 @@ ei_add_test(sparse_llt "" "${SPARSE_LIBS}") ei_add_test(sparse_ldlt "" "${SPARSE_LIBS}") ei_add_test(sparse_lu_legacy "" "${SPARSE_LIBS}") +ei_add_test(simplicial_cholesky) +ei_add_test(conjugate_gradient) + if(UMFPACK_FOUND) ei_add_test(umfpack_support "" "${UMFPACK_ALL_LIBS}") endif() diff --git a/unsupported/test/conjugate_gradient.cpp b/unsupported/test/conjugate_gradient.cpp new file mode 100644 index 000000000..47cafaac3 --- /dev/null +++ b/unsupported/test/conjugate_gradient.cpp @@ -0,0 +1,47 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#include "sparse_llt.h" +#include + +template void test_conjugate_gradient_T() +{ + ConjugateGradient, Lower> cg_colmajor_lower_diag; + ConjugateGradient, Upper> cg_colmajor_upper_diag; + ConjugateGradient, Lower, IdentityPreconditioner> cg_colmajor_lower_I; + ConjugateGradient, Upper, IdentityPreconditioner> cg_colmajor_upper_I; + + sparse_llt(cg_colmajor_lower_diag); + sparse_llt(cg_colmajor_upper_diag); + sparse_llt(cg_colmajor_lower_I); + sparse_llt(cg_colmajor_upper_I); +} + +void test_conjugate_gradient() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1(test_conjugate_gradient_T()); + CALL_SUBTEST_2(test_conjugate_gradient_T >()); + } +} diff --git a/unsupported/test/simplicial_cholesky.cpp b/unsupported/test/simplicial_cholesky.cpp new file mode 100644 index 000000000..00a5e1096 --- /dev/null +++ b/unsupported/test/simplicial_cholesky.cpp @@ -0,0 +1,50 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#include "sparse_llt.h" + +template void test_simplicial_cholesky_T() +{ + SimplicialCholesky, Lower> chol_colmajor_lower; + SimplicialCholesky, Upper> chol_colmajor_upper; + SimplicialLLt, Lower> llt_colmajor_lower; + SimplicialLDLt, Upper> llt_colmajor_upper; + SimplicialLDLt, Lower> ldlt_colmajor_lower; + SimplicialLDLt, Upper> ldlt_colmajor_upper; + + sparse_llt(chol_colmajor_lower); + sparse_llt(chol_colmajor_upper); + sparse_llt(llt_colmajor_lower); + sparse_llt(llt_colmajor_upper); + sparse_llt(ldlt_colmajor_lower); + sparse_llt(ldlt_colmajor_upper); +} + +void test_simplicial_cholesky() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1(test_simplicial_cholesky_T()); + CALL_SUBTEST_2(test_simplicial_cholesky_T >()); + } +} diff --git a/unsupported/test/sparse_llt.h b/unsupported/test/sparse_llt.h new file mode 100644 index 000000000..0c3cbe1cd --- /dev/null +++ b/unsupported/test/sparse_llt.h @@ -0,0 +1,92 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#include "sparse.h" +#include + +template +void check_sllt(LLtSolver& llt, const typename LLtSolver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db) +{ + typedef typename LLtSolver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + + DenseRhs refX = dA.ldlt().solve(db); + //Scalar refDet = dA.determinant(); + + Rhs x(b.rows(), b.cols()); + Rhs oldb = b; + + llt.compute(A); + if (llt.info() != Success) + { + std::cerr << "sparse LLt: factorization failed\n"; + return; + } + x = llt.solve(b); + if (llt.info() != Success) + { + std::cerr << "sparse LLt: solving failed\n"; + return; + } + VERIFY(oldb.isApprox(b) && "sparse LLt: the rhs should not be modified!"); + + VERIFY(refX.isApprox(x,test_precision())); + + if(A.cols()<30) + { + //std::cout << refDet << " == " << lu.determinant() << "\n"; + //VERIFY_IS_APPROX(refDet,llt.determinant()); + } +} + +template void sparse_llt(LLtSolver& llt) +{ + typedef typename LLtSolver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + std::vector zeroCoords; + std::vector nonzeroCoords; + + int size = internal::random(1,300); + double density = (std::max)(8./(size*size), 0.01); + //int rhsSize = internal::random(1,10); + + Mat m2(size, size); + DenseMatrix refMat2(size, size); + + DenseVector b = DenseVector::Random(size); + DenseVector refX(size), x(size); + + initSparse(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords); + + Mat m3 = m2 * m2.adjoint(), m3_half(size,size); + DenseMatrix refMat3 = refMat2 * refMat2.adjoint(); + + m3_half.template selfadjointView().rankUpdate(m2,0); + + check_sllt(llt, m3, b, refMat3, b); + check_sllt(llt, m3_half, b, refMat3, b); +}