sp_solver.cpp
Go to the documentation of this file.
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>
17 using namespace std;
18 using namespace Eigen;
19 
20 int main(int argc, char **args)
21 {
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;
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
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 }
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Index cols() const
Definition: SparseMatrix.h:140
bool loadMarketVector(VectorType &vec, const std::string &filename)
Definition: MarketIO.h:200
Scalar * b
Definition: benchVecAdd.cpp:17
A versatible sparse matrix representation.
Definition: SparseMatrix.h:96
Index rows() const
Definition: SparseMatrix.h:138
Modified Incomplete Cholesky with dual threshold.
Definition: pytypes.h:2012
int n
const Solve< ConjugateGradient< _MatrixType, _UpLo, _Preconditioner >, Rhs > solve(const MatrixBase< Rhs > &b) const
double value(int TIMER=CPU_TIMER) const
Definition: BenchTimer.h:104
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: BFloat16.h:88
BiCGSTAB< SparseMatrix< double > > solver
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
bool getMarketHeader(const std::string &filename, int &sym, bool &iscomplex, bool &isvector)
Definition: MarketIO.h:109
int main(int argc, char **args)
Definition: sp_solver.cpp:20
ConjugateGradient< _MatrixType, _UpLo, _Preconditioner > & compute(const EigenBase< MatrixDerived > &A)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
ConstSelfAdjointViewReturnType< UpLo >::Type selfadjointView() const
bool loadMarket(SparseMatrixType &mat, const std::string &filename)
Definition: MarketIO.h:134
The matrix class, also used for vectors and row-vectors.
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
static BenchTimer timer


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:53