1
2 //g++ -O3 -g0 -DNDEBUG sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.005 -DSIZE=10000 && ./a.out
3 //g++ -O3 -g0 -DNDEBUG sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out
4 // -DNOGMM -DNOMTL -DCSPARSE
5 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
6 #ifndef SIZE
7 #define SIZE 650000
8 #endif
9
10 #ifndef DENSITY
11 #define DENSITY 0.01
12 #endif
13
14 #ifndef REPEAT
15 #define REPEAT 1
16 #endif
17
18 #include "BenchSparseUtil.h"
19
20 #ifndef MINDENSITY
21 #define MINDENSITY 0.0004
22 #endif
23
24 #ifndef NBTRIES
25 #define NBTRIES 10
26 #endif
27
28 #define BENCH(X) \
29 timer.reset(); \
30 for (int _j=0; _j<NBTRIES; ++_j) { \
31 timer.start(); \
32 for (int _k=0; _k<REPEAT; ++_k) { \
33 X \
34 } timer.stop(); }
35
36
37 #ifdef CSPARSE
cs_sorted_multiply(const cs * a,const cs * b)38 cs* cs_sorted_multiply(const cs* a, const cs* b)
39 {
40 cs* A = cs_transpose (a, 1) ;
41 cs* B = cs_transpose (b, 1) ;
42 cs* D = cs_multiply (B,A) ; /* D = B'*A' */
43 cs_spfree (A) ;
44 cs_spfree (B) ;
45 cs_dropzeros (D) ; /* drop zeros from D */
46 cs* C = cs_transpose (D, 1) ; /* C = D', so that C is sorted */
47 cs_spfree (D) ;
48 return C;
49 }
50 #endif
51
main(int argc,char * argv[])52 int main(int argc, char *argv[])
53 {
54 int rows = SIZE;
55 int cols = SIZE;
56 float density = DENSITY;
57
58 EigenSparseMatrix sm1(rows,cols);
59 DenseVector v1(cols), v2(cols);
60 v1.setRandom();
61
62 BenchTimer timer;
63 for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
64 {
65 //fillMatrix(density, rows, cols, sm1);
66 fillMatrix2(7, rows, cols, sm1);
67
68 // dense matrices
69 #ifdef DENSEMATRIX
70 {
71 std::cout << "Eigen Dense\t" << density*100 << "%\n";
72 DenseMatrix m1(rows,cols);
73 eiToDense(sm1, m1);
74
75 timer.reset();
76 timer.start();
77 for (int k=0; k<REPEAT; ++k)
78 v2 = m1 * v1;
79 timer.stop();
80 std::cout << " a * v:\t" << timer.best() << " " << double(REPEAT)/timer.best() << " * / sec " << endl;
81
82 timer.reset();
83 timer.start();
84 for (int k=0; k<REPEAT; ++k)
85 v2 = m1.transpose() * v1;
86 timer.stop();
87 std::cout << " a' * v:\t" << timer.best() << endl;
88 }
89 #endif
90
91 // eigen sparse matrices
92 {
93 std::cout << "Eigen sparse\t" << sm1.nonZeros()/float(sm1.rows()*sm1.cols())*100 << "%\n";
94
95 BENCH(asm("#myc"); v2 = sm1 * v1; asm("#myd");)
96 std::cout << " a * v:\t" << timer.best()/REPEAT << " " << double(REPEAT)/timer.best(REAL_TIMER) << " * / sec " << endl;
97
98
99 BENCH( { asm("#mya"); v2 = sm1.transpose() * v1; asm("#myb"); })
100
101 std::cout << " a' * v:\t" << timer.best()/REPEAT << endl;
102 }
103
104 // {
105 // DynamicSparseMatrix<Scalar> m1(sm1);
106 // std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/float(m1.rows()*m1.cols())*100 << "%\n";
107 //
108 // BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1 * v1;)
109 // std::cout << " a * v:\t" << timer.value() << endl;
110 //
111 // BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1.transpose() * v1;)
112 // std::cout << " a' * v:\t" << timer.value() << endl;
113 // }
114
115 // GMM++
116 #ifndef NOGMM
117 {
118 std::cout << "GMM++ sparse\t" << density*100 << "%\n";
119 //GmmDynSparse gmmT3(rows,cols);
120 GmmSparse m1(rows,cols);
121 eiToGmm(sm1, m1);
122
123 std::vector<Scalar> gmmV1(cols), gmmV2(cols);
124 Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1;
125 Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2;
126
127 BENCH( asm("#myx"); gmm::mult(m1, gmmV1, gmmV2); asm("#myy"); )
128 std::cout << " a * v:\t" << timer.value() << endl;
129
130 BENCH( gmm::mult(gmm::transposed(m1), gmmV1, gmmV2); )
131 std::cout << " a' * v:\t" << timer.value() << endl;
132 }
133 #endif
134
135 #ifndef NOUBLAS
136 {
137 std::cout << "ublas sparse\t" << density*100 << "%\n";
138 UBlasSparse m1(rows,cols);
139 eiToUblas(sm1, m1);
140
141 boost::numeric::ublas::vector<Scalar> uv1, uv2;
142 eiToUblasVec(v1,uv1);
143 eiToUblasVec(v2,uv2);
144
145 // std::vector<Scalar> gmmV1(cols), gmmV2(cols);
146 // Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1;
147 // Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2;
148
149 BENCH( uv2 = boost::numeric::ublas::prod(m1, uv1); )
150 std::cout << " a * v:\t" << timer.value() << endl;
151
152 // BENCH( boost::ublas::prod(gmm::transposed(m1), gmmV1, gmmV2); )
153 // std::cout << " a' * v:\t" << timer.value() << endl;
154 }
155 #endif
156
157 // MTL4
158 #ifndef NOMTL
159 {
160 std::cout << "MTL4\t" << density*100 << "%\n";
161 MtlSparse m1(rows,cols);
162 eiToMtl(sm1, m1);
163 mtl::dense_vector<Scalar> mtlV1(cols, 1.0);
164 mtl::dense_vector<Scalar> mtlV2(cols, 1.0);
165
166 timer.reset();
167 timer.start();
168 for (int k=0; k<REPEAT; ++k)
169 mtlV2 = m1 * mtlV1;
170 timer.stop();
171 std::cout << " a * v:\t" << timer.value() << endl;
172
173 timer.reset();
174 timer.start();
175 for (int k=0; k<REPEAT; ++k)
176 mtlV2 = trans(m1) * mtlV1;
177 timer.stop();
178 std::cout << " a' * v:\t" << timer.value() << endl;
179 }
180 #endif
181
182 std::cout << "\n\n";
183 }
184
185 return 0;
186 }
187
188