1// Small bench routine for Eigen available in Eigen 2// (C) Desire NUENTSA WAKAM, INRIA 3 4#include <iostream> 5#include <fstream> 6#include <iomanip> 7#include <unsupported/Eigen/SparseExtra> 8#include <Eigen/SparseLU> 9#include <bench/BenchTimer.h> 10#ifdef EIGEN_METIS_SUPPORT 11#include <Eigen/MetisSupport> 12#endif 13 14using namespace std; 15using namespace Eigen; 16 17int main(int argc, char **args) 18{ 19// typedef complex<double> scalar; 20 typedef double scalar; 21 SparseMatrix<scalar, ColMajor> A; 22 typedef SparseMatrix<scalar, ColMajor>::Index Index; 23 typedef Matrix<scalar, Dynamic, Dynamic> DenseMatrix; 24 typedef Matrix<scalar, Dynamic, 1> DenseRhs; 25 Matrix<scalar, Dynamic, 1> b, x, tmp; 26// SparseLU<SparseMatrix<scalar, ColMajor>, AMDOrdering<int> > solver; 27// #ifdef EIGEN_METIS_SUPPORT 28// SparseLU<SparseMatrix<scalar, ColMajor>, MetisOrdering<int> > solver; 29// std::cout<< "ORDERING : METIS\n"; 30// #else 31 SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<int> > solver; 32 std::cout<< "ORDERING : COLAMD\n"; 33// #endif 34 35 ifstream matrix_file; 36 string line; 37 int n; 38 BenchTimer timer; 39 40 // Set parameters 41 /* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */ 42 if (argc < 2) assert(false && "please, give the matrix market file "); 43 loadMarket(A, args[1]); 44 cout << "End charging matrix " << endl; 45 bool iscomplex=false, isvector=false; 46 int sym; 47 getMarketHeader(args[1], sym, iscomplex, isvector); 48// if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; } 49 if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;} 50 if (sym != 0) { // symmetric matrices, only the lower part is stored 51 SparseMatrix<scalar, ColMajor> temp; 52 temp = A; 53 A = temp.selfadjointView<Lower>(); 54 } 55 n = A.cols(); 56 /* Fill the right hand side */ 57 58 if (argc > 2) 59 loadMarketVector(b, args[2]); 60 else 61 { 62 b.resize(n); 63 tmp.resize(n); 64// tmp.setRandom(); 65 for (int i = 0; i < n; i++) tmp(i) = i; 66 b = A * tmp ; 67 } 68 69 /* Compute the factorization */ 70// solver.isSymmetric(true); 71 timer.start(); 72// solver.compute(A); 73 solver.analyzePattern(A); 74 timer.stop(); 75 cout << "Time to analyze " << timer.value() << std::endl; 76 timer.reset(); 77 timer.start(); 78 solver.factorize(A); 79 timer.stop(); 80 cout << "Factorize Time " << timer.value() << std::endl; 81 timer.reset(); 82 timer.start(); 83 x = solver.solve(b); 84 timer.stop(); 85 cout << "solve time " << timer.value() << std::endl; 86 /* Check the accuracy */ 87 Matrix<scalar, Dynamic, 1> tmp2 = b - A*x; 88 scalar tempNorm = tmp2.norm()/b.norm(); 89 cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; 90 cout << "Number of nonzeros in the factor : " << solver.nnzL() + solver.nnzU() << std::endl; 91 92 return 0; 93}