ground/gcs/src/libs/eigen/bench/sparse_cholesky.cpp

   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;
  42 typedef SparseMatrix<Scalar,SelfAdjoint|LowerTriangular> EigenSparseSelfAdjointMatrix;
  43
  44 void 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
  63 template<int Backend>
  64 void 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
  77 int 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