test_sparseLU.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 <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;
24  typedef Matrix<scalar, Dynamic, 1> DenseRhs;
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
32  std::cout<< "ORDERING : COLAMD\n";
33 // #endif
34 
35  ifstream matrix_file;
36  string line;
37  int n;
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
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 }
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 nnzU() const
Definition: SparseLU.h:456
Definition: pytypes.h:2012
int n
const Solve< SparseLU< _MatrixType, _OrderingType >, 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
Definition: BFloat16.h:88
BiCGSTAB< SparseMatrix< double > > solver
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
void factorize(const MatrixType &matrix)
Definition: SparseLU.h:595
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
bool getMarketHeader(const std::string &filename, int &sym, bool &iscomplex, bool &isvector)
Definition: MarketIO.h:109
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Index nnzL() const
Definition: SparseLU.h:455
ConstSelfAdjointViewReturnType< UpLo >::Type selfadjointView() const
int main(int argc, char **args)
bool loadMarket(SparseMatrixType &mat, const std::string &filename)
Definition: MarketIO.h:134
void analyzePattern(const MatrixType &matrix)
Definition: SparseLU.h:510
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
mxArray * scalar(mxClassID classid)
Definition: matlab.h:82
static BenchTimer timer


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:37:46