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 <Eigen/Jacobi> 8#include <Eigen/Householder> 9#include <Eigen/IterativeLinearSolvers> 10#include <Eigen/LU> 11#include <unsupported/Eigen/SparseExtra> 12//#include <Eigen/SparseLU> 13#include <Eigen/SuperLUSupport> 14// #include <unsupported/Eigen/src/IterativeSolvers/Scaling.h> 15#include <bench/BenchTimer.h> 16#include <unsupported/Eigen/IterativeSolvers> 17using namespace std; 18using namespace Eigen; 19 20int main(int argc, char **args) 21{ 22 SparseMatrix<double, ColMajor> A; 23 typedef SparseMatrix<double, ColMajor>::Index Index; 24 typedef Matrix<double, Dynamic, Dynamic> DenseMatrix; 25 typedef Matrix<double, Dynamic, 1> DenseRhs; 26 VectorXd b, x, tmp; 27 BenchTimer timer,totaltime; 28 //SparseLU<SparseMatrix<double, ColMajor> > solver; 29// SuperLU<SparseMatrix<double, ColMajor> > solver; 30 ConjugateGradient<SparseMatrix<double, ColMajor>, Lower,IncompleteCholesky<double,Lower> > solver; 31 ifstream matrix_file; 32 string line; 33 int n; 34 // Set parameters 35// solver.iparm(IPARM_THREAD_NBR) = 4; 36 /* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */ 37 if (argc < 2) assert(false && "please, give the matrix market file "); 38 39 timer.start(); 40 totaltime.start(); 41 loadMarket(A, args[1]); 42 cout << "End charging matrix " << endl; 43 bool iscomplex=false, isvector=false; 44 int sym; 45 getMarketHeader(args[1], sym, iscomplex, isvector); 46 if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; } 47 if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;} 48 if (sym != 0) { // symmetric matrices, only the lower part is stored 49 SparseMatrix<double, ColMajor> temp; 50 temp = A; 51 A = temp.selfadjointView<Lower>(); 52 } 53 timer.stop(); 54 55 n = A.cols(); 56 // ====== TESTS FOR SPARSE TUTORIAL ====== 57// cout<< "OuterSize " << A.outerSize() << " inner " << A.innerSize() << endl; 58// SparseMatrix<double, RowMajor> mat1(A); 59// SparseMatrix<double, RowMajor> mat2; 60// cout << " norm of A " << mat1.norm() << endl; ; 61// PermutationMatrix<Dynamic, Dynamic, int> perm(n); 62// perm.resize(n,1); 63// perm.indices().setLinSpaced(n, 0, n-1); 64// mat2 = perm * mat1; 65// mat.subrows(); 66// mat2.resize(n,n); 67// mat2.reserve(10); 68// mat2.setConstant(); 69// std::cout<< "NORM " << mat1.squaredNorm()<< endl; 70 71 cout<< "Time to load the matrix " << timer.value() <<endl; 72 /* Fill the right hand side */ 73 74// solver.set_restart(374); 75 if (argc > 2) 76 loadMarketVector(b, args[2]); 77 else 78 { 79 b.resize(n); 80 tmp.resize(n); 81// tmp.setRandom(); 82 for (int i = 0; i < n; i++) tmp(i) = i; 83 b = A * tmp ; 84 } 85// Scaling<SparseMatrix<double> > scal; 86// scal.computeRef(A); 87// b = scal.LeftScaling().cwiseProduct(b); 88 89 /* Compute the factorization */ 90 cout<< "Starting the factorization "<< endl; 91 timer.reset(); 92 timer.start(); 93 cout<< "Size of Input Matrix "<< b.size()<<"\n\n"; 94 cout<< "Rows and columns "<< A.rows() <<" " <<A.cols() <<"\n"; 95 solver.compute(A); 96// solver.analyzePattern(A); 97// solver.factorize(A); 98 if (solver.info() != Success) { 99 std::cout<< "The solver failed \n"; 100 return -1; 101 } 102 timer.stop(); 103 float time_comp = timer.value(); 104 cout <<" Compute Time " << time_comp<< endl; 105 106 timer.reset(); 107 timer.start(); 108 x = solver.solve(b); 109// x = scal.RightScaling().cwiseProduct(x); 110 timer.stop(); 111 float time_solve = timer.value(); 112 cout<< " Time to solve " << time_solve << endl; 113 114 /* Check the accuracy */ 115 VectorXd tmp2 = b - A*x; 116 double tempNorm = tmp2.norm()/b.norm(); 117 cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; 118// cout << "Iterations : " << solver.iterations() << "\n"; 119 120 totaltime.stop(); 121 cout << "Total time " << totaltime.value() << "\n"; 122// std::cout<<x.transpose()<<"\n"; 123 124 return 0; 125}