vlsvrs 1.0.0

VLSV reader with C and Python bindings
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132


// g++ -I.. sparse_lu.cpp -O3 -g0 -I /usr/include/superlu/ -lsuperlu -lgfortran -DSIZE=1000 -DDENSITY=.05 && ./a.out

#define EIGEN_SUPERLU_SUPPORT
#define EIGEN_UMFPACK_SUPPORT
#include <Eigen/Sparse>

#define NOGMM
#define NOMTL

#ifndef SIZE
#define SIZE 10
#endif

#ifndef DENSITY
#define DENSITY 0.01
#endif

#ifndef REPEAT
#define REPEAT 1
#endif

#include "BenchSparseUtil.h"

#ifndef MINDENSITY
#define MINDENSITY 0.0004
#endif

#ifndef NBTRIES
#define NBTRIES 10
#endif

#define BENCH(X) \
  timer.reset(); \
  for (int _j=0; _j<NBTRIES; ++_j) { \
    timer.start(); \
    for (int _k=0; _k<REPEAT; ++_k) { \
        X  \
  } timer.stop(); }

typedef Matrix<Scalar,Dynamic,1> VectorX;

#include <Eigen/LU>

template<int Backend>
void doEigen(const char* name, const EigenSparseMatrix& sm1, const VectorX& b, VectorX& x, int flags = 0)
{
  std::cout << name << "..." << std::flush;
  BenchTimer timer; timer.start();
  SparseLU<EigenSparseMatrix,Backend> lu(sm1, flags);
  timer.stop();
  if (lu.succeeded())
    std::cout << ":\t" << timer.value() << endl;
  else
  {
    std::cout << ":\t FAILED" << endl;
    return;
  }

  bool ok;
  timer.reset(); timer.start();
  ok = lu.solve(b,&x);
  timer.stop();
  if (ok)
    std::cout << "  solve:\t" << timer.value() << endl;
  else
    std::cout << "  solve:\t" << " FAILED" << endl;

  //std::cout << x.transpose() << "\n";
}

int main(int argc, char *argv[])
{
  int rows = SIZE;
  int cols = SIZE;
  float density = DENSITY;
  BenchTimer timer;

  VectorX b = VectorX::Random(cols);
  VectorX x = VectorX::Random(cols);

  bool densedone = false;

  //for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
//   float density = 0.5;
  {
    EigenSparseMatrix sm1(rows, cols);
    fillMatrix(density, rows, cols, sm1);

    // dense matrices
    #ifdef DENSEMATRIX
    if (!densedone)
    {
      densedone = true;
      std::cout << "Eigen Dense\t" << density*100 << "%\n";
      DenseMatrix m1(rows,cols);
      eiToDense(sm1, m1);

      BenchTimer timer;
      timer.start();
      FullPivLU<DenseMatrix> lu(m1);
      timer.stop();
      std::cout << "Eigen/dense:\t" << timer.value() << endl;

      timer.reset();
      timer.start();
      lu.solve(b,&x);
      timer.stop();
      std::cout << "  solve:\t" << timer.value() << endl;
//       std::cout << b.transpose() << "\n";
//       std::cout << x.transpose() << "\n";
    }
    #endif

    #ifdef EIGEN_UMFPACK_SUPPORT
    x.setZero();
    doEigen<Eigen::UmfPack>("Eigen/UmfPack (auto)", sm1, b, x, 0);
    #endif

    #ifdef EIGEN_SUPERLU_SUPPORT
    x.setZero();
    doEigen<Eigen::SuperLU>("Eigen/SuperLU (nat)", sm1, b, x, Eigen::NaturalOrdering);
//     doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD AT+A)", sm1, b, x, Eigen::MinimumDegree_AT_PLUS_A);
//     doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD ATA)", sm1, b, x, Eigen::MinimumDegree_ATA);
    doEigen<Eigen::SuperLU>("Eigen/SuperLU (COLAMD)", sm1, b, x, Eigen::ColApproxMinimumDegree);
    #endif

  }

  return 0;
}