test_sparseLU.cpp revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba 
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}