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}