mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

modify the unit tests of sparse linear solvers to enable tests on real matrices, from MatrixMarket for instance

Browse Source
This commit is contained in:
Desire NUENTSA 2012-03-29 14:32:54 +02:00
parent ada9e79145
commit f804a319c8
14 changed files with 177 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists.txt
View File

@ -338,6 +338,13 @@ if(EIGEN_BUILD_BTL)
add_subdirectory(bench/btl EXCLUDE_FROM_ALL)
endif(EIGEN_BUILD_BTL)
if(TEST_REAL_CASES)
if(NOT WIN32)
add_subdirectory(bench/spbench EXCLUDE_FROM_ALL)
set(ENV(EIGEN_MATRIX_DIR) ${TEST_REAL_CASES})
endif(NOT WIN32)
endif(TEST_REAL_CASES)
ei_testing_print_summary()
message(STATUS "")

16
Eigen/src/PaStiXSupport/PaStiXSupport.h
View File

@ -136,10 +136,10 @@ namespace internal
// WARNING It is assumed here that successive calls to this routine are done
// with matrices having the same pattern.
template <typename MatrixType>
void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans)
void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans, bool& hasTranspose)
{
eigen_assert(In.cols()==In.rows() && " Can only symmetrize the graph of a square matrix");
if (StrMatTrans.outerSize() == 0)
if (!hasTranspose)
{ //First call to this routine, need to compute the structural pattern of In^T
StrMatTrans = In.transpose();
// Set the elements of the matrix to zero
@ -148,6 +148,7 @@ namespace internal
for (typename MatrixType::InnerIterator it(StrMatTrans, i); it; ++it)
it.valueRef() = 0.0;
}
hasTranspose = true;
}
Out = (StrMatTrans + In).eval();
}
@ -172,6 +173,7 @@ class PastixBase
PastixBase():m_initisOk(false),m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false)
{
m_pastixdata = 0;
m_hasTranspose = false;
PastixInit();
}
@ -314,7 +316,6 @@ class PastixBase
m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
}
Derived& compute (MatrixType& mat);
@ -328,6 +329,7 @@ class PastixBase
mutable SparseMatrix<Scalar, ColMajor> m_mat_null; // An input null matrix
mutable Matrix<Scalar, Dynamic,1> m_vec_null; // An input null vector
mutable SparseMatrix<Scalar, ColMajor> m_StrMatTrans; // The transpose pattern of the input matrix
mutable bool m_hasTranspose; // The transpose of the current matrix has already been computed
mutable int m_comm; // The MPI communicator identifier
mutable Matrix<Index,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
@ -357,6 +359,7 @@ void PastixBase<Derived>::PastixInit()
PastixDestroy();
m_pastixdata = 0;
m_iparm(IPARM_MODIFY_PARAMETER) = API_YES;
m_hasTranspose = false;
}
m_iparm(IPARM_START_TASK) = API_TASK_INIT;
@ -537,7 +540,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
temp = matrix;
else
{
internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans);
internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans, m_hasTranspose);
}
m_iparm[IPARM_SYM] = API_SYM_NO;
m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
@ -559,7 +562,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
temp = matrix;
else
{
internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans);
internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
}
m_iparm(IPARM_SYM) = API_SYM_NO;
@ -581,7 +584,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
temp = matrix;
else
{
internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans);
internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
}
m_iparm(IPARM_SYM) = API_SYM_NO;
m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
@ -591,6 +594,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
using Base::m_iparm;
using Base::m_dparm;
using Base::m_StrMatTrans;
using Base::m_hasTranspose;
};
/** \ingroup PaStiXSupport_Module

8
Eigen/src/UmfPackSupport/UmfPackSupport.h
View File

@ -124,9 +124,11 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N
* \brief A sparse LU factorization and solver based on UmfPack
*
* This class allows to solve for A.X = B sparse linear problems via a LU factorization
* using the UmfPack library. The sparse matrix A must be column-major, squared and full rank.
* using the UmfPack library. The sparse matrix A must be in a compressed column-major form, squared and full rank.
* The vectors or matrices X and B can be either dense or sparse.
*
* WARNING The Eigen column-major SparseMatrix is not always in compressed form.
* The user should call makeCompressed() to get a matrix in CSC suitable for UMFPACK
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
*
* \sa \ref TutorialSparseDirectSolvers
@ -198,7 +200,9 @@ class UmfPackLU
return m_q;
}
/** Computes the sparse Cholesky decomposition of \a matrix */
/** Computes the sparse Cholesky decomposition of \a matrix
* Note that the matrix should be in compressed format. Please, use makeCompressed() to get it !!
*/
void compute(const MatrixType& matrix)
{
analyzePattern(matrix);

9
test/CMakeLists.txt
View File

@ -18,6 +18,14 @@ set(LAPACK_FOUND TRUE)
set(BLAS_LIBRARIES eigen_blas)
set(LAPACK_LIBRARIES eigen_lapack)
if(TEST_REAL_CASES)
if(NOT WIN32)
add_definitions( -DEIGEN_MATRIXDIR="${TEST_REAL_CASES}" )
else(NOT WIN32)
message(STATUS, "REAL CASES CAN NOT BE CURRENTLY TESTED ON WIN32")
endif(NOT WIN32)
endif(TEST_REAL_CASES)
set(SPARSE_LIBS " ")
find_package(Cholmod)
@ -192,6 +200,7 @@ ei_add_test(vectorwiseop)
ei_add_test(simplicial_cholesky)
ei_add_test(conjugate_gradient)
ei_add_test(bicgstab)
ei_add_test(gmres)
if(UMFPACK_FOUND)

6
test/bicgstab.cpp
View File

@ -40,8 +40,6 @@ template<typename T> void test_bicgstab_T()
void test_bicgstab()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(test_bicgstab_T<double>());
CALL_SUBTEST_2(test_bicgstab_T<std::complex<double> >());
}
CALL_SUBTEST_1(test_bicgstab_T<double>());
CALL_SUBTEST_2(test_bicgstab_T<std::complex<double> >());
}

18
test/cholmod_support.cpp
View File

@ -42,18 +42,16 @@ template<typename T> void test_cholmod_T()
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);
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_cholmod_support()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(test_cholmod_T<double>());
CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >());
}
CALL_SUBTEST_1(test_cholmod_T<double>());
CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >());
}

6
test/conjugate_gradient.cpp
View File

@ -40,8 +40,6 @@ template<typename T> void test_conjugate_gradient_T()
void test_conjugate_gradient()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(test_conjugate_gradient_T<double>());
CALL_SUBTEST_2(test_conjugate_gradient_T<std::complex<double> >());
}
CALL_SUBTEST_1(test_conjugate_gradient_T<double>());
CALL_SUBTEST_2(test_conjugate_gradient_T<std::complex<double> >());
}

10
test/pardiso_support.cpp
View File

@ -22,10 +22,8 @@ template<typename T> void test_pardiso_T()
void test_pardiso_support()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(test_pardiso_T<float>());
CALL_SUBTEST_2(test_pardiso_T<double>());
CALL_SUBTEST_3(test_pardiso_T< std::complex<float> >());
CALL_SUBTEST_4(test_pardiso_T< std::complex<double> >());
}
CALL_SUBTEST_1(test_pardiso_T<float>());
CALL_SUBTEST_2(test_pardiso_T<double>());
CALL_SUBTEST_3(test_pardiso_T< std::complex<float> >());
CALL_SUBTEST_4(test_pardiso_T< std::complex<double> >());
}

10
test/pastix_support.cpp
View File

@ -52,10 +52,8 @@ template<typename T> void test_pastix_T_LU()
void test_pastix_support()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(test_pastix_T<float>());
CALL_SUBTEST_2(test_pastix_T<double>());
CALL_SUBTEST_3( (test_pastix_T_LU<std::complex<float> >()) );
CALL_SUBTEST_4(test_pastix_T_LU<std::complex<double> >());
}
CALL_SUBTEST_1(test_pastix_T<float>());
CALL_SUBTEST_2(test_pastix_T<double>());
CALL_SUBTEST_3( (test_pastix_T_LU<std::complex<float> >()) );
CALL_SUBTEST_4(test_pastix_T_LU<std::complex<double> >());
}

6
test/simplicial_cholesky.cpp
View File

@ -50,8 +50,6 @@ template<typename T> void test_simplicial_cholesky_T()
void test_simplicial_cholesky()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(test_simplicial_cholesky_T<double>());
CALL_SUBTEST_2(test_simplicial_cholesky_T<std::complex<double> >());
}
CALL_SUBTEST_1(test_simplicial_cholesky_T<double>());
CALL_SUBTEST_2(test_simplicial_cholesky_T<std::complex<double> >());
}

1
test/sparse.h
View File

@ -193,4 +193,5 @@ initSparse(double density,
}
}
#include <unsupported/Eigen/SparseExtra>
#endif // EIGEN_TESTSPARSE_H

123
test/sparse_solver.h
View File

@ -74,6 +74,56 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A,
VERIFY(x.isApprox(refX,test_precision<Scalar>()));
}
template<typename Scalar>
inline std::string get_matrixfolder()
{
std::string mat_folder = EIGEN_MATRIXDIR;
if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value )
mat_folder = mat_folder + static_cast<string>("/complex/");
else
mat_folder = mat_folder + static_cast<string>("/real/");
return mat_folder;
}
template<typename Solver, typename Rhs>
void check_sparse_solving_real_cases(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const Rhs& refX)
{
typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef typename Mat::RealScalar RealScalar;
Rhs x(b.rows(), b.cols());
solver.compute(A);
if (solver.info() != Success)
{
std::cerr << "sparse solver testing: factorization failed (check_sparse_solving_real_cases)\n";
exit(0);
return;
}
x = solver.solve(b);
if (solver.info() != Success)
{
std::cerr << "sparse solver testing: solving failed\n";
return;
}
RealScalar res_error;
// Compute the norm of the relative error
if(refX.size() != 0)
res_error = (refX - x).norm()/refX.norm();
else
{
// Compute the relative residual norm
res_error = (b - A * x).norm()/b.norm();
}
if (res_error > test_precision<Scalar>() ){
std::cerr << "Test " << g_test_stack.back() << " failed in "EI_PP_MAKE_STRING(__FILE__)
<< " (" << EI_PP_MAKE_STRING(__LINE__) << ")" << std::endl << std::endl;
abort();
}
}
template<typename Solver, typename DenseMat>
void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
{
@ -121,6 +171,7 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver)
{
typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar;
typedef typename Mat::Index Index;
typedef SparseMatrix<Scalar,ColMajor> SpMat;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
@ -137,13 +188,37 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver)
DenseVector b = DenseVector::Random(size);
DenseMatrix dB(size,rhsCols);
initSparse<Scalar>(density, dB, B);
for (int i = 0; i < g_repeat; i++) {
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);
}
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);
// First, get the folder
#ifdef EIGEN_MATRIXDIR
if (internal::is_same<Scalar, float>::value
|| internal::is_same<Scalar, std::complex<float> >::value)
return ;
std::string mat_folder = get_matrixfolder<Scalar>();
MatrixMarketIterator<Scalar> it(mat_folder);
for (; it; ++it)
{
if (it.sym() == SPD){
Mat halfA;
PermutationMatrix<Dynamic, Dynamic, Index> pnull;
halfA.template selfadjointView<Solver::UpLo>() = it.matrix().template triangularView<Eigen::Lower>().twistedBy(pnull);
std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n";
check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX());
check_sparse_solving_real_cases(solver, halfA, it.rhs(), it.refX());
}
}
#endif
}
template<typename Solver> void check_sparse_spd_determinant(Solver& solver)
@ -156,9 +231,11 @@ template<typename Solver> void check_sparse_spd_determinant(Solver& solver)
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 );
for (int i = 0; i < g_repeat; i++) {
check_sparse_determinant(solver, A, dA);
check_sparse_determinant(solver, halfA, dA );
}
}
template<typename Solver, typename DenseMat>
@ -194,9 +271,27 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver)
DenseVector b = DenseVector::Random(size);
DenseMatrix dB = DenseMatrix::Random(size,rhsCols);
A.makeCompressed();
for (int i = 0; i < g_repeat; i++) {
check_sparse_solving(solver, A, b, dA, b);
check_sparse_solving(solver, A, dB, dA, dB);
}
// First, get the folder
#ifdef EIGEN_MATRIXDIR
if (internal::is_same<Scalar, float>::value
|| internal::is_same<Scalar, std::complex<float> >::value)
return ;
std::string mat_folder = get_matrixfolder<Scalar>();
MatrixMarketIterator<Scalar> it(mat_folder);
for (; it; ++it)
{
std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n";
check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX());
}
#endif
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)
@ -209,6 +304,8 @@ template<typename Solver> void check_sparse_square_determinant(Solver& solver)
Mat A;
DenseMatrix dA;
generate_sparse_square_problem(solver, A, dA, 30);
check_sparse_determinant(solver, A, dA);
A.makeCompressed();
for (int i = 0; i < g_repeat; i++) {
check_sparse_determinant(solver, A, dA);
}
}

14
test/superlu_support.cpp
View File

@ -28,12 +28,10 @@
void test_superlu_support()
{
for(int i = 0; i < g_repeat; i++) {
SuperLU<SparseMatrix<double> > superlu_double_colmajor;
SuperLU<SparseMatrix<std::complex<double> > > 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) );
}
SuperLU<SparseMatrix<double> > superlu_double_colmajor;
SuperLU<SparseMatrix<std::complex<double> > > 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) );
}

14
test/umfpack_support.cpp
View File

@ -28,12 +28,10 @@
void test_umfpack_support()
{
for(int i = 0; i < g_repeat; i++) {
UmfPackLU<SparseMatrix<double> > umfpack_double_colmajor;
UmfPackLU<SparseMatrix<std::complex<double> > > 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));
}
UmfPackLU<SparseMatrix<double, ColMajor> > umfpack_double_colmajor;
UmfPackLU<SparseMatrix<std::complex<double> > > 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));
}