diff --git a/unsupported/test/conjugate_gradient.cpp b/unsupported/test/conjugate_gradient.cpp index 47cafaac3..91a205a77 100644 --- a/unsupported/test/conjugate_gradient.cpp +++ b/unsupported/test/conjugate_gradient.cpp @@ -22,7 +22,7 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#include "sparse_llt.h" +#include "sparse_solver.h" #include template void test_conjugate_gradient_T() @@ -32,10 +32,10 @@ template void test_conjugate_gradient_T() 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); + CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_diag) ); + CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_diag) ); + CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_I) ); + CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_I) ); } void test_conjugate_gradient() diff --git a/unsupported/test/simplicial_cholesky.cpp b/unsupported/test/simplicial_cholesky.cpp index 00a5e1096..b9a55ef36 100644 --- a/unsupported/test/simplicial_cholesky.cpp +++ b/unsupported/test/simplicial_cholesky.cpp @@ -22,7 +22,7 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#include "sparse_llt.h" +#include "sparse_solver.h" template void test_simplicial_cholesky_T() { @@ -33,12 +33,19 @@ template void test_simplicial_cholesky_T() 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); + check_sparse_spd_solving(chol_colmajor_lower); + check_sparse_spd_solving(chol_colmajor_upper); + check_sparse_spd_solving(llt_colmajor_lower); + check_sparse_spd_solving(llt_colmajor_upper); + check_sparse_spd_solving(ldlt_colmajor_lower); + check_sparse_spd_solving(ldlt_colmajor_upper); + + check_sparse_spd_determinant(chol_colmajor_lower); + check_sparse_spd_determinant(chol_colmajor_upper); + check_sparse_spd_determinant(llt_colmajor_lower); + check_sparse_spd_determinant(llt_colmajor_upper); + check_sparse_spd_determinant(ldlt_colmajor_lower); + check_sparse_spd_determinant(ldlt_colmajor_upper); } void test_simplicial_cholesky() diff --git a/unsupported/test/sparse_ldlt.cpp b/unsupported/test/sparse_ldlt.cpp index c387f78ac..4cb61bf6c 100644 --- a/unsupported/test/sparse_ldlt.cpp +++ b/unsupported/test/sparse_ldlt.cpp @@ -132,17 +132,17 @@ template void sparse_ldlt(int rows, int cols) // with a sparse rhs -// SparseMatrixType spB(rows,cols), spX(rows,cols); -// B.diagonal().array() += 1; -// spB = B.sparseView(0.5,1); + SparseMatrixType spB(rows,cols), spX(rows,cols); + B.diagonal().array() += 1; + spB = B.sparseView(0.5,1); + + ref_X = refMat3.template selfadjointView().llt().solve(DenseMatrix(spB)); + + spX = SimplicialCholesky(m3).solve(spB); + VERIFY(ref_X.isApprox(spX.toDense(),test_precision()) && "LLT: SimplicialCholesky solve, multiple sparse rhs"); // -// ref_X = refMat3.template selfadjointView().llt().solve(DenseMatrix(spB)); -// -// spX = SimplicialCholesky(m3).solve(spB); -// VERIFY(ref_X.isApprox(spX.toDense(),test_precision()) && "LLT: cholmod solve, multiple sparse rhs"); -// -// spX = SimplicialCholesky(m3).solve(spB); -// VERIFY(ref_X.isApprox(spX.toDense(),test_precision()) && "LLT: cholmod solve, multiple sparse rhs"); + spX = SimplicialCholesky(m3).solve(spB); + VERIFY(ref_X.isApprox(spX.toDense(),test_precision()) && "LLT: SimplicialCholesky solve, multiple sparse rhs"); } diff --git a/unsupported/test/sparse_llt.h b/unsupported/test/sparse_llt.h deleted file mode 100644 index 0c3cbe1cd..000000000 --- a/unsupported/test/sparse_llt.h +++ /dev/null @@ -1,92 +0,0 @@ -// 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); -} diff --git a/unsupported/test/sparse_lu.h b/unsupported/test/sparse_lu.h deleted file mode 100644 index 579d7de66..000000000 --- a/unsupported/test/sparse_lu.h +++ /dev/null @@ -1,90 +0,0 @@ -// 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_slu(LUSolver& lu, const typename LUSolver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db) -{ - typedef typename LUSolver::MatrixType Mat; - typedef typename Mat::Scalar Scalar; - - DenseRhs refX = dA.lu().solve(db); - Scalar refDet = dA.determinant(); - - Rhs x(b.rows(), b.cols()); - Rhs oldb = b; - - lu.compute(A); - if (lu.info() != Success) - { - std::cerr << "sparse LU: factorization failed\n"; - return; - } - x = lu.solve(b); - if (lu.info() != Success) - { - std::cerr << "sparse LU: solving failed\n"; - return; - } - VERIFY(oldb.isApprox(b) && "sparse LU: 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,lu.determinant()); - } -} - -template void sparse_lu(LUSolver& lu) -{ - typedef typename LUSolver::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); - check_slu(lu, m2, b, refMat2, b); - - refMat2.setZero(); - m2.setZero(); - initSparse(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords); - check_slu(lu, m2, b, refMat2, b); -} diff --git a/unsupported/test/sparse_solver.h b/unsupported/test/sparse_solver.h new file mode 100644 index 000000000..d803ef539 --- /dev/null +++ b/unsupported/test/sparse_solver.h @@ -0,0 +1,193 @@ +// 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_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + + DenseRhs refX = dA.lu().solve(db); + + Rhs x(b.rows(), b.cols()); + Rhs oldb = b; + + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse SPD: factorization failed (check_sparse_solving)\n"; + exit(0); + return; + } + x = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "sparse SPD: solving failed\n"; + return; + } + VERIFY(oldb.isApprox(b) && "sparse SPD: the rhs should not be modified!"); + + VERIFY(x.isApprox(refX,test_precision())); +} + +template +void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef typename Mat::RealScalar RealScalar; + + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse SPD: factorization failed (check_sparse_determinant)\n"; + return; + } + + Scalar refDet = dA.determinant(); + VERIFY_IS_APPROX(refDet,solver.determinant()); +} + + +template +int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix DenseMatrix; + + int size = internal::random(1,maxSize); + double density = (std::max)(8./(size*size), 0.01); + + Mat M(size, size); + DenseMatrix dM(size, size); + + initSparse(density, dM, M, ForceNonZeroDiag); + + A = M * M.adjoint(); + dA = dM * dM.adjoint(); + + halfA.resize(size,size); + halfA.template selfadjointView().rankUpdate(M); + + return size; +} + +template void check_sparse_spd_solving(Solver& solver) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + // generate the problem + Mat A, halfA; + DenseMatrix dA; + int size = generate_sparse_spd_problem(solver, A, halfA, dA); + + // generate the right hand sides + int rhsCols = internal::random(1,16); + double density = (std::max)(8./(size*rhsCols), 0.1); + Mat B(size,rhsCols); + DenseVector b = DenseVector::Random(size); + DenseMatrix dB(size,rhsCols); + initSparse(density, dB, B); + + check_sparse_solving(solver, A, b, dA, b); + check_sparse_solving(solver, halfA, b, dA, b); + check_sparse_solving(solver, A, dB, dA, dB); + check_sparse_solving(solver, halfA, dB, dA, dB); + check_sparse_solving(solver, A, B, dA, dB); + check_sparse_solving(solver, halfA, B, dA, dB); +} + +template void check_sparse_spd_determinant(Solver& solver) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix DenseMatrix; + + // generate the problem + Mat A, halfA; + DenseMatrix dA; + generate_sparse_spd_problem(solver, A, halfA, dA, 30); + + check_sparse_determinant(solver, A, dA); + check_sparse_determinant(solver, halfA, dA ); +} + +template +int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix DenseMatrix; + + int size = internal::random(1,maxSize); + double density = (std::max)(8./(size*size), 0.01); + + A.resize(size,size); + dA.resize(size,size); + + initSparse(density, dA, A, ForceNonZeroDiag); + + return size; +} + +template void check_sparse_square_solving(Solver& solver) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + + int rhsCols = internal::random(1,16); + + Mat A; + DenseMatrix dA; + int size = generate_sparse_square_problem(solver, A, dA); + + DenseVector b = DenseVector::Random(size); + DenseMatrix dB = DenseMatrix::Random(size,rhsCols); + + check_sparse_solving(solver, A, b, dA, b); + check_sparse_solving(solver, A, dB, dA, dB); +} + +template void check_sparse_square_determinant(Solver& solver) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix DenseMatrix; + + // generate the problem + Mat A; + DenseMatrix dA; + generate_sparse_square_problem(solver, A, dA, 30); + + check_sparse_determinant(solver, A, dA); +} diff --git a/unsupported/test/superlu_support.cpp b/unsupported/test/superlu_support.cpp index 31d623614..e54d72c1c 100644 --- a/unsupported/test/superlu_support.cpp +++ b/unsupported/test/superlu_support.cpp @@ -22,7 +22,7 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#include "sparse_lu.h" +#include "sparse_solver.h" #ifdef EIGEN_SUPERLU_SUPPORT #include @@ -33,7 +33,9 @@ void test_superlu_support() for(int i = 0; i < g_repeat; i++) { SuperLU > superlu_double_colmajor; SuperLU > > superlu_cplxdouble_colmajor; - CALL_SUBTEST_1(sparse_lu(superlu_double_colmajor)); - CALL_SUBTEST_1(sparse_lu(superlu_cplxdouble_colmajor)); + CALL_SUBTEST_1( check_sparse_square_solving(superlu_double_colmajor) ); + CALL_SUBTEST_2( check_sparse_square_solving(superlu_cplxdouble_colmajor) ); + CALL_SUBTEST_1( check_sparse_square_determinant(superlu_double_colmajor) ); + CALL_SUBTEST_2( check_sparse_square_determinant(superlu_cplxdouble_colmajor) ); } } diff --git a/unsupported/test/umfpack_support.cpp b/unsupported/test/umfpack_support.cpp index a2dde2196..16f688302 100644 --- a/unsupported/test/umfpack_support.cpp +++ b/unsupported/test/umfpack_support.cpp @@ -22,7 +22,7 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#include "sparse_lu.h" +#include "sparse_solver.h" #ifdef EIGEN_UMFPACK_SUPPORT #include @@ -33,7 +33,9 @@ void test_umfpack_support() for(int i = 0; i < g_repeat; i++) { UmfPackLU > umfpack_double_colmajor; UmfPackLU > > umfpack_cplxdouble_colmajor; - CALL_SUBTEST_1(sparse_lu(umfpack_double_colmajor)); - CALL_SUBTEST_1(sparse_lu(umfpack_cplxdouble_colmajor)); + CALL_SUBTEST_1(check_sparse_square_solving(umfpack_double_colmajor)); + CALL_SUBTEST_2(check_sparse_square_solving(umfpack_cplxdouble_colmajor)); + CALL_SUBTEST_1(check_sparse_square_determinant(umfpack_double_colmajor)); + CALL_SUBTEST_2(check_sparse_square_determinant(umfpack_cplxdouble_colmajor)); } }