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 } 555