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 217
