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