clean sparse LU tests

This commit is contained in:
Gael Guennebaud 2011-09-24 17:15:37 +02:00
parent b2988375e8
commit 9bba0e7ba1
5 changed files with 182 additions and 29 deletions

View File

@ -25,6 +25,7 @@ if(UMFPACK_FOUND AND BLAS_FOUND)
add_definitions("-DEIGEN_UMFPACK_SUPPORT")
include_directories(${UMFPACK_INCLUDES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${UMFPACK_LIBRARIES} ${BLAS_LIBRARIES})
set(UMFPACK_ALL_LIBS ${UMFPACK_LIBRARIES} ${BLAS_LIBRARIES})
ei_add_property(EIGEN_TESTED_BACKENDS "UmfPack, ")
else()
ei_add_property(EIGEN_MISSING_BACKENDS "UmfPack, ")
@ -35,6 +36,7 @@ if(SUPERLU_FOUND AND BLAS_FOUND)
add_definitions("-DEIGEN_SUPERLU_SUPPORT")
include_directories(${SUPERLU_INCLUDES})
set(SPARSE_LIBS ${SPARSE_LIBS} ${SUPERLU_LIBRARIES} ${BLAS_LIBRARIES})
set(SUPERLU_ALL_LIBS ${SUPERLU_LIBRARIES} ${BLAS_LIBRARIES})
ei_add_property(EIGEN_TESTED_BACKENDS "SuperLU, ")
else()
ei_add_property(EIGEN_MISSING_BACKENDS "SuperLU, ")
@ -88,7 +90,16 @@ endif()
ei_add_test(sparse_llt "" "${SPARSE_LIBS}")
ei_add_test(sparse_ldlt "" "${SPARSE_LIBS}")
ei_add_test(sparse_lu "" "${SPARSE_LIBS}")
ei_add_test(sparse_lu_legacy "" "${SPARSE_LIBS}")
if(UMFPACK_FOUND)
ei_add_test(umfpack_support "" "${UMFPACK_ALL_LIBS}")
endif()
if(SUPERLU_FOUND)
ei_add_test(superlu_support "" "${SUPERLU_ALL_LIBS}")
endif()
ei_add_test(sparse_extra "" "")
find_package(FFTW)

View File

@ -0,0 +1,86 @@
// 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);
}

View File

@ -33,7 +33,8 @@
#include <Eigen/SuperLUSupport>
#endif
template<typename Scalar> void sparse_lu(int rows, int cols)
template<typename Scalar> void sparse_lu_legacy(int rows, int cols)
{
double density = (std::max)(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
@ -110,39 +111,16 @@ template<typename Scalar> void sparse_lu(int rows, int cols)
else
std::cerr << "super lu factorize failed\n";
}
// New API
{
x.setZero();
SuperLU<SparseMatrix<Scalar> > slu(m2);
if (slu.info() == Success)
{
DenseVector oldb = b;
x = slu.solve(b);
VERIFY(oldb.isApprox(b) && "the rhs should not be modified!");
if (slu.info() == Success) {
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "SuperLU");
}
else
std::cerr << "super lu solving failed\n";
if (!NumTraits<Scalar>::IsComplex) {
VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
}
}
else
std::cerr << "super lu factorize failed\n";
}
#endif
}
void test_sparse_lu()
void test_sparse_lu_legacy()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(sparse_lu<double>(8, 8) );
CALL_SUBTEST_1(sparse_lu_legacy<double>(8, 8) );
int s = internal::random<int>(1,300);
CALL_SUBTEST_2(sparse_lu<std::complex<double> >(s,s) );
CALL_SUBTEST_1(sparse_lu<double>(s,s) );
CALL_SUBTEST_1(sparse_lu_legacy<std::complex<double> >(s,s) );
CALL_SUBTEST_1(sparse_lu_legacy<double>(s,s) );
}
}

View File

@ -0,0 +1,39 @@
// 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_lu.h"
#ifdef EIGEN_SUPERLU_SUPPORT
#include <Eigen/SuperLUSupport>
#endif
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(sparse_lu(superlu_double_colmajor));
CALL_SUBTEST_1(sparse_lu(superlu_cplxdouble_colmajor));
}
}

View File

@ -0,0 +1,39 @@
// 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_lu.h"
#ifdef EIGEN_UMFPACK_SUPPORT
#include <Eigen/UmfPackSupport>
#endif
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(sparse_lu(umfpack_double_colmajor));
CALL_SUBTEST_1(sparse_lu(umfpack_cplxdouble_colmajor));
}
}