1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) 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 25 #ifdef EIGEN_CHOLMOD_SUPPORT 26 #include <Eigen/CholmodSupport> 27 #endif 28 29 #ifdef EIGEN_UMFPACK_SUPPORT 30 #include <Eigen/UmfPackSupport> 31 #endif 32 33 #ifdef EIGEN_PARDISO_SUPPORT 34 #include <Eigen/PardisoSupport> 35 #endif 36 37 #ifdef EIGEN_SUPERLU_SUPPORT 38 #include <Eigen/SuperLUSupport> 39 #endif 40 41 #ifdef EIGEN_PASTIX_SUPPORT 42 #include <Eigen/PaStiXSupport> 43 #endif 44 45 // CONSTANTS 46 #define EIGEN_UMFPACK 0 47 #define EIGEN_SUPERLU 1 48 #define EIGEN_PASTIX 2 49 #define EIGEN_PARDISO 3 50 #define EIGEN_BICGSTAB 4 51 #define EIGEN_BICGSTAB_ILUT 5 52 #define EIGEN_GMRES 6 53 #define EIGEN_GMRES_ILUT 7 54 #define EIGEN_SIMPLICIAL_LDLT 8 55 #define EIGEN_CHOLMOD_LDLT 9 56 #define EIGEN_PASTIX_LDLT 10 57 #define EIGEN_PARDISO_LDLT 11 58 #define EIGEN_SIMPLICIAL_LLT 12 59 #define EIGEN_CHOLMOD_SUPERNODAL_LLT 13 60 #define EIGEN_CHOLMOD_SIMPLICIAL_LLT 14 61 #define EIGEN_PASTIX_LLT 15 62 #define EIGEN_PARDISO_LLT 16 63 #define EIGEN_CG 17 64 #define EIGEN_CG_PRECOND 18 65 #define EIGEN_ALL_SOLVERS 19 66 67 using namespace Eigen; 68 using namespace std; 69 70 struct Stats{ 71 ComputationInfo info; 72 double total_time; 73 double compute_time; 74 double solve_time; 75 double rel_error; 76 int memory_used; 77 int iterations; 78 int isavail; 79 int isIterative; 80 }; 81 82 // Global variables for input parameters 83 int MaximumIters; // Maximum number of iterations 84 double RelErr; // Relative error of the computed solution 85 86 template<typename T> inline typename NumTraits<T>::Real test_precision() { return NumTraits<T>::dummy_precision(); } 87 template<> inline float test_precision<float>() { return 1e-3f; } 88 template<> inline double test_precision<double>() { return 1e-6; } 89 template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); } 90 template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); } 91 92 void printStatheader(std::ofstream& out) 93 { 94 int LUcnt = 0; 95 string LUlist =" ", LLTlist = "<TH > LLT", LDLTlist = "<TH > LDLT "; 96 97 #ifdef EIGEN_UMFPACK_SUPPORT 98 LUlist += "<TH > UMFPACK "; LUcnt++; 99 #endif 100 #ifdef EIGEN_SUPERLU_SUPPORT 101 LUlist += "<TH > SUPERLU "; LUcnt++; 102 #endif 103 #ifdef EIGEN_CHOLMOD_SUPPORT 104 LLTlist += "<TH > CHOLMOD SP LLT<TH > CHOLMOD LLT"; 105 LDLTlist += "<TH>CHOLMOD LDLT"; 106 #endif 107 #ifdef EIGEN_PARDISO_SUPPORT 108 LUlist += "<TH > PARDISO LU"; LUcnt++; 109 LLTlist += "<TH > PARDISO LLT"; 110 LDLTlist += "<TH > PARDISO LDLT"; 111 #endif 112 #ifdef EIGEN_PASTIX_SUPPORT 113 LUlist += "<TH > PASTIX LU"; LUcnt++; 114 LLTlist += "<TH > PASTIX LLT"; 115 LDLTlist += "<TH > PASTIX LDLT"; 116 #endif 117 118 out << "<TABLE border=\"1\" >\n "; 119 out << "<TR><TH>Matrix <TH> N <TH> NNZ <TH> "; 120 if (LUcnt) out << LUlist; 121 out << " <TH >BiCGSTAB <TH >BiCGSTAB+ILUT"<< "<TH >GMRES+ILUT" <<LDLTlist << LLTlist << "<TH> CG "<< std::endl; 122 } 123 124 125 template<typename Solver, typename Scalar> 126 Stats call_solver(Solver &solver, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX) 127 { 128 Stats stat; 129 Matrix<Scalar, Dynamic, 1> x; 130 BenchTimer timer; 131 timer.reset(); 132 timer.start(); 133 solver.compute(A); 134 if (solver.info() != Success) 135 { 136 stat.info = NumericalIssue; 137 std::cerr << "Solver failed ... \n"; 138 return stat; 139 } 140 timer.stop(); 141 stat.compute_time = timer.value(); 142 143 timer.reset(); 144 timer.start(); 145 x = solver.solve(b); 146 if (solver.info() == NumericalIssue) 147 { 148 stat.info = NumericalIssue; 149 std::cerr << "Solver failed ... \n"; 150 return stat; 151 } 152 153 timer.stop(); 154 stat.solve_time = timer.value(); 155 stat.total_time = stat.solve_time + stat.compute_time; 156 stat.memory_used = 0; 157 // Verify the relative error 158 if(refX.size() != 0) 159 stat.rel_error = (refX - x).norm()/refX.norm(); 160 else 161 { 162 // Compute the relative residual norm 163 Matrix<Scalar, Dynamic, 1> temp; 164 temp = A * x; 165 stat.rel_error = (b-temp).norm()/b.norm(); 166 } 167 if ( stat.rel_error > RelErr ) 168 { 169 stat.info = NoConvergence; 170 return stat; 171 } 172 else 173 { 174 stat.info = Success; 175 return stat; 176 } 177 } 178 179 template<typename Solver, typename Scalar> 180 Stats call_directsolver(Solver& solver, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX) 181 { 182 Stats stat; 183 stat = call_solver(solver, A, b, refX); 184 return stat; 185 } 186 187 template<typename Solver, typename Scalar> 188 Stats call_itersolver(Solver &solver, const typename Solver::MatrixType& A, const Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX) 189 { 190 Stats stat; 191 solver.setTolerance(RelErr); 192 solver.setMaxIterations(MaximumIters); 193 194 stat = call_solver(solver, A, b, refX); 195 stat.iterations = solver.iterations(); 196 return stat; 197 } 198 199 inline void printStatItem(Stats *stat, int solver_id, int& best_time_id, double& best_time_val) 200 { 201 stat[solver_id].isavail = 1; 202 203 if (stat[solver_id].info == NumericalIssue) 204 { 205 cout << " SOLVER FAILED ... Probably a numerical issue \n"; 206 return; 207 } 208 if (stat[solver_id].info == NoConvergence){ 209 cout << "REL. ERROR " << stat[solver_id].rel_error; 210 if(stat[solver_id].isIterative == 1) 211 cout << " (" << stat[solver_id].iterations << ") \n"; 212 return; 213 } 214 215 // Record the best CPU time 216 if (!best_time_val) 217 { 218 best_time_val = stat[solver_id].total_time; 219 best_time_id = solver_id; 220 } 221 else if (stat[solver_id].total_time < best_time_val) 222 { 223 best_time_val = stat[solver_id].total_time; 224 best_time_id = solver_id; 225 } 226 // Print statistics to standard output 227 if (stat[solver_id].info == Success){ 228 cout<< "COMPUTE TIME : " << stat[solver_id].compute_time<< " \n"; 229 cout<< "SOLVE TIME : " << stat[solver_id].solve_time<< " \n"; 230 cout<< "TOTAL TIME : " << stat[solver_id].total_time<< " \n"; 231 cout << "REL. ERROR : " << stat[solver_id].rel_error ; 232 if(stat[solver_id].isIterative == 1) { 233 cout << " (" << stat[solver_id].iterations << ") "; 234 } 235 cout << std::endl; 236 } 237 238 } 239 240 241 /* Print the results from all solvers corresponding to a particular matrix 242 * The best CPU time is printed in bold 243 */ 244 inline void printHtmlStatLine(Stats *stat, int best_time_id, string& statline) 245 { 246 247 string markup; 248 ostringstream compute,solve,total,error; 249 for (int i = 0; i < EIGEN_ALL_SOLVERS; i++) 250 { 251 if (stat[i].isavail == 0) continue; 252 if(i == best_time_id) 253 markup = "<TD style=\"background-color:red\">"; 254 else 255 markup = "<TD>"; 256 257 if (stat[i].info == Success){ 258 compute << markup << stat[i].compute_time; 259 solve << markup << stat[i].solve_time; 260 total << markup << stat[i].total_time; 261 error << " <TD> " << stat[i].rel_error; 262 if(stat[i].isIterative == 1) { 263 error << " (" << stat[i].iterations << ") "; 264 } 265 } 266 else { 267 compute << " <TD> -" ; 268 solve << " <TD> -" ; 269 total << " <TD> -" ; 270 if(stat[i].info == NoConvergence){ 271 error << " <TD> "<< stat[i].rel_error ; 272 if(stat[i].isIterative == 1) 273 error << " (" << stat[i].iterations << ") "; 274 } 275 else error << " <TD> - "; 276 } 277 } 278 279 statline = "<TH>Compute Time " + compute.str() + "\n" 280 + "<TR><TH>Solve Time " + solve.str() + "\n" 281 + "<TR><TH>Total Time " + total.str() + "\n" 282 +"<TR><TH>Error(Iter)" + error.str() + "\n"; 283 284 } 285 286 template <typename Scalar> 287 int SelectSolvers(const SparseMatrix<Scalar>&A, unsigned int sym, Matrix<Scalar, Dynamic, 1>& b, const Matrix<Scalar, Dynamic, 1>& refX, Stats *stat) 288 { 289 typedef SparseMatrix<Scalar, ColMajor> SpMat; 290 // First, deal with Nonsymmetric and symmetric matrices 291 int best_time_id = 0; 292 double best_time_val = 0.0; 293 //UMFPACK 294 #ifdef EIGEN_UMFPACK_SUPPORT 295 { 296 cout << "Solving with UMFPACK LU ... \n"; 297 UmfPackLU<SpMat> solver; 298 stat[EIGEN_UMFPACK] = call_directsolver(solver, A, b, refX); 299 printStatItem(stat, EIGEN_UMFPACK, best_time_id, best_time_val); 300 } 301 #endif 302 //SuperLU 303 #ifdef EIGEN_SUPERLU_SUPPORT 304 { 305 cout << "\nSolving with SUPERLU ... \n"; 306 SuperLU<SpMat> solver; 307 stat[EIGEN_SUPERLU] = call_directsolver(solver, A, b, refX); 308 printStatItem(stat, EIGEN_SUPERLU, best_time_id, best_time_val); 309 } 310 #endif 311 312 // PaStix LU 313 #ifdef EIGEN_PASTIX_SUPPORT 314 { 315 cout << "\nSolving with PASTIX LU ... \n"; 316 PastixLU<SpMat> solver; 317 stat[EIGEN_PASTIX] = call_directsolver(solver, A, b, refX) ; 318 printStatItem(stat, EIGEN_PASTIX, best_time_id, best_time_val); 319 } 320 #endif 321 322 //PARDISO LU 323 #ifdef EIGEN_PARDISO_SUPPORT 324 { 325 cout << "\nSolving with PARDISO LU ... \n"; 326 PardisoLU<SpMat> solver; 327 stat[EIGEN_PARDISO] = call_directsolver(solver, A, b, refX); 328 printStatItem(stat, EIGEN_PARDISO, best_time_id, best_time_val); 329 } 330 #endif 331 332 333 334 //BiCGSTAB 335 { 336 cout << "\nSolving with BiCGSTAB ... \n"; 337 BiCGSTAB<SpMat> solver; 338 stat[EIGEN_BICGSTAB] = call_itersolver(solver, A, b, refX); 339 stat[EIGEN_BICGSTAB].isIterative = 1; 340 printStatItem(stat, EIGEN_BICGSTAB, best_time_id, best_time_val); 341 } 342 //BiCGSTAB+ILUT 343 { 344 cout << "\nSolving with BiCGSTAB and ILUT ... \n"; 345 BiCGSTAB<SpMat, IncompleteLUT<Scalar> > solver; 346 stat[EIGEN_BICGSTAB_ILUT] = call_itersolver(solver, A, b, refX); 347 stat[EIGEN_BICGSTAB_ILUT].isIterative = 1; 348 printStatItem(stat, EIGEN_BICGSTAB_ILUT, best_time_id, best_time_val); 349 } 350 351 352 //GMRES 353 // { 354 // cout << "\nSolving with GMRES ... \n"; 355 // GMRES<SpMat> solver; 356 // stat[EIGEN_GMRES] = call_itersolver(solver, A, b, refX); 357 // stat[EIGEN_GMRES].isIterative = 1; 358 // printStatItem(stat, EIGEN_GMRES, best_time_id, best_time_val); 359 // } 360 //GMRES+ILUT 361 { 362 cout << "\nSolving with GMRES and ILUT ... \n"; 363 GMRES<SpMat, IncompleteLUT<Scalar> > solver; 364 stat[EIGEN_GMRES_ILUT] = call_itersolver(solver, A, b, refX); 365 stat[EIGEN_GMRES_ILUT].isIterative = 1; 366 printStatItem(stat, EIGEN_GMRES_ILUT, best_time_id, best_time_val); 367 } 368 369 // Hermitian and not necessarily positive-definites 370 if (sym != NonSymmetric) 371 { 372 // Internal Cholesky 373 { 374 cout << "\nSolving with Simplicial LDLT ... \n"; 375 SimplicialLDLT<SpMat, Lower> solver; 376 stat[EIGEN_SIMPLICIAL_LDLT] = call_directsolver(solver, A, b, refX); 377 printStatItem(stat, EIGEN_SIMPLICIAL_LDLT, best_time_id, best_time_val); 378 } 379 380 // CHOLMOD 381 #ifdef EIGEN_CHOLMOD_SUPPORT 382 { 383 cout << "\nSolving with CHOLMOD LDLT ... \n"; 384 CholmodDecomposition<SpMat, Lower> solver; 385 solver.setMode(CholmodLDLt); 386 stat[EIGEN_CHOLMOD_LDLT] = call_directsolver(solver, A, b, refX); 387 printStatItem(stat,EIGEN_CHOLMOD_LDLT, best_time_id, best_time_val); 388 } 389 #endif 390 391 //PASTIX LLT 392 #ifdef EIGEN_PASTIX_SUPPORT 393 { 394 cout << "\nSolving with PASTIX LDLT ... \n"; 395 PastixLDLT<SpMat, Lower> solver; 396 stat[EIGEN_PASTIX_LDLT] = call_directsolver(solver, A, b, refX); 397 printStatItem(stat,EIGEN_PASTIX_LDLT, best_time_id, best_time_val); 398 } 399 #endif 400 401 //PARDISO LLT 402 #ifdef EIGEN_PARDISO_SUPPORT 403 { 404 cout << "\nSolving with PARDISO LDLT ... \n"; 405 PardisoLDLT<SpMat, Lower> solver; 406 stat[EIGEN_PARDISO_LDLT] = call_directsolver(solver, A, b, refX); 407 printStatItem(stat,EIGEN_PARDISO_LDLT, best_time_id, best_time_val); 408 } 409 #endif 410 } 411 412 // Now, symmetric POSITIVE DEFINITE matrices 413 if (sym == SPD) 414 { 415 416 //Internal Sparse Cholesky 417 { 418 cout << "\nSolving with SIMPLICIAL LLT ... \n"; 419 SimplicialLLT<SpMat, Lower> solver; 420 stat[EIGEN_SIMPLICIAL_LLT] = call_directsolver(solver, A, b, refX); 421 printStatItem(stat,EIGEN_SIMPLICIAL_LLT, best_time_id, best_time_val); 422 } 423 424 // CHOLMOD 425 #ifdef EIGEN_CHOLMOD_SUPPORT 426 { 427 // CholMOD SuperNodal LLT 428 cout << "\nSolving with CHOLMOD LLT (Supernodal)... \n"; 429 CholmodDecomposition<SpMat, Lower> solver; 430 solver.setMode(CholmodSupernodalLLt); 431 stat[EIGEN_CHOLMOD_SUPERNODAL_LLT] = call_directsolver(solver, A, b, refX); 432 printStatItem(stat,EIGEN_CHOLMOD_SUPERNODAL_LLT, best_time_id, best_time_val); 433 // CholMod Simplicial LLT 434 cout << "\nSolving with CHOLMOD LLT (Simplicial) ... \n"; 435 solver.setMode(CholmodSimplicialLLt); 436 stat[EIGEN_CHOLMOD_SIMPLICIAL_LLT] = call_directsolver(solver, A, b, refX); 437 printStatItem(stat,EIGEN_CHOLMOD_SIMPLICIAL_LLT, best_time_id, best_time_val); 438 } 439 #endif 440 441 //PASTIX LLT 442 #ifdef EIGEN_PASTIX_SUPPORT 443 { 444 cout << "\nSolving with PASTIX LLT ... \n"; 445 PastixLLT<SpMat, Lower> solver; 446 stat[EIGEN_PASTIX_LLT] = call_directsolver(solver, A, b, refX); 447 printStatItem(stat,EIGEN_PASTIX_LLT, best_time_id, best_time_val); 448 } 449 #endif 450 451 //PARDISO LLT 452 #ifdef EIGEN_PARDISO_SUPPORT 453 { 454 cout << "\nSolving with PARDISO LLT ... \n"; 455 PardisoLLT<SpMat, Lower> solver; 456 stat[EIGEN_PARDISO_LLT] = call_directsolver(solver, A, b, refX); 457 printStatItem(stat,EIGEN_PARDISO_LLT, best_time_id, best_time_val); 458 } 459 #endif 460 461 // Internal CG 462 { 463 cout << "\nSolving with CG ... \n"; 464 ConjugateGradient<SpMat, Lower> solver; 465 stat[EIGEN_CG] = call_itersolver(solver, A, b, refX); 466 stat[EIGEN_CG].isIterative = 1; 467 printStatItem(stat,EIGEN_CG, best_time_id, best_time_val); 468 } 469 //CG+IdentityPreconditioner 470 // { 471 // cout << "\nSolving with CG and IdentityPreconditioner ... \n"; 472 // ConjugateGradient<SpMat, Lower, IdentityPreconditioner> solver; 473 // stat[EIGEN_CG_PRECOND] = call_itersolver(solver, A, b, refX); 474 // stat[EIGEN_CG_PRECOND].isIterative = 1; 475 // printStatItem(stat,EIGEN_CG_PRECOND, best_time_id, best_time_val); 476 // } 477 } // End SPD matrices 478 479 return best_time_id; 480 } 481 482 /* Browse all the matrices available in the specified folder 483 * and solve the associated linear system. 484 * The results of each solve are printed in the standard output 485 * and optionally in the provided html file 486 */ 487 template <typename Scalar> 488 void Browse_Matrices(const string folder, bool statFileExists, std::string& statFile, int maxiters, double tol) 489 { 490 MaximumIters = maxiters; // Maximum number of iterations, global variable 491 RelErr = tol; //Relative residual error as stopping criterion for iterative solvers 492 MatrixMarketIterator<Scalar> it(folder); 493 Stats stat[EIGEN_ALL_SOLVERS]; 494 for ( ; it; ++it) 495 { 496 for (int i = 0; i < EIGEN_ALL_SOLVERS; i++) 497 { 498 stat[i].isavail = 0; 499 stat[i].isIterative = 0; 500 } 501 502 int best_time_id; 503 cout<< "\n\n===================================================== \n"; 504 cout<< " ====== SOLVING WITH MATRIX " << it.matname() << " ====\n"; 505 cout<< " =================================================== \n\n"; 506 Matrix<Scalar, Dynamic, 1> refX; 507 if(it.hasrefX()) refX = it.refX(); 508 best_time_id = SelectSolvers<Scalar>(it.matrix(), it.sym(), it.rhs(), refX, &stat[0]); 509 510 if(statFileExists) 511 { 512 string statline; 513 printHtmlStatLine(&stat[0], best_time_id, statline); 514 std::ofstream statbuf(statFile.c_str(), std::ios::app); 515 statbuf << "<TR><TH rowspan=\"4\">" << it.matname() << " <TD rowspan=\"4\"> " 516 << it.matrix().rows() << " <TD rowspan=\"4\"> " << it.matrix().nonZeros()<< " "<< statline ; 517 statbuf.close(); 518 } 519 } 520 } 521 522 bool get_options(int argc, char **args, string option, string* value=0) 523 { 524 int idx = 1, found=false; 525 while (idx<argc && !found){ 526 if (option.compare(args[idx]) == 0){ 527 found = true; 528 if(value) *value = args[idx+1]; 529 } 530 idx+=2; 531 } 532 return found; 533 } 534
