sparse_cholesky.cpp
Go to the documentation of this file.
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;
43 
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;
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;
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  {
93  std::cout << "Generate sparse matrix (might take a while)...\n";
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";
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 
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)
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();
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 
doEigen
void doEigen(const char *name, const EigenSparseSelfAdjointMatrix &sm1, int flags=0)
Definition: sparse_cholesky.cpp:64
perm
idx_t idx_t idx_t idx_t idx_t * perm
Definition: include/metis.h:223
gtsam.examples.SFMExample_bal.stdout
stdout
Definition: SFMExample_bal.py:24
name
Annotation for function names.
Definition: attr.h:51
Eigen::SparseMatrix
A versatible sparse matrix representation.
Definition: SparseMatrix.h:96
eiToDense
void eiToDense(const EigenSparseMatrix &src, DenseMatrix &dst)
Definition: BenchSparseUtil.h:62
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
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
m1
Matrix3d m1
Definition: IOFormat.cpp:2
DENSITY
#define DENSITY
Definition: sparse_cholesky.cpp:16
Eigen::SparseMatrix::nonZeros
Index nonZeros() const
Definition: SparseCompressedBase.h:56
Eigen::internal::isMuchSmallerThan
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
Definition: Eigen/src/Core/MathFunctions.h:1940
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
Eigen::BenchTimer::value
double value(int TIMER=CPU_TIMER) const
Definition: BenchTimer.h:104
fillSpdMatrix
void fillSpdMatrix(float density, int rows, int cols, EigenSparseSelfAdjointMatrix &dst)
Definition: sparse_cholesky.cpp:44
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::BenchTimer::start
void start()
Definition: BenchTimer.h:81
Eigen::LLT::matrixL
Traits::MatrixL matrixL() const
Definition: LLT.h:135
Eigen::BenchTimer
Definition: BenchTimer.h:59
density
Definition: testGaussianConditional.cpp:127
L
MatrixXd L
Definition: LLT_example.cpp:6
Eigen::BenchTimer::reset
void reset()
Definition: BenchTimer.h:75
main
int main(int argc, char *argv[])
Definition: sparse_cholesky.cpp:77
Eigen::LLT
Standard Cholesky decomposition (LL^T) of a matrix and associated features.
Definition: LLT.h:66
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
Eigen::BenchTimer::stop
void stop()
Definition: BenchTimer.h:86
abs
#define abs(x)
Definition: datatypes.h:17
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
SIZE
#define SIZE
Definition: sparse_cholesky.cpp:12
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
EigenSparseSelfAdjointMatrix
SparseMatrix< Scalar, SelfAdjoint|LowerTriangular > EigenSparseSelfAdjointMatrix
Definition: sparse_cholesky.cpp:42
BenchSparseUtil.h
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:33