1// #define EIGEN_TAUCS_SUPPORT
2// #define EIGEN_CHOLMOD_SUPPORT
3#include <iostream>
4#include <Eigen/Sparse>
5
6// g++ -DSIZE=10000 -DDENSITY=0.001  sparse_cholesky.cpp -I.. -DDENSEMATRI -O3 -g0 -DNDEBUG   -DNBTRIES=1 -I /home/gael/Coding/LinearAlgebra/taucs_full/src/ -I/home/gael/Coding/LinearAlgebra/taucs_full/build/linux/  -L/home/gael/Coding/LinearAlgebra/taucs_full/lib/linux/ -ltaucs /home/gael/Coding/LinearAlgebra/GotoBLAS/libgoto.a -lpthread -I /home/gael/Coding/LinearAlgebra/SuiteSparse/CHOLMOD/Include/ $CHOLLIB -I /home/gael/Coding/LinearAlgebra/SuiteSparse/UFconfig/ /home/gael/Coding/LinearAlgebra/SuiteSparse/CCOLAMD/Lib/libccolamd.a   /home/gael/Coding/LinearAlgebra/SuiteSparse/CHOLMOD/Lib/libcholmod.a -lmetis /home/gael/Coding/LinearAlgebra/SuiteSparse/AMD/Lib/libamd.a  /home/gael/Coding/LinearAlgebra/SuiteSparse/CAMD/Lib/libcamd.a   /home/gael/Coding/LinearAlgebra/SuiteSparse/CCOLAMD/Lib/libccolamd.a  /home/gael/Coding/LinearAlgebra/SuiteSparse/COLAMD/Lib/libcolamd.a -llapack && ./a.out
7
8#define NOGMM
9#define NOMTL
10
11#ifndef SIZE
12#define SIZE 10
13#endif
14
15#ifndef DENSITY
16#define DENSITY 0.01
17#endif
18
19#ifndef REPEAT
20#define REPEAT 1
21#endif
22
23#include "BenchSparseUtil.h"
24
25#ifndef MINDENSITY
26#define MINDENSITY 0.0004
27#endif
28
29#ifndef NBTRIES
30#define NBTRIES 10
31#endif
32
33#define BENCH(X) \
34  timer.reset(); \
35  for (int _j=0; _j<NBTRIES; ++_j) { \
36    timer.start(); \
37    for (int _k=0; _k<REPEAT; ++_k) { \
38        X  \
39  } timer.stop(); }
40
41// typedef SparseMatrix<Scalar,UpperTriangular> EigenSparseTriMatrix;
42typedef SparseMatrix<Scalar,SelfAdjoint|LowerTriangular> EigenSparseSelfAdjointMatrix;
43
44void fillSpdMatrix(float density, int rows, int cols,  EigenSparseSelfAdjointMatrix& dst)
45{
46  dst.startFill(rows*cols*density);
47  for(int j = 0; j < cols; j++)
48  {
49    dst.fill(j,j) = internal::random<Scalar>(10,20);
50    for(int i = j+1; i < rows; i++)
51    {
52      Scalar v = (internal::random<float>(0,1) < density) ? internal::random<Scalar>() : 0;
53      if (v!=0)
54        dst.fill(i,j) = v;
55    }
56
57  }
58  dst.endFill();
59}
60
61#include <Eigen/Cholesky>
62
63template<int Backend>
64void doEigen(const char* name, const EigenSparseSelfAdjointMatrix& sm1, int flags = 0)
65{
66  std::cout << name << "..." << std::flush;
67  BenchTimer timer;
68  timer.start();
69  SparseLLT<EigenSparseSelfAdjointMatrix,Backend> chol(sm1, flags);
70  timer.stop();
71  std::cout << ":\t" << timer.value() << endl;
72
73  std::cout << "  nnz: " << sm1.nonZeros() << " => " << chol.matrixL().nonZeros() << "\n";
74//   std::cout << "sparse\n" << chol.matrixL() << "%\n";
75}
76
77int main(int argc, char *argv[])
78{
79  int rows = SIZE;
80  int cols = SIZE;
81  float density = DENSITY;
82  BenchTimer timer;
83
84  VectorXf b = VectorXf::Random(cols);
85  VectorXf x = VectorXf::Random(cols);
86
87  bool densedone = false;
88
89  //for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
90//   float density = 0.5;
91  {
92    EigenSparseSelfAdjointMatrix sm1(rows, cols);
93    std::cout << "Generate sparse matrix (might take a while)...\n";
94    fillSpdMatrix(density, rows, cols, sm1);
95    std::cout << "DONE\n\n";
96
97    // dense matrices
98    #ifdef DENSEMATRIX
99    if (!densedone)
100    {
101      densedone = true;
102      std::cout << "Eigen Dense\t" << density*100 << "%\n";
103      DenseMatrix m1(rows,cols);
104      eiToDense(sm1, m1);
105      m1 = (m1 + m1.transpose()).eval();
106      m1.diagonal() *= 0.5;
107
108//       BENCH(LLT<DenseMatrix> chol(m1);)
109//       std::cout << "dense:\t" << timer.value() << endl;
110
111      BenchTimer timer;
112      timer.start();
113      LLT<DenseMatrix> chol(m1);
114      timer.stop();
115      std::cout << "dense:\t" << timer.value() << endl;
116      int count = 0;
117      for (int j=0; j<cols; ++j)
118        for (int i=j; i<rows; ++i)
119          if (!internal::isMuchSmallerThan(internal::abs(chol.matrixL()(i,j)), 0.1))
120            count++;
121      std::cout << "dense: " << "nnz = " << count << "\n";
122//       std::cout << "dense:\n" << m1 << "\n\n" << chol.matrixL() << endl;
123    }
124    #endif
125
126    // eigen sparse matrices
127    doEigen<Eigen::DefaultBackend>("Eigen/Sparse", sm1, Eigen::IncompleteFactorization);
128
129    #ifdef EIGEN_CHOLMOD_SUPPORT
130    doEigen<Eigen::Cholmod>("Eigen/Cholmod", sm1, Eigen::IncompleteFactorization);
131    #endif
132
133    #ifdef EIGEN_TAUCS_SUPPORT
134    doEigen<Eigen::Taucs>("Eigen/Taucs", sm1, Eigen::IncompleteFactorization);
135    #endif
136
137    #if 0
138    // TAUCS
139    {
140      taucs_ccs_matrix A = sm1.asTaucsMatrix();
141
142      //BENCH(taucs_ccs_matrix* chol = taucs_ccs_factor_llt(&A, 0, 0);)
143//       BENCH(taucs_supernodal_factor_to_ccs(taucs_ccs_factor_llt_ll(&A));)
144//       std::cout << "taucs:\t" << timer.value() << endl;
145
146      taucs_ccs_matrix* chol = taucs_ccs_factor_llt(&A, 0, 0);
147
148      for (int j=0; j<cols; ++j)
149      {
150        for (int i=chol->colptr[j]; i<chol->colptr[j+1]; ++i)
151          std::cout << chol->values.d[i] << " ";
152      }
153    }
154
155    // CHOLMOD
156    #ifdef EIGEN_CHOLMOD_SUPPORT
157    {
158      cholmod_common c;
159      cholmod_start (&c);
160      cholmod_sparse A;
161      cholmod_factor *L;
162
163      A = sm1.asCholmodMatrix();
164      BenchTimer timer;
165//       timer.reset();
166      timer.start();
167      std::vector<int> perm(cols);
168//       std::vector<int> set(ncols);
169      for (int i=0; i<cols; ++i)
170        perm[i] = i;
171//       c.nmethods = 1;
172//       c.method[0] = 1;
173
174      c.nmethods = 1;
175      c.method [0].ordering = CHOLMOD_NATURAL;
176      c.postorder = 0;
177      c.final_ll = 1;
178
179      L = cholmod_analyze_p(&A, &perm[0], &perm[0], cols, &c);
180      timer.stop();
181      std::cout << "cholmod/analyze:\t" << timer.value() << endl;
182      timer.reset();
183      timer.start();
184      cholmod_factorize(&A, L, &c);
185      timer.stop();
186      std::cout << "cholmod/factorize:\t" << timer.value() << endl;
187
188      cholmod_sparse* cholmat = cholmod_factor_to_sparse(L, &c);
189
190      cholmod_print_factor(L, "Factors", &c);
191
192      cholmod_print_sparse(cholmat, "Chol", &c);
193      cholmod_write_sparse(stdout, cholmat, 0, 0, &c);
194//
195//       cholmod_print_sparse(&A, "A", &c);
196//       cholmod_write_sparse(stdout, &A, 0, 0, &c);
197
198
199//       for (int j=0; j<cols; ++j)
200//       {
201//           for (int i=chol->colptr[j]; i<chol->colptr[j+1]; ++i)
202//             std::cout << chol->values.s[i] << " ";
203//       }
204    }
205    #endif
206
207    #endif
208
209
210
211  }
212
213
214  return 0;
215}
216
217