diff options
| -rw-r--r-- | bench/spbench/CMakeLists.txt | 65 | ||||
| -rw-r--r-- | bench/spbench/spbenchsolver.cpp | 79 | ||||
| -rw-r--r-- | bench/spbench/spbenchsolver.h | 539 | ||||
| -rw-r--r-- | unsupported/Eigen/SparseExtra | 2 | ||||
| -rw-r--r-- | unsupported/Eigen/src/SparseExtra/MarketIO.h | 11 | ||||
| -rw-r--r-- | unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h | 235 |
6 files changed, 926 insertions, 5 deletions
diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt new file mode 100644 index 000000000..079912266 --- /dev/null +++ b/bench/spbench/CMakeLists.txt @@ -0,0 +1,65 @@ + + +set(BLAS_FOUND TRUE) +set(LAPACK_FOUND TRUE) +set(BLAS_LIBRARIES eigen_blas_static) +set(LAPACK_LIBRARIES eigen_lapack_static) + +set(SPARSE_LIBS "") + +# find_library(PARDISO_LIBRARIES pardiso412-GNU450-X86-64) +# if(PARDISO_LIBRARIES) +# add_definitions("-DEIGEN_PARDISO_SUPPORT") +# set(SPARSE_LIBS ${SPARSE_LIBS} ${PARDISO_LIBRARIES}) +# endif(PARDISO_LIBRARIES) + +find_package(Cholmod) +if(CHOLMOD_FOUND AND BLAS_FOUND AND LAPACK_FOUND) + add_definitions("-DEIGEN_CHOLMOD_SUPPORT") + include_directories(${CHOLMOD_INCLUDES}) + set(SPARSE_LIBS ${SPARSE_LIBS} ${CHOLMOD_LIBRARIES} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) + set(CHOLMOD_ALL_LIBS ${CHOLMOD_LIBRARIES} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) +endif() + +find_package(Umfpack) +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}) +endif() + +find_package(SuperLU) +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}) +endif() + + +find_package(Pastix) +find_package(Scotch) +find_package(Metis) +if(PASTIX_FOUND AND BLAS_FOUND) + add_definitions("-DEIGEN_PASTIX_SUPPORT") + include_directories(${PASTIX_INCLUDES}) + if(SCOTCH_FOUND) + include_directories(${SCOTCH_INCLUDES}) + set(PASTIX_LIBRARIES ${PASTIX_LIBRARIES} ${SCOTCH_LIBRARIES}) + elseif(METIS_FOUND) + include_directories(${METIS_INCLUDES}) + set(PASTIX_LIBRARIES ${PASTIX_LIBRARIES} ${METIS_LIBRARIES}) + endif(SCOTCH_FOUND) + set(SPARSE_LIBS ${SPARSE_LIBS} ${PASTIX_LIBRARIES} ${ORDERING_LIBRARIES} ${BLAS_LIBRARIES}) + set(PASTIX_ALL_LIBS ${PASTIX_LIBRARIES} ${BLAS_LIBRARIES}) +endif(PASTIX_FOUND AND BLAS_FOUND) + +find_library(RT_LIBRARY rt) +if(RT_LIBRARY) + set(SPARSE_LIBS ${SPARSE_LIBS} ${RT_LIBRARY}) +endif(RT_LIBRARY) + +add_executable(spbenchsolver spbenchsolver.cpp) +target_link_libraries (spbenchsolver ${SPARSE_LIBS}) + diff --git a/bench/spbench/spbenchsolver.cpp b/bench/spbench/spbenchsolver.cpp new file mode 100644 index 000000000..bcbe0d84f --- /dev/null +++ b/bench/spbench/spbenchsolver.cpp @@ -0,0 +1,79 @@ +#include <bench/spbench/spbenchsolver.h> + +void bench_printhelp() +{ + cout<< " \nbenchsolver : performs a benchmark of all the solvers available in Eigen \n\n"; + cout<< " MATRIX FOLDER : \n"; + cout<< " The matrices for the benchmark should be collected in a folder specified with an environment variable EIGEN_MATRIXDIR \n"; + cout<< " This folder should contain the subfolders real/ and complex/ : \n"; + cout<< " The matrices are stored using the matrix market coordinate format \n"; + cout<< " The matrix and associated right-hand side (rhs) files are named respectively \n"; + cout<< " as MatrixName.mtx and MatrixName_b.mtx. If the rhs does not exist, a random one is generated. \n"; + cout<< " If a matrix is SPD, the matrix should be named as MatrixName_SPD.mtx \n"; + cout<< " If a true solution exists, it should be named as MatrixName_x.mtx; \n" ; + cout<< " it will be used to compute the norm of the error relative to the computed solutions\n\n"; + cout<< " OPTIONS : \n"; + cout<< " -h or --help \n print this help and return\n\n"; + cout<< " -d matrixdir \n Use matrixdir as the matrix folder instead of the one specified in the environment variable EIGEN_MATRIXDIR\n\n"; + cout<< " -o outputfile.html \n Output the statistics to a html file \n\n"; + +} +int main(int argc, char ** args) +{ + + bool help = ( get_options(argc, args, "-h") || get_options(argc, args, "--help") ); + if(help) { + bench_printhelp(); + return 0; + } + + // Get the location of the test matrices + string matrix_dir; + if (!get_options(argc, args, "-d", &matrix_dir)) + { + if(getenv("EIGEN_MATRIXDIR") == NULL){ + std::cerr << "Please, specify the location of the matrices with -d mat_folder or the environment variable EIGEN_MATRIXDIR \n"; + std::cerr << " Run with --help to see the list of all the available options \n"; + return -1; + } + matrix_dir = getenv("EIGEN_MATRIXDIR"); + } + + std::ofstream statbuf; + string statFile ; + + // Get the file to write the statistics + bool statFileExists = get_options(argc, args, "-o", &statFile); + if(statFileExists) + { + statbuf.open(statFile.c_str(), std::ios::out); + if(statbuf.good()){ + statFileExists = true; + printStatheader(statbuf); + statbuf.close(); + } + else + std::cerr << "Unable to open the provided file for writting... \n"; + } + + string current_dir; + // Test the matrices in %EIGEN_MATRIXDIR/real + current_dir = matrix_dir + "/real"; + Browse_Matrices<double>(current_dir, statFileExists, statFile); + + // Test the matrices in %EIGEN_MATRIXDIR/complex + current_dir = matrix_dir + "/complex"; + Browse_Matrices<std::complex<double> >(current_dir, statFileExists, statFile); + + if(statFileExists) + { + statbuf.open(statFile.c_str(), std::ios::app); + statbuf << "</TABLE> \n"; + cout << "\n Output written in " << statFile << " ...\n"; + statbuf.close(); + } + + return 0; +} + + diff --git a/bench/spbench/spbenchsolver.h b/bench/spbench/spbenchsolver.h new file mode 100644 index 000000000..2dffb3b4b --- /dev/null +++ b/bench/spbench/spbenchsolver.h @@ -0,0 +1,539 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.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 <iostream> +#include <fstream> +#include "Eigen/SparseCore" +#include <bench/BenchTimer.h> +#include <cstdlib> +#include <string> +#include <Eigen/Cholesky> +#include <Eigen/Jacobi> +#include <Eigen/Householder> +#include <Eigen/IterativeLinearSolvers> +#include <Eigen/LU> +#include <unsupported/Eigen/SparseExtra> + +#ifdef EIGEN_CHOLMOD_SUPPORT +#include <Eigen/CholmodSupport> +#endif + +#ifdef EIGEN_UMFPACK_SUPPORT +#include <Eigen/UmfPackSupport> +#endif + +#ifdef EIGEN_PARDISO_SUPPORT +#include <Eigen/PardisoSupport> +#endif + +#ifdef EIGEN_SUPERLU_SUPPORT +#include <Eigen/SuperLUSupport> +#endif + +#ifdef EIGEN_PASTIX_SUPPORT +#include <Eigen/PaStiXSupport> +#endif + +// CONSTANTS +#define EIGEN_UMFPACK 0 +#define EIGEN_SUPERLU 1 +#define EIGEN_PASTIX 2 +#define EIGEN_PARDISO 3 +#define EIGEN_BICGSTAB 4 +#define EIGEN_BICGSTAB_ILUT 5 +#define EIGEN_GMRES 6 +#define EIGEN_GMRES_ILUT 7 +#define EIGEN_SIMPLICIAL_LDLT 8 +#define EIGEN_CHOLMOD_LDLT 9 +#define EIGEN_PASTIX_LDLT 10 +#define EIGEN_PARDISO_LDLT 11 +#define EIGEN_SIMPLICIAL_LLT 12 +#define EIGEN_CHOLMOD_SUPERNODAL_LLT 13 +#define EIGEN_CHOLMOD_SIMPLICIAL_LLT 14 +#define EIGEN_PASTIX_LLT 15 +#define EIGEN_PARDISO_LLT 16 +#define EIGEN_CG 17 +#define EIGEN_CG_PRECOND 18 +#define EIGEN_ALL_SOLVERS 19 + +using namespace Eigen; +using namespace std; + +struct Stats{ + ComputationInfo info; + double total_time; + double compute_time; + double solve_time; + double rel_error; + int memory_used; + int iterations; + int isavail; + int isIterative; +}; + +template<typename T> inline typename NumTraits<T>::Real test_precision() { return NumTraits<T>::dummy_precision(); } +template<> inline float test_precision<float>() { return 1e-3f; } +template<> inline double test_precision<double>() { return 1e-6; } +template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); } +template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); } + +void printStatheader(std::ofstream& out) +{ + int LUcnt = 0; + string LUlist =" ", LLTlist = "<TH > LLT", LDLTlist = "<TH > LDLT "; + +#ifdef EIGEN_UMFPACK_SUPPORT + LUlist += "<TH > UMFPACK "; LUcnt++; +#endif +#ifdef EIGEN_SUPERLU_SUPPORT + LUlist += "<TH > SUPERLU "; LUcnt++; +#endif +#ifdef EIGEN_CHOLMOD_SUPPORT + LLTlist += "<TH > CHOLMOD SP LLT<TH > CHOLMOD LLT"; + LDLTlist += "<TH>CHOLMOD LDLT"; +#endif +#ifdef EIGEN_PARDISO_SUPPORT + LUlist += "<TH > PARDISO LU"; LUcnt++; + LLTlist += "<TH > PARDISO LLT"; + LDLTlist += "<TH > PARDISO LDLT"; +#endif +#ifdef EIGEN_PASTIX_SUPPORT + LUlist += "<TH > PASTIX LU"; LUcnt++; + LLTlist += "<TH > PASTIX LLT"; + LDLTlist += "<TH > PASTIX LDLT"; +#endif + + out << "<TABLE border=\"1\" >\n "; + out << "<TR><TH>Matrix <TH> N <TH> NNZ <TH> "; + if (LUcnt) out << LUlist; + out << " <TH >BiCGSTAB <TH >BiCGSTAB+ILUT"<< "<TH >GMRES+ILUT" <<LDLTlist << LLTlist << "<TH> CG "<< std::endl; +} + + +template<typename Solver, typename Scalar> +Stats call_solver(Solver &solver, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX) +{ + Stats stat; + Matrix<Scalar, Dynamic, 1> x; + BenchTimer timer; + timer.reset(); + timer.start(); + solver.compute(A); + if (solver.info() != Success) + { + stat.info = NumericalIssue; + std::cerr << "Solver failed ... \n"; + return stat; + } + timer.stop(); + stat.compute_time = timer.value(); + + timer.reset(); + timer.start(); + x = solver.solve(b); + if (solver.info() == NumericalIssue) + { + stat.info = NumericalIssue; + std::cerr << "Solver failed ... \n"; + return stat; + } + + timer.stop(); + stat.solve_time = timer.value(); + stat.total_time = stat.solve_time + stat.compute_time; + stat.memory_used = 0; + // Verify the relative error + if(refX.size() != 0) + stat.rel_error = (refX - x).norm()/refX.norm(); + else + { + // Compute the relative residual norm + Matrix<Scalar, Dynamic, 1> temp; + temp = A * x; + stat.rel_error = (b-temp).norm()/b.norm(); + } + if ( stat.rel_error > test_precision<Scalar>() ) + { + stat.info = NoConvergence; + return stat; + } + else + { + stat.info = Success; + return stat; + } +} + +template<typename Solver, typename Scalar> +Stats call_directsolver(Solver& solver, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX) +{ + Stats stat; + stat = call_solver(solver, A, b, refX); + return stat; +} + +template<typename Solver, typename Scalar> +Stats call_itersolver(Solver &solver, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX) +{ + Stats stat; + solver.setTolerance(1e-10); + stat = call_solver(solver, A, b, refX); + stat.iterations = solver.iterations(); + return stat; +} + +inline void printStatItem(Stats *stat, int solver_id, int& best_time_id, double& best_time_val) +{ + stat[solver_id].isavail = 1; + + if (stat[solver_id].info == NumericalIssue) + { + cout << " SOLVER FAILED ... Probably a numerical issue \n"; + return; + } + if (stat[solver_id].info == NoConvergence){ + cout << "REL. ERROR " << stat[solver_id].rel_error; + if(stat[solver_id].isIterative == 1) + cout << " (" << stat[solver_id].iterations << ") \n"; + return; + } + + // Record the best CPU time + if (!best_time_val) + { + best_time_val = stat[solver_id].total_time; + best_time_id = solver_id; + } + else if (stat[solver_id].total_time < best_time_val) + { + best_time_val = stat[solver_id].total_time; + best_time_id = solver_id; + } + // Print statistics to standard output + if (stat[solver_id].info == Success){ + cout<< "COMPUTE TIME : " << stat[solver_id].compute_time<< " \n"; + cout<< "SOLVE TIME : " << stat[solver_id].solve_time<< " \n"; + cout<< "TOTAL TIME : " << stat[solver_id].total_time<< " \n"; + cout << "REL. ERROR : " << stat[solver_id].rel_error ; + if(stat[solver_id].isIterative == 1) { + cout << " (" << stat[solver_id].iterations << ") "; + } + cout << std::endl; + } + +} + + +/* Print the results from all solvers corresponding to a particular matrix + * The best CPU time is printed in bold + */ +inline void printHtmlStatLine(Stats *stat, int best_time_id, string& statline) +{ + + string markup; + ostringstream compute,solve,total,error; + for (int i = 0; i < EIGEN_ALL_SOLVERS; i++) + { + if (stat[i].isavail == 0) continue; + if(i == best_time_id) + markup = "<TD style=\"background-color:red\">"; + else + markup = "<TD>"; + + if (stat[i].info == Success){ + compute << markup << stat[i].compute_time; + solve << markup << stat[i].solve_time; + total << markup << stat[i].total_time; + error << " <TD> " << stat[i].rel_error; + if(stat[i].isIterative == 1) { + error << " (" << stat[i].iterations << ") "; + } + } + else { + compute << " <TD> -" ; + solve << " <TD> -" ; + total << " <TD> -" ; + if(stat[i].info == NoConvergence){ + error << " <TD> "<< stat[i].rel_error ; + if(stat[i].isIterative == 1) + error << " (" << stat[i].iterations << ") "; + } + else error << " <TD> - "; + } + } + + statline = "<TH>Compute Time " + compute.str() + "\n" + + "<TR><TH>Solve Time " + solve.str() + "\n" + + "<TR><TH>Total Time " + total.str() + "\n" + +"<TR><TH>Error(Iter)" + error.str() + "\n"; + +} + +template <typename Scalar> +int SelectSolvers(const SparseMatrix<Scalar>&A, unsigned int sym, Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX, Stats *stat) +{ + typedef SparseMatrix<Scalar, ColMajor> SpMat; + // First, deal with Nonsymmetric and symmetric matrices + int best_time_id = 0; + double best_time_val = 0.0; + //UMFPACK + #ifdef EIGEN_UMFPACK_SUPPORT + { + cout << "Solving with UMFPACK LU ... \n"; + UmfPackLU<SpMat> solver; + stat[EIGEN_UMFPACK] = call_directsolver(solver, A, b, refX); + printStatItem(stat, EIGEN_UMFPACK, best_time_id, best_time_val); + } + #endif + //SuperLU + #ifdef EIGEN_SUPERLU_SUPPORT + { + cout << "\nSolving with SUPERLU ... \n"; + SuperLU<SpMat> solver; + stat[EIGEN_SUPERLU] = call_directsolver(solver, A, b, refX); + printStatItem(stat, EIGEN_SUPERLU, best_time_id, best_time_val); + } + #endif + + // PaStix LU + #ifdef EIGEN_PASTIX_SUPPORT + { + cout << "\nSolving with PASTIX LU ... \n"; + PastixLU<SpMat> solver; + stat[EIGEN_PASTIX] = call_directsolver(solver, A, b, refX) ; + printStatItem(stat, EIGEN_PASTIX, best_time_id, best_time_val); + } + #endif + + //PARDISO LU + #ifdef EIGEN_PARDISO_SUPPORT + { + cout << "\nSolving with PARDISO LU ... \n"; + PardisoLU<SpMat> solver; + stat[EIGEN_PARDISO] = call_directsolver(solver, A, b, refX); + printStatItem(stat, EIGEN_PARDISO, best_time_id, best_time_val); + } + #endif + + + + //BiCGSTAB + { + cout << "\nSolving with BiCGSTAB ... \n"; + BiCGSTAB<SpMat> solver; + stat[EIGEN_BICGSTAB] = call_itersolver(solver, A, b, refX); + stat[EIGEN_BICGSTAB].isIterative = 1; + printStatItem(stat, EIGEN_BICGSTAB, best_time_id, best_time_val); + } + //BiCGSTAB+ILUT + { + cout << "\nSolving with BiCGSTAB and ILUT ... \n"; + BiCGSTAB<SpMat, IncompleteLUT<Scalar> > solver; + stat[EIGEN_BICGSTAB_ILUT] = call_itersolver(solver, A, b, refX); + stat[EIGEN_BICGSTAB_ILUT].isIterative = 1; + printStatItem(stat, EIGEN_BICGSTAB_ILUT, best_time_id, best_time_val); + } + + + //GMRES +// { +// cout << "\nSolving with GMRES ... \n"; +// GMRES<SpMat> solver; +// stat[EIGEN_GMRES] = call_itersolver(solver, A, b, refX); +// stat[EIGEN_GMRES].isIterative = 1; +// printStatItem(stat, EIGEN_GMRES, best_time_id, best_time_val); +// } + //GMRES+ILUT + { + cout << "\nSolving with GMRES and ILUT ... \n"; + GMRES<SpMat, IncompleteLUT<Scalar> > solver; + stat[EIGEN_GMRES_ILUT] = call_itersolver(solver, A, b, refX); + stat[EIGEN_GMRES_ILUT].isIterative = 1; + printStatItem(stat, EIGEN_GMRES_ILUT, best_time_id, best_time_val); + } + + // Symmetric and not necessarily positive-definites + if ( (sym == Symmetric) || (sym == SPD) ) + { + // Internal Cholesky + { + cout << "\nSolving with Simplicial LDLT ... \n"; + SimplicialLDLT<SpMat, Lower> solver; + stat[EIGEN_SIMPLICIAL_LDLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat, EIGEN_SIMPLICIAL_LDLT, best_time_id, best_time_val); + } + + // CHOLMOD + #ifdef EIGEN_CHOLMOD_SUPPORT + { + cout << "\nSolving with CHOLMOD LDLT ... \n"; + CholmodDecomposition<SpMat, Lower> solver; + solver.setMode(CholmodLDLt); + stat[EIGEN_CHOLMOD_LDLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_CHOLMOD_LDLT, best_time_id, best_time_val); + } + #endif + + //PASTIX LLT + #ifdef EIGEN_PASTIX_SUPPORT + { + cout << "\nSolving with PASTIX LDLT ... \n"; + PastixLDLT<SpMat, Lower> solver; + stat[EIGEN_PASTIX_LDLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_PASTIX_LDLT, best_time_id, best_time_val); + } + #endif + + //PARDISO LLT + #ifdef EIGEN_PARDISO_SUPPORT + { + cout << "\nSolving with PARDISO LDLT ... \n"; + PardisoLDLT<SpMat, Lower> solver; + stat[EIGEN_PARDISO_LDLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_PARDISO_LDLT, best_time_id, best_time_val); + } + #endif + } + + // Now, symmetric POSITIVE DEFINITE matrices + if (sym == SPD) + { + + //Internal Sparse Cholesky + { + cout << "\nSolving with SIMPLICIAL LLT ... \n"; + SimplicialLLT<SpMat, Lower> solver; + stat[EIGEN_SIMPLICIAL_LLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_SIMPLICIAL_LLT, best_time_id, best_time_val); + } + + // CHOLMOD + #ifdef EIGEN_CHOLMOD_SUPPORT + { + // CholMOD SuperNodal LLT + cout << "\nSolving with CHOLMOD LLT (Supernodal)... \n"; + CholmodDecomposition<SpMat, Lower> solver; + solver.setMode(CholmodSupernodalLLt); + stat[EIGEN_CHOLMOD_SUPERNODAL_LLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_CHOLMOD_SUPERNODAL_LLT, best_time_id, best_time_val); + // CholMod Simplicial LLT + cout << "\nSolving with CHOLMOD LLT (Simplicial) ... \n"; + solver.setMode(CholmodSimplicialLLt); + stat[EIGEN_CHOLMOD_SIMPLICIAL_LLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_CHOLMOD_SIMPLICIAL_LLT, best_time_id, best_time_val); + } + #endif + + //PASTIX LLT + #ifdef EIGEN_PASTIX_SUPPORT + { + cout << "\nSolving with PASTIX LLT ... \n"; + PastixLLT<SpMat, Lower> solver; + stat[EIGEN_PASTIX_LLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_PASTIX_LLT, best_time_id, best_time_val); + } + #endif + + //PARDISO LLT + #ifdef EIGEN_PARDISO_SUPPORT + { + cout << "\nSolving with PARDISO LLT ... \n"; + PardisoLLT<SpMat, Lower> solver; + stat[EIGEN_PARDISO_LLT] = call_directsolver(solver, A, b, refX); + printStatItem(stat,EIGEN_PARDISO_LLT, best_time_id, best_time_val); + } + #endif + + // Internal CG + { + cout << "\nSolving with CG ... \n"; + ConjugateGradient<SpMat, Lower> solver; + stat[EIGEN_CG] = call_itersolver(solver, A, b, refX); + stat[EIGEN_CG].isIterative = 1; + printStatItem(stat,EIGEN_CG, best_time_id, best_time_val); + } + //CG+IdentityPreconditioner +// { +// cout << "\nSolving with CG and IdentityPreconditioner ... \n"; +// ConjugateGradient<SpMat, Lower, IdentityPreconditioner> solver; +// stat[EIGEN_CG_PRECOND] = call_itersolver(solver, A, b, refX); +// stat[EIGEN_CG_PRECOND].isIterative = 1; +// printStatItem(stat,EIGEN_CG_PRECOND, best_time_id, best_time_val); +// } + } // End SPD matrices + + return best_time_id; +} + +/* Browse all the matrices available in the specified folder + * and solve the associated linear system. + * The results of each solve are printed in the standard output + * and optionally in the provided html file + */ +template <typename Scalar> +void Browse_Matrices(const string folder, bool statFileExists, std::string& statFile) +{ + MatrixMarketIterator<Scalar> it(folder); + Stats stat[EIGEN_ALL_SOLVERS]; + for ( ; it; ++it) + { + for (int i = 0; i < EIGEN_ALL_SOLVERS; i++) + { + stat[i].isavail = 0; + stat[i].isIterative = 0; + } + + int best_time_id; + cout<< "\n\n===================================================== \n"; + cout<< " ====== SOLVING WITH MATRIX " << it.matname() << " ====\n"; + cout<< " =================================================== \n\n"; + Matrix<Scalar, Dynamic, 1> refX; + if(it.hasrefX()) refX = it.refX(); + best_time_id = SelectSolvers<Scalar>(it.matrix(), it.sym(), it.rhs(), refX, &stat[0]); + + if(statFileExists) + { + string statline; + printHtmlStatLine(&stat[0], best_time_id, statline); + std::ofstream statbuf(statFile.c_str(), std::ios::app); + statbuf << "<TR><TH rowspan=\"4\">" << it.matname() << " <TD rowspan=\"4\"> " + << it.matrix().rows() << " <TD rowspan=\"4\"> " << it.matrix().nonZeros()<< " "<< statline ; + statbuf.close(); + } + } +} + +bool get_options(int argc, char **args, string option, string* value=0) +{ + int idx = 1, found=false; + while (idx<argc && !found){ + if (option.compare(args[idx]) == 0){ + found = true; + if(value) *value = args[idx+1]; + } + idx+=2; + } + return found; +} diff --git a/unsupported/Eigen/SparseExtra b/unsupported/Eigen/SparseExtra index a09bfb7e5..f73a085be 100644 --- a/unsupported/Eigen/SparseExtra +++ b/unsupported/Eigen/SparseExtra @@ -38,7 +38,7 @@ namespace Eigen { #include "src/SparseExtra/RandomSetter.h" #include "src/SparseExtra/MarketIO.h" - +#include "src/SparseExtra/MatrixMarketIterator.h" } // namespace Eigen #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h" diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h index 4e2d5dc4e..3b103914a 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -166,10 +166,13 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename) if(!readsizes) { line >> M >> N >> NNZ; - readsizes = true; - std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; - mat.resize(M,N); - mat.reserve(NNZ); + if(M > 0 && N > 0 && NNZ > 0) + { + readsizes = true; + std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; + mat.resize(M,N); + mat.reserve(NNZ); + } } else { diff --git a/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h new file mode 100644 index 000000000..fe8406ea4 --- /dev/null +++ b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h @@ -0,0 +1,235 @@ + +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 +// +// 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/>. +#ifndef EIGEN_BROWSE_MATRICES_H +#define EIGEN_BROWSE_MATRICES_H + +#include <dirent.h> +#include <unsupported/Eigen/SparseExtra> +using namespace Eigen; +using std::string; + +enum { + SPD = 0x100, + NonSymmetric = 0x0 +}; + +/** + * @brief Iterator to browse matrices from a specified folder + * + * This is used to load all the matrices from a folder. + * The matrices should be in Matrix Market format + * It is assumed that the matrices are named as matname.mtx + * and matname_SPD.mtx if the matrix is Symmetric and positive definite (or Hermitian) + * The right hand side vectors are loaded as well, if they exist. + * They should be named as matname_b.mtx. + * Note that the right hand side for a SPD matrix is named as matname_SPD_b.mtx + * + * Sometimes a reference solution is available. In this case, it should be named as matname_x.mtx + * + * Sample code + * \code + * + * \endcode + * + * \tparam Scalar The scalar type + */ +template <typename Scalar> +class MatrixMarketIterator +{ + public: + typedef Matrix<Scalar,Dynamic,1> VectorType; + typedef SparseMatrix<Scalar,ColMajor> MatrixType; + + public: + MatrixMarketIterator(const string folder):m_sym(0),m_isvalid(false),m_matIsLoaded(false),m_hasRhs(false),m_hasrefX(false),m_folder(folder) + { + m_folder_id = opendir(folder.c_str()); + if (!m_folder_id){ + m_isvalid = false; + std::cerr << "The provided Matrix folder could not be opened \n\n"; + abort(); + } + Getnextvalidmatrix(); + } + + ~MatrixMarketIterator() + { + if (m_folder_id) closedir(m_folder_id); + } + + inline MatrixMarketIterator& operator++() + { + m_matIsLoaded = false; + m_hasrefX = false; + m_hasRhs = false; + Getnextvalidmatrix(); + return *this; + } + inline operator bool() { return m_isvalid;} + + /** Return the sparse matrix corresponding to the current file */ + inline MatrixType& matrix() + { + // Read the matrix + if (m_matIsLoaded) return m_mat; + + string matrix_file = m_folder + "/" + m_matname + ".mtx"; + if ( !loadMarket(m_mat, matrix_file)) + { + m_matIsLoaded = false; + return m_mat; + } + m_matIsLoaded = true; + + if (m_sym != NonSymmetric) + { // Store the upper part of the matrix. It is needed by the solvers dealing with nonsymmetric matrices ?? + MatrixType B; + B = m_mat; + m_mat = B.template selfadjointView<Lower>(); + } + return m_mat; + } + + /** Return the right hand side corresponding to the current matrix. + * If the rhs file is not provided, a random rhs is generated + */ + inline VectorType& rhs() + { + // Get the right hand side + if (m_hasRhs) return m_rhs; + + string rhs_file; + rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx + m_hasRhs = Fileexists(rhs_file); + if (m_hasRhs) + { + m_rhs.resize(m_mat.cols()); + m_hasRhs = loadMarketVector(m_rhs, rhs_file); + } + if (!m_hasRhs) + { + // Generate a random right hand side + if (!m_matIsLoaded) this->matrix(); + m_refX.resize(m_mat.cols()); + m_refX.setRandom(); + m_rhs = m_mat * m_refX; + m_hasrefX = true; + m_hasRhs = true; + } + return m_rhs; + } + + /** Return a reference solution + * If it is not provided and if the right hand side is not available + * then refX is randomly generated such that A*refX = b + * where A and b are the matrix and the rhs. + * Note that when a rhs is provided, refX is not available + */ + inline VectorType& refX() + { + // Check if a reference solution is provided + if (m_hasrefX) return m_refX; + + string lhs_file; + lhs_file = m_folder + "/" + m_matname + "_x.mtx"; + m_hasrefX = Fileexists(lhs_file); + if (m_hasrefX) + { + m_refX.resize(m_mat.cols()); + m_hasrefX = loadMarketVector(m_refX, lhs_file); + } + return m_refX; + } + + inline string& matname() { return m_matname; } + + inline int sym() { return m_sym; } + + inline bool hasRhs() {return m_hasRhs; } + inline bool hasrefX() {return m_hasrefX; } + + protected: + + inline bool Fileexists(string file) + { + std::ifstream file_id(file.c_str()); + if (!file_id.good() ) + { + return false; + } + else + { + file_id.close(); + return true; + } + } + + void Getnextvalidmatrix( ) + { + // Here, we return with the next valid matrix in the folder + while ( (m_curs_id = readdir(m_folder_id)) != NULL) { + m_isvalid = false; + string curfile; + curfile = m_folder + "/" + m_curs_id->d_name; + // Discard if it is a folder + if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems +// struct stat st_buf; +// stat (curfile.c_str(), &st_buf); +// if (S_ISDIR(st_buf.st_mode)) continue; + + // Determine from the header if it is a matrix or a right hand side + bool isvector,iscomplex; + if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue; + if(isvector) continue; + + // Get the matrix name + string filename = m_curs_id->d_name; + m_matname = filename.substr(0, filename.length()-4); + + // Find if the matrix is SPD + size_t found = m_matname.find("SPD"); + if( (found!=string::npos) && (m_sym == Symmetric) ) + m_sym = SPD; + + m_isvalid = true; + break; + } + } + int m_sym; // Symmetry of the matrix + MatrixType m_mat; // Current matrix + VectorType m_rhs; // Current vector + VectorType m_refX; // The reference solution, if exists + string m_matname; // Matrix Name + bool m_isvalid; + bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file + bool m_hasRhs; // The right hand side exists + bool m_hasrefX; // A reference solution is provided + string m_folder; + DIR * m_folder_id; + struct dirent *m_curs_id; + +}; + +#endif
\ No newline at end of file |