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