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
fillSpdMatrix(float density,int rows,int cols,EigenSparseSelfAdjointMatrix & dst)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>
doEigen(const char * name,const EigenSparseSelfAdjointMatrix & sm1,int flags=0)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
main(int argc,char * argv[])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
217