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