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 75using namespace Eigen; 76using namespace std; 77 78 79// Global variables for input parameters 80int MaximumIters; // Maximum number of iterations 81double RelErr; // Relative error of the computed solution 82double best_time_val; // Current best time overall solvers 83int best_time_id; // id of the best solver for the current system 84 85template<typename T> inline typename NumTraits<T>::Real test_precision() { return NumTraits<T>::dummy_precision(); } 86template<> inline float test_precision<float>() { return 1e-3f; } 87template<> inline double test_precision<double>() { return 1e-6; } 88template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); } 89template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); } 90 91void 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 213template<typename Solver, typename Scalar> 214void 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 278template<typename Solver, typename Scalar> 279void 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 288template<typename Solver, typename Scalar> 289void 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 304template <typename Scalar> 305void 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 */ 489template <typename Scalar> 490void 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 543bool 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