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 
14 using namespace std;
15 using namespace Eigen;
16 
17 int 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 }