GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-13 20:26:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

refactor sparse solving unit tests

Browse Source
This commit is contained in:
Gael Guennebaud 2011-10-11 11:32:26 +02:00
parent 4f237f035c
commit 3172749f32
8 changed files with 232 additions and 210 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
unsupported/test/conjugate_gradient.cpp
View File

@ -22,7 +22,7 @@
// License and a copy of the GNU General Public License along with // License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>. // Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "sparse_llt.h" #include "sparse_solver.h"
#include <Eigen/IterativeSolvers> #include <Eigen/IterativeSolvers>
template<typename T> void test_conjugate_gradient_T() template<typename T> void test_conjugate_gradient_T()
@ -32,10 +32,10 @@ template<typename T> void test_conjugate_gradient_T()
ConjugateGradient<SparseMatrix<T>, Lower, IdentityPreconditioner> cg_colmajor_lower_I; ConjugateGradient<SparseMatrix<T>, Lower, IdentityPreconditioner> cg_colmajor_lower_I;
ConjugateGradient<SparseMatrix<T>, Upper, IdentityPreconditioner> cg_colmajor_upper_I; ConjugateGradient<SparseMatrix<T>, Upper, IdentityPreconditioner> cg_colmajor_upper_I;
sparse_llt(cg_colmajor_lower_diag); CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_diag) );
sparse_llt(cg_colmajor_upper_diag); CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_diag) );
sparse_llt(cg_colmajor_lower_I); CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_lower_I) );
sparse_llt(cg_colmajor_upper_I); CALL_SUBTEST( check_sparse_spd_solving(cg_colmajor_upper_I) );
} }
void test_conjugate_gradient() void test_conjugate_gradient()

21
unsupported/test/simplicial_cholesky.cpp
View File

@ -22,7 +22,7 @@
// License and a copy of the GNU General Public License along with // License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>. // Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "sparse_llt.h" #include "sparse_solver.h"
template<typename T> void test_simplicial_cholesky_T() template<typename T> void test_simplicial_cholesky_T()
{ {
@ -33,12 +33,19 @@ template<typename T> void test_simplicial_cholesky_T()
SimplicialLDLt<SparseMatrix<T>, Lower> ldlt_colmajor_lower; SimplicialLDLt<SparseMatrix<T>, Lower> ldlt_colmajor_lower;
SimplicialLDLt<SparseMatrix<T>, Upper> ldlt_colmajor_upper; SimplicialLDLt<SparseMatrix<T>, Upper> ldlt_colmajor_upper;
sparse_llt(chol_colmajor_lower); check_sparse_spd_solving(chol_colmajor_lower);
sparse_llt(chol_colmajor_upper); check_sparse_spd_solving(chol_colmajor_upper);
sparse_llt(llt_colmajor_lower); check_sparse_spd_solving(llt_colmajor_lower);
sparse_llt(llt_colmajor_upper); check_sparse_spd_solving(llt_colmajor_upper);
sparse_llt(ldlt_colmajor_lower); check_sparse_spd_solving(ldlt_colmajor_lower);
sparse_llt(ldlt_colmajor_upper); 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() void test_simplicial_cholesky()

20
unsupported/test/sparse_ldlt.cpp
View File

@ -132,17 +132,17 @@ template<typename Scalar,typename Index> void sparse_ldlt(int rows, int cols)
// with a sparse rhs // with a sparse rhs
// SparseMatrixType spB(rows,cols), spX(rows,cols); SparseMatrixType spB(rows,cols), spX(rows,cols);
// B.diagonal().array() += 1; B.diagonal().array() += 1;
// spB = B.sparseView(0.5,1); spB = B.sparseView(0.5,1);
ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB));
spX = SimplicialCholesky<SparseMatrixType, Lower>(m3).solve(spB);
VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: SimplicialCholesky solve, multiple sparse rhs");
// //
// ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB)); spX = SimplicialCholesky<SparseMatrixType, Upper>(m3).solve(spB);
// VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: SimplicialCholesky solve, multiple sparse rhs");
// spX = SimplicialCholesky<SparseMatrixType, Lower>(m3).solve(spB);
// VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
//
// spX = SimplicialCholesky<SparseMatrixType, Upper>(m3).solve(spB);
// VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
} }

92
unsupported/test/sparse_llt.h
View File

@ -1,92 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
//
// 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 <http://www.gnu.org/licenses/>.
#include "sparse.h"
#include <Eigen/SparseExtra>
template<typename LLtSolver, typename Rhs, typename DenseMat, typename DenseRhs>
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<Scalar>()));
if(A.cols()<30)
{
//std::cout << refDet << " == " << lu.determinant() << "\n";
//VERIFY_IS_APPROX(refDet,llt.determinant());
}
}
template<typename LLtSolver> void sparse_llt(LLtSolver& llt)
{
typedef typename LLtSolver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
std::vector<Vector2i> zeroCoords;
std::vector<Vector2i> nonzeroCoords;
int size = internal::random<int>(1,300);
double density = (std::max)(8./(size*size), 0.01);
//int rhsSize = internal::random<int>(1,10);
Mat m2(size, size);
DenseMatrix refMat2(size, size);
DenseVector b = DenseVector::Random(size);
DenseVector refX(size), x(size);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
Mat m3 = m2 * m2.adjoint(), m3_half(size,size);
DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
m3_half.template selfadjointView<LLtSolver::UpLo>().rankUpdate(m2,0);
check_sllt(llt, m3, b, refMat3, b);
check_sllt(llt, m3_half, b, refMat3, b);
}

90
unsupported/test/sparse_lu.h
View File

@ -1,90 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
//
// 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 <http://www.gnu.org/licenses/>.
#include "sparse.h"
#include <Eigen/SparseExtra>
template<typename LUSolver, typename Rhs, typename DenseMat, typename DenseRhs>
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<Scalar>()));
if(A.cols()<30)
{
//std::cout << refDet << " == " << lu.determinant() << "\n";
VERIFY_IS_APPROX(refDet,lu.determinant());
}
}
template<typename LUSolver> void sparse_lu(LUSolver& lu)
{
typedef typename LUSolver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
std::vector<Vector2i> zeroCoords;
std::vector<Vector2i> nonzeroCoords;
int size = internal::random<int>(1,300);
double density = (std::max)(8./(size*size), 0.01);
//int rhsSize = internal::random<int>(1,10);
Mat m2(size, size);
DenseMatrix refMat2(size, size);
DenseVector b = DenseVector::Random(size);
DenseVector refX(size), x(size);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
check_slu(lu, m2, b, refMat2, b);
refMat2.setZero();
m2.setZero();
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
check_slu(lu, m2, b, refMat2, b);
}

193
unsupported/test/sparse_solver.h Normal file
View File

@ -0,0 +1,193 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
//
// 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 <http://www.gnu.org/licenses/>.
#include "sparse.h"
#include <Eigen/SparseExtra>
template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs>
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<Scalar>()));
}
template<typename Solver, typename DenseMat>
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<typename Solver, typename DenseMat>
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<Scalar,Dynamic,Dynamic> DenseMatrix;
int size = internal::random<int>(1,maxSize);
double density = (std::max)(8./(size*size), 0.01);
Mat M(size, size);
DenseMatrix dM(size, size);
initSparse<Scalar>(density, dM, M, ForceNonZeroDiag);
A = M * M.adjoint();
dA = dM * dM.adjoint();
halfA.resize(size,size);
halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M);
return size;
}
template<typename Solver> void check_sparse_spd_solving(Solver& solver)
{
typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> 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<int>(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<Scalar>(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<typename Solver> void check_sparse_spd_determinant(Solver& solver)
{
typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> 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<typename Solver, typename DenseMat>
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<Scalar,Dynamic,Dynamic> DenseMatrix;
int size = internal::random<int>(1,maxSize);
double density = (std::max)(8./(size*size), 0.01);
A.resize(size,size);
dA.resize(size,size);
initSparse<Scalar>(density, dA, A, ForceNonZeroDiag);
return size;
}
template<typename Solver> void check_sparse_square_solving(Solver& solver)
{
typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
int rhsCols = internal::random<int>(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<typename Solver> void check_sparse_square_determinant(Solver& solver)
{
typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
// generate the problem
Mat A;
DenseMatrix dA;
generate_sparse_square_problem(solver, A, dA, 30);
check_sparse_determinant(solver, A, dA);
}

8
unsupported/test/superlu_support.cpp
View File

@ -22,7 +22,7 @@
// License and a copy of the GNU General Public License along with // License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>. // Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "sparse_lu.h" #include "sparse_solver.h"
#ifdef EIGEN_SUPERLU_SUPPORT #ifdef EIGEN_SUPERLU_SUPPORT
#include <Eigen/SuperLUSupport> #include <Eigen/SuperLUSupport>
@ -33,7 +33,9 @@ void test_superlu_support()
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
SuperLU<SparseMatrix<double> > superlu_double_colmajor; SuperLU<SparseMatrix<double> > superlu_double_colmajor;
SuperLU<SparseMatrix<std::complex<double> > > superlu_cplxdouble_colmajor; SuperLU<SparseMatrix<std::complex<double> > > superlu_cplxdouble_colmajor;
CALL_SUBTEST_1(sparse_lu(superlu_double_colmajor)); CALL_SUBTEST_1( check_sparse_square_solving(superlu_double_colmajor) );
CALL_SUBTEST_1(sparse_lu(superlu_cplxdouble_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) );
} }
} }

8
unsupported/test/umfpack_support.cpp
View File

@ -22,7 +22,7 @@
// License and a copy of the GNU General Public License along with // License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>. // Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "sparse_lu.h" #include "sparse_solver.h"
#ifdef EIGEN_UMFPACK_SUPPORT #ifdef EIGEN_UMFPACK_SUPPORT
#include <Eigen/UmfPackSupport> #include <Eigen/UmfPackSupport>
@ -33,7 +33,9 @@ void test_umfpack_support()
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
UmfPackLU<SparseMatrix<double> > umfpack_double_colmajor; UmfPackLU<SparseMatrix<double> > umfpack_double_colmajor;
UmfPackLU<SparseMatrix<std::complex<double> > > umfpack_cplxdouble_colmajor; UmfPackLU<SparseMatrix<std::complex<double> > > umfpack_cplxdouble_colmajor;
CALL_SUBTEST_1(sparse_lu(umfpack_double_colmajor)); CALL_SUBTEST_1(check_sparse_square_solving(umfpack_double_colmajor));
CALL_SUBTEST_1(sparse_lu(umfpack_cplxdouble_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));
} }
} }