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 */
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 }
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::SparseMatrix
A versatible sparse matrix representation.
Definition: SparseMatrix.h:96
timer
static BenchTimer timer
Definition: benchmark-blocking-sizes.cpp:31
b
Scalar * b
Definition: benchVecAdd.cpp:17
x
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
Definition: gnuplot_common_settings.hh:12
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
solver
BiCGSTAB< SparseMatrix< double > > solver
Definition: BiCGSTAB_simple.cpp:5
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
Eigen::loadMarketVector
bool loadMarketVector(VectorType &vec, const std::string &filename)
Definition: MarketIO.h:200
Eigen::BenchTimer::value
double value(int TIMER=CPU_TIMER) const
Definition: BenchTimer.h:104
n
int n
Definition: BiCGSTAB_simple.cpp:1
Eigen::BenchTimer::start
void start()
Definition: BenchTimer.h:81
scalar
mxArray * scalar(mxClassID classid)
Definition: matlab.h:82
Eigen::BenchTimer
Definition: BenchTimer.h:59
Eigen::Lower
@ Lower
Definition: Constants.h:209
Eigen::BenchTimer::reset
void reset()
Definition: BenchTimer.h:75
Eigen::COLAMDOrdering
Definition: 3rdparty/Eigen/Eigen/src/OrderingMethods/Ordering.h:114
std
Definition: BFloat16.h:88
args
Definition: pytypes.h:2210
BenchTimer.h
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
Eigen::getMarketHeader
bool getMarketHeader(const std::string &filename, int &sym, bool &iscomplex, bool &isvector)
Definition: MarketIO.h:109
main
int main(int argc, char **args)
Definition: test_sparseLU.cpp:17
Eigen::BenchTimer::stop
void stop()
Definition: BenchTimer.h:86
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Eigen::SparseLU
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
Eigen::loadMarket
bool loadMarket(SparseMatrixType &mat, const std::string &filename)
Definition: MarketIO.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:07:08