cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

spbenchsolver.h (18167B)


      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 
     11 #include <iostream>
     12 #include <fstream>
     13 #include <Eigen/SparseCore>
     14 #include <bench/BenchTimer.h>
     15 #include <cstdlib>
     16 #include <string>
     17 #include <Eigen/Cholesky>
     18 #include <Eigen/Jacobi>
     19 #include <Eigen/Householder>
     20 #include <Eigen/IterativeLinearSolvers>
     21 #include <unsupported/Eigen/IterativeSolvers>
     22 #include <Eigen/LU>
     23 #include <unsupported/Eigen/SparseExtra>
     24 #include <Eigen/SparseLU>
     25 
     26 #include "spbenchstyle.h"
     27 
     28 #ifdef EIGEN_METIS_SUPPORT
     29 #include <Eigen/MetisSupport>
     30 #endif
     31 
     32 #ifdef EIGEN_CHOLMOD_SUPPORT
     33 #include <Eigen/CholmodSupport>
     34 #endif
     35 
     36 #ifdef EIGEN_UMFPACK_SUPPORT
     37 #include <Eigen/UmfPackSupport>
     38 #endif
     39 
     40 #ifdef EIGEN_KLU_SUPPORT
     41 #include <Eigen/KLUSupport>
     42 #endif
     43 
     44 #ifdef EIGEN_PARDISO_SUPPORT
     45 #include <Eigen/PardisoSupport>
     46 #endif
     47 
     48 #ifdef EIGEN_SUPERLU_SUPPORT
     49 #include <Eigen/SuperLUSupport>
     50 #endif
     51 
     52 #ifdef EIGEN_PASTIX_SUPPORT
     53 #include <Eigen/PaStiXSupport>
     54 #endif
     55 
     56 // CONSTANTS
     57 #define EIGEN_UMFPACK  10
     58 #define EIGEN_KLU  11
     59 #define EIGEN_SUPERLU  20
     60 #define EIGEN_PASTIX  30
     61 #define EIGEN_PARDISO  40
     62 #define EIGEN_SPARSELU_COLAMD 50
     63 #define EIGEN_SPARSELU_METIS 51
     64 #define EIGEN_BICGSTAB  60
     65 #define EIGEN_BICGSTAB_ILUT  61
     66 #define EIGEN_GMRES 70
     67 #define EIGEN_GMRES_ILUT 71
     68 #define EIGEN_SIMPLICIAL_LDLT  80
     69 #define EIGEN_CHOLMOD_LDLT  90
     70 #define EIGEN_PASTIX_LDLT  100
     71 #define EIGEN_PARDISO_LDLT  110
     72 #define EIGEN_SIMPLICIAL_LLT  120
     73 #define EIGEN_CHOLMOD_SUPERNODAL_LLT  130
     74 #define EIGEN_CHOLMOD_SIMPLICIAL_LLT  140
     75 #define EIGEN_PASTIX_LLT  150
     76 #define EIGEN_PARDISO_LLT  160
     77 #define EIGEN_CG  170
     78 #define EIGEN_CG_PRECOND  180
     79 
     80 using namespace Eigen;
     81 using namespace std; 
     82 
     83 
     84 // Global variables for input parameters
     85 int MaximumIters; // Maximum number of iterations
     86 double RelErr; // Relative error of the computed solution
     87 double best_time_val; // Current best time overall solvers 
     88 int best_time_id; //  id of the best solver for the current system 
     89 
     90 template<typename T> inline typename NumTraits<T>::Real test_precision() { return NumTraits<T>::dummy_precision(); }
     91 template<> inline float test_precision<float>() { return 1e-3f; }                                                             
     92 template<> inline double test_precision<double>() { return 1e-6; }                                                            
     93 template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); }
     94 template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); }
     95 
     96 void printStatheader(std::ofstream& out)
     97 {
     98   // Print XML header
     99   // NOTE It would have been much easier to write these XML documents using external libraries like tinyXML or Xerces-C++.
    100   
    101   out << "<?xml version='1.0' encoding='UTF-8'?> \n";
    102   out << "<?xml-stylesheet type='text/xsl' href='#stylesheet' ?> \n"; 
    103   out << "<!DOCTYPE BENCH  [\n<!ATTLIST xsl:stylesheet\n id\t ID  #REQUIRED>\n]>";
    104   out << "\n\n<!-- Generated by the Eigen library -->\n"; 
    105   
    106   out << "\n<BENCH> \n" ; //root XML element 
    107   // Print the xsl style section
    108   printBenchStyle(out); 
    109   // List all available solvers 
    110   out << " <AVAILSOLVER> \n";
    111 #ifdef EIGEN_UMFPACK_SUPPORT
    112   out <<"  <SOLVER ID='" << EIGEN_UMFPACK << "'>\n"; 
    113   out << "   <TYPE> LU </TYPE> \n";
    114   out << "   <PACKAGE> UMFPACK </PACKAGE> \n"; 
    115   out << "  </SOLVER> \n"; 
    116 #endif
    117 #ifdef EIGEN_KLU_SUPPORT
    118   out <<"  <SOLVER ID='" << EIGEN_KLU << "'>\n";
    119   out << "   <TYPE> LU </TYPE> \n";
    120   out << "   <PACKAGE> KLU </PACKAGE> \n";
    121   out << "  </SOLVER> \n";
    122 #endif
    123 #ifdef EIGEN_SUPERLU_SUPPORT
    124   out <<"  <SOLVER ID='" << EIGEN_SUPERLU << "'>\n"; 
    125   out << "   <TYPE> LU </TYPE> \n";
    126   out << "   <PACKAGE> SUPERLU </PACKAGE> \n"; 
    127   out << "  </SOLVER> \n"; 
    128 #endif
    129 #ifdef EIGEN_CHOLMOD_SUPPORT
    130   out <<"  <SOLVER ID='" << EIGEN_CHOLMOD_SIMPLICIAL_LLT << "'>\n"; 
    131   out << "   <TYPE> LLT SP</TYPE> \n";
    132   out << "   <PACKAGE> CHOLMOD </PACKAGE> \n";
    133   out << "  </SOLVER> \n"; 
    134   
    135   out <<"  <SOLVER ID='" << EIGEN_CHOLMOD_SUPERNODAL_LLT << "'>\n"; 
    136   out << "   <TYPE> LLT</TYPE> \n";
    137   out << "   <PACKAGE> CHOLMOD </PACKAGE> \n";
    138   out << "  </SOLVER> \n";
    139   
    140   out <<"  <SOLVER ID='" << EIGEN_CHOLMOD_LDLT << "'>\n"; 
    141   out << "   <TYPE> LDLT </TYPE> \n";
    142   out << "   <PACKAGE> CHOLMOD </PACKAGE> \n";  
    143   out << "  </SOLVER> \n"; 
    144 #endif
    145 #ifdef EIGEN_PARDISO_SUPPORT
    146   out <<"  <SOLVER ID='" << EIGEN_PARDISO << "'>\n"; 
    147   out << "   <TYPE> LU </TYPE> \n";
    148   out << "   <PACKAGE> PARDISO </PACKAGE> \n"; 
    149   out << "  </SOLVER> \n"; 
    150   
    151   out <<"  <SOLVER ID='" << EIGEN_PARDISO_LLT << "'>\n"; 
    152   out << "   <TYPE> LLT </TYPE> \n";
    153   out << "   <PACKAGE> PARDISO </PACKAGE> \n"; 
    154   out << "  </SOLVER> \n"; 
    155   
    156   out <<"  <SOLVER ID='" << EIGEN_PARDISO_LDLT << "'>\n"; 
    157   out << "   <TYPE> LDLT </TYPE> \n";
    158   out << "   <PACKAGE> PARDISO </PACKAGE> \n"; 
    159   out << "  </SOLVER> \n"; 
    160 #endif
    161 #ifdef EIGEN_PASTIX_SUPPORT
    162   out <<"  <SOLVER ID='" << EIGEN_PASTIX << "'>\n"; 
    163   out << "   <TYPE> LU </TYPE> \n";
    164   out << "   <PACKAGE> PASTIX </PACKAGE> \n"; 
    165   out << "  </SOLVER> \n"; 
    166   
    167   out <<"  <SOLVER ID='" << EIGEN_PASTIX_LLT << "'>\n"; 
    168   out << "   <TYPE> LLT </TYPE> \n";
    169   out << "   <PACKAGE> PASTIX </PACKAGE> \n"; 
    170   out << "  </SOLVER> \n"; 
    171   
    172   out <<"  <SOLVER ID='" << EIGEN_PASTIX_LDLT << "'>\n"; 
    173   out << "   <TYPE> LDLT </TYPE> \n";
    174   out << "   <PACKAGE> PASTIX </PACKAGE> \n"; 
    175   out << "  </SOLVER> \n"; 
    176 #endif
    177   
    178   out <<"  <SOLVER ID='" << EIGEN_BICGSTAB << "'>\n"; 
    179   out << "   <TYPE> BICGSTAB </TYPE> \n";
    180   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    181   out << "  </SOLVER> \n"; 
    182   
    183   out <<"  <SOLVER ID='" << EIGEN_BICGSTAB_ILUT << "'>\n"; 
    184   out << "   <TYPE> BICGSTAB_ILUT </TYPE> \n";
    185   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    186   out << "  </SOLVER> \n"; 
    187   
    188   out <<"  <SOLVER ID='" << EIGEN_GMRES_ILUT << "'>\n"; 
    189   out << "   <TYPE> GMRES_ILUT </TYPE> \n";
    190   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    191   out << "  </SOLVER> \n"; 
    192   
    193   out <<"  <SOLVER ID='" << EIGEN_SIMPLICIAL_LDLT << "'>\n"; 
    194   out << "   <TYPE> LDLT </TYPE> \n";
    195   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    196   out << "  </SOLVER> \n"; 
    197   
    198   out <<"  <SOLVER ID='" << EIGEN_SIMPLICIAL_LLT << "'>\n"; 
    199   out << "   <TYPE> LLT </TYPE> \n";
    200   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    201   out << "  </SOLVER> \n"; 
    202   
    203   out <<"  <SOLVER ID='" << EIGEN_CG << "'>\n"; 
    204   out << "   <TYPE> CG </TYPE> \n";
    205   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    206   out << "  </SOLVER> \n"; 
    207   
    208   out <<"  <SOLVER ID='" << EIGEN_SPARSELU_COLAMD << "'>\n"; 
    209   out << "   <TYPE> LU_COLAMD </TYPE> \n";
    210   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    211   out << "  </SOLVER> \n"; 
    212   
    213 #ifdef EIGEN_METIS_SUPPORT
    214   out <<"  <SOLVER ID='" << EIGEN_SPARSELU_METIS << "'>\n"; 
    215   out << "   <TYPE> LU_METIS </TYPE> \n";
    216   out << "   <PACKAGE> EIGEN </PACKAGE> \n"; 
    217   out << "  </SOLVER> \n"; 
    218 #endif
    219   out << " </AVAILSOLVER> \n"; 
    220   
    221 }
    222 
    223 
    224 template<typename Solver, typename Scalar>
    225 void call_solver(Solver &solver, const int solver_id, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX,std::ofstream& statbuf)
    226 {
    227   
    228   double total_time;
    229   double compute_time;
    230   double solve_time; 
    231   double rel_error;
    232   Matrix<Scalar, Dynamic, 1> x; 
    233   BenchTimer timer; 
    234   timer.reset();
    235   timer.start();
    236   solver.compute(A); 
    237   if (solver.info() != Success)
    238   {
    239     std::cerr << "Solver failed ... \n";
    240     return;
    241   }
    242   timer.stop();
    243   compute_time = timer.value();
    244   statbuf << "    <TIME>\n"; 
    245   statbuf << "     <COMPUTE> " << timer.value() << "</COMPUTE>\n";
    246   std::cout<< "COMPUTE TIME : " << timer.value() <<std::endl; 
    247     
    248   timer.reset();
    249   timer.start();
    250   x = solver.solve(b); 
    251   if (solver.info() == NumericalIssue)
    252   {
    253     std::cerr << "Solver failed ... \n";
    254     return;
    255   }
    256   timer.stop();
    257   solve_time = timer.value();
    258   statbuf << "     <SOLVE> " << timer.value() << "</SOLVE>\n"; 
    259   std::cout<< "SOLVE TIME : " << timer.value() <<std::endl; 
    260   
    261   total_time = solve_time + compute_time;
    262   statbuf << "     <TOTAL> " << total_time << "</TOTAL>\n"; 
    263   std::cout<< "TOTAL TIME : " << total_time <<std::endl; 
    264   statbuf << "    </TIME>\n"; 
    265   
    266   // Verify the relative error
    267   if(refX.size() != 0)
    268     rel_error = (refX - x).norm()/refX.norm();
    269   else 
    270   {
    271     // Compute the relative residual norm
    272     Matrix<Scalar, Dynamic, 1> temp; 
    273     temp = A * x; 
    274     rel_error = (b-temp).norm()/b.norm();
    275   }
    276   statbuf << "    <ERROR> " << rel_error << "</ERROR>\n"; 
    277   std::cout<< "REL. ERROR : " << rel_error << "\n\n" ;
    278   if ( rel_error <= RelErr )
    279   {
    280     // check the best time if convergence
    281     if(!best_time_val || (best_time_val > total_time))
    282     {
    283       best_time_val = total_time;
    284       best_time_id = solver_id;
    285     }
    286   }
    287 }
    288 
    289 template<typename Solver, typename Scalar>
    290 void call_directsolver(Solver& solver, const int solver_id, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX, std::string& statFile)
    291 {
    292     std::ofstream statbuf(statFile.c_str(), std::ios::app);
    293     statbuf << "   <SOLVER_STAT ID='" << solver_id <<"'>\n"; 
    294     call_solver(solver, solver_id, A, b, refX,statbuf);
    295     statbuf << "   </SOLVER_STAT>\n";
    296     statbuf.close();
    297 }
    298 
    299 template<typename Solver, typename Scalar>
    300 void call_itersolver(Solver &solver, const int solver_id, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX, std::string& statFile)
    301 {
    302   solver.setTolerance(RelErr); 
    303   solver.setMaxIterations(MaximumIters);
    304   
    305   std::ofstream statbuf(statFile.c_str(), std::ios::app);
    306   statbuf << " <SOLVER_STAT ID='" << solver_id <<"'>\n"; 
    307   call_solver(solver, solver_id, A, b, refX,statbuf); 
    308   statbuf << "   <ITER> "<< solver.iterations() << "</ITER>\n";
    309   statbuf << " </SOLVER_STAT>\n";
    310   std::cout << "ITERATIONS : " << solver.iterations() <<"\n\n\n"; 
    311   
    312 }
    313 
    314 
    315 template <typename Scalar>
    316 void SelectSolvers(const SparseMatrix<Scalar>&A, unsigned int sym, Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX, std::string& statFile)
    317 {
    318   typedef SparseMatrix<Scalar, ColMajor> SpMat; 
    319   // First, deal with Nonsymmetric and symmetric matrices
    320   best_time_id = 0; 
    321   best_time_val = 0.0;
    322   //UMFPACK
    323   #ifdef EIGEN_UMFPACK_SUPPORT
    324   {
    325     cout << "Solving with UMFPACK LU ... \n"; 
    326     UmfPackLU<SpMat> solver; 
    327     call_directsolver(solver, EIGEN_UMFPACK, A, b, refX,statFile); 
    328   }
    329   #endif
    330   //KLU
    331   #ifdef EIGEN_KLU_SUPPORT
    332   {
    333     cout << "Solving with KLU LU ... \n";
    334     KLU<SpMat> solver;
    335     call_directsolver(solver, EIGEN_KLU, A, b, refX,statFile);
    336   }
    337   #endif
    338     //SuperLU
    339   #ifdef EIGEN_SUPERLU_SUPPORT
    340   {
    341     cout << "\nSolving with SUPERLU ... \n"; 
    342     SuperLU<SpMat> solver;
    343     call_directsolver(solver, EIGEN_SUPERLU, A, b, refX,statFile); 
    344   }
    345   #endif
    346     
    347    // PaStix LU
    348   #ifdef EIGEN_PASTIX_SUPPORT
    349   {
    350     cout << "\nSolving with PASTIX LU ... \n"; 
    351     PastixLU<SpMat> solver; 
    352     call_directsolver(solver, EIGEN_PASTIX, A, b, refX,statFile) ;
    353   }
    354   #endif
    355 
    356    //PARDISO LU
    357   #ifdef EIGEN_PARDISO_SUPPORT
    358   {
    359     cout << "\nSolving with PARDISO LU ... \n"; 
    360     PardisoLU<SpMat>  solver; 
    361     call_directsolver(solver, EIGEN_PARDISO, A, b, refX,statFile);
    362   }
    363   #endif
    364   
    365   // Eigen SparseLU METIS
    366   cout << "\n Solving with Sparse LU AND COLAMD ... \n";
    367   SparseLU<SpMat, COLAMDOrdering<int> >   solver;
    368   call_directsolver(solver, EIGEN_SPARSELU_COLAMD, A, b, refX, statFile); 
    369   // Eigen SparseLU METIS
    370   #ifdef EIGEN_METIS_SUPPORT
    371   {
    372     cout << "\n Solving with Sparse LU AND METIS ... \n";
    373     SparseLU<SpMat, MetisOrdering<int> >   solver;
    374     call_directsolver(solver, EIGEN_SPARSELU_METIS, A, b, refX, statFile); 
    375   }
    376   #endif
    377   
    378   //BiCGSTAB
    379   {
    380     cout << "\nSolving with BiCGSTAB ... \n"; 
    381     BiCGSTAB<SpMat> solver; 
    382     call_itersolver(solver, EIGEN_BICGSTAB, A, b, refX,statFile);
    383   }
    384   //BiCGSTAB+ILUT
    385   {
    386     cout << "\nSolving with BiCGSTAB and ILUT ... \n"; 
    387     BiCGSTAB<SpMat, IncompleteLUT<Scalar> > solver; 
    388     call_itersolver(solver, EIGEN_BICGSTAB_ILUT, A, b, refX,statFile); 
    389   }
    390   
    391    
    392   //GMRES
    393 //   {
    394 //     cout << "\nSolving with GMRES ... \n"; 
    395 //     GMRES<SpMat> solver; 
    396 //     call_itersolver(solver, EIGEN_GMRES, A, b, refX,statFile); 
    397 //   }
    398   //GMRES+ILUT
    399   {
    400     cout << "\nSolving with GMRES and ILUT ... \n"; 
    401     GMRES<SpMat, IncompleteLUT<Scalar> > solver; 
    402     call_itersolver(solver, EIGEN_GMRES_ILUT, A, b, refX,statFile);
    403   }
    404   
    405   // Hermitian and not necessarily positive-definites
    406   if (sym != NonSymmetric)
    407   {
    408     // Internal Cholesky
    409     {
    410       cout << "\nSolving with Simplicial LDLT ... \n"; 
    411       SimplicialLDLT<SpMat, Lower> solver;
    412       call_directsolver(solver, EIGEN_SIMPLICIAL_LDLT, A, b, refX,statFile); 
    413     }
    414     
    415     // CHOLMOD
    416     #ifdef EIGEN_CHOLMOD_SUPPORT
    417     {
    418       cout << "\nSolving with CHOLMOD LDLT ... \n"; 
    419       CholmodDecomposition<SpMat, Lower> solver;
    420       solver.setMode(CholmodLDLt);
    421        call_directsolver(solver,EIGEN_CHOLMOD_LDLT, A, b, refX,statFile);
    422     }
    423     #endif
    424     
    425     //PASTIX LLT
    426     #ifdef EIGEN_PASTIX_SUPPORT
    427     {
    428       cout << "\nSolving with PASTIX LDLT ... \n"; 
    429       PastixLDLT<SpMat, Lower> solver; 
    430       call_directsolver(solver,EIGEN_PASTIX_LDLT, A, b, refX,statFile); 
    431     }
    432     #endif
    433     
    434     //PARDISO LLT
    435     #ifdef EIGEN_PARDISO_SUPPORT
    436     {
    437       cout << "\nSolving with PARDISO LDLT ... \n"; 
    438       PardisoLDLT<SpMat, Lower> solver; 
    439       call_directsolver(solver,EIGEN_PARDISO_LDLT, A, b, refX,statFile); 
    440     }
    441     #endif
    442   }
    443 
    444    // Now, symmetric POSITIVE DEFINITE matrices
    445   if (sym == SPD)
    446   {
    447     
    448     //Internal Sparse Cholesky
    449     {
    450       cout << "\nSolving with SIMPLICIAL LLT ... \n"; 
    451       SimplicialLLT<SpMat, Lower> solver; 
    452       call_directsolver(solver,EIGEN_SIMPLICIAL_LLT, A, b, refX,statFile); 
    453     }
    454     
    455     // CHOLMOD
    456     #ifdef EIGEN_CHOLMOD_SUPPORT
    457     {
    458       // CholMOD SuperNodal LLT
    459       cout << "\nSolving with CHOLMOD LLT (Supernodal)... \n"; 
    460       CholmodDecomposition<SpMat, Lower> solver;
    461       solver.setMode(CholmodSupernodalLLt);
    462        call_directsolver(solver,EIGEN_CHOLMOD_SUPERNODAL_LLT, A, b, refX,statFile);
    463       // CholMod Simplicial LLT
    464       cout << "\nSolving with CHOLMOD LLT (Simplicial) ... \n"; 
    465       solver.setMode(CholmodSimplicialLLt);
    466       call_directsolver(solver,EIGEN_CHOLMOD_SIMPLICIAL_LLT, A, b, refX,statFile);
    467     }
    468     #endif
    469     
    470     //PASTIX LLT
    471     #ifdef EIGEN_PASTIX_SUPPORT
    472     {
    473       cout << "\nSolving with PASTIX LLT ... \n"; 
    474       PastixLLT<SpMat, Lower> solver; 
    475       call_directsolver(solver,EIGEN_PASTIX_LLT, A, b, refX,statFile);
    476     }
    477     #endif
    478     
    479     //PARDISO LLT
    480     #ifdef EIGEN_PARDISO_SUPPORT
    481     {
    482       cout << "\nSolving with PARDISO LLT ... \n"; 
    483       PardisoLLT<SpMat, Lower> solver; 
    484       call_directsolver(solver,EIGEN_PARDISO_LLT, A, b, refX,statFile); 
    485     }
    486     #endif
    487     
    488     // Internal CG
    489     {
    490       cout << "\nSolving with CG ... \n"; 
    491       ConjugateGradient<SpMat, Lower> solver; 
    492       call_itersolver(solver,EIGEN_CG, A, b, refX,statFile);
    493     }
    494     //CG+IdentityPreconditioner
    495 //     {
    496 //       cout << "\nSolving with CG and IdentityPreconditioner ... \n"; 
    497 //       ConjugateGradient<SpMat, Lower, IdentityPreconditioner> solver; 
    498 //       call_itersolver(solver,EIGEN_CG_PRECOND, A, b, refX,statFile);
    499 //     }
    500   } // End SPD matrices 
    501 }
    502 
    503 /* Browse all the matrices available in the specified folder 
    504  * and solve the associated linear system.
    505  * The results of each solve are printed in the standard output
    506  * and optionally in the provided html file
    507  */
    508 template <typename Scalar>
    509 void Browse_Matrices(const string folder, bool statFileExists, std::string& statFile, int maxiters, double tol)
    510 {
    511   MaximumIters = maxiters; // Maximum number of iterations, global variable 
    512   RelErr = tol;  //Relative residual error  as stopping criterion for iterative solvers
    513   MatrixMarketIterator<Scalar> it(folder);
    514   for ( ; it; ++it)
    515   {
    516     //print the infos for this linear system 
    517     if(statFileExists)
    518     {
    519       std::ofstream statbuf(statFile.c_str(), std::ios::app);
    520       statbuf << "<LINEARSYSTEM> \n";
    521       statbuf << "   <MATRIX> \n";
    522       statbuf << "     <NAME> " << it.matname() << " </NAME>\n"; 
    523       statbuf << "     <SIZE> " << it.matrix().rows() << " </SIZE>\n"; 
    524       statbuf << "     <ENTRIES> " << it.matrix().nonZeros() << "</ENTRIES>\n";
    525       if (it.sym()!=NonSymmetric)
    526       {
    527         statbuf << "     <SYMMETRY> Symmetric </SYMMETRY>\n" ; 
    528         if (it.sym() == SPD) 
    529           statbuf << "     <POSDEF> YES </POSDEF>\n"; 
    530         else 
    531           statbuf << "     <POSDEF> NO </POSDEF>\n"; 
    532           
    533       }
    534       else
    535       {
    536         statbuf << "     <SYMMETRY> NonSymmetric </SYMMETRY>\n" ; 
    537         statbuf << "     <POSDEF> NO </POSDEF>\n"; 
    538       }
    539       statbuf << "   </MATRIX> \n";
    540       statbuf.close();
    541     }
    542     
    543     cout<< "\n\n===================================================== \n";
    544     cout<< " ======  SOLVING WITH MATRIX " << it.matname() << " ====\n";
    545     cout<< " =================================================== \n\n";
    546     Matrix<Scalar, Dynamic, 1> refX;
    547     if(it.hasrefX()) refX = it.refX();
    548     // Call all suitable solvers for this linear system 
    549     SelectSolvers<Scalar>(it.matrix(), it.sym(), it.rhs(), refX, statFile);
    550     
    551     if(statFileExists)
    552     {
    553       std::ofstream statbuf(statFile.c_str(), std::ios::app);
    554       statbuf << "  <BEST_SOLVER ID='"<< best_time_id
    555               << "'></BEST_SOLVER>\n"; 
    556       statbuf << " </LINEARSYSTEM> \n"; 
    557       statbuf.close();
    558     }
    559   } 
    560 } 
    561 
    562 bool get_options(int argc, char **args, string option, string* value=0)
    563 {
    564   int idx = 1, found=false; 
    565   while (idx<argc && !found){
    566     if (option.compare(args[idx]) == 0){
    567       found = true; 
    568       if(value) *value = args[idx+1];
    569     }
    570     idx+=2;
    571   }
    572   return found; 
    573 }