• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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