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
7#include <typeinfo>
8
9#ifndef SIZE
10#define SIZE 1000000
11#endif
12
13#ifndef NNZPERCOL
14#define NNZPERCOL 6
15#endif
16
17#ifndef REPEAT
18#define REPEAT 1
19#endif
20
21#include <algorithm>
22#include "BenchTimer.h"
23#include "BenchUtil.h"
24#include "BenchSparseUtil.h"
25
26#ifndef NBTRIES
27#define NBTRIES 1
28#endif
29
30#define BENCH(X) \
31  timer.reset(); \
32  for (int _j=0; _j<NBTRIES; ++_j) { \
33    timer.start(); \
34    for (int _k=0; _k<REPEAT; ++_k) { \
35        X  \
36  } timer.stop(); }
37
38// #ifdef MKL
39//
40// #include "mkl_types.h"
41// #include "mkl_spblas.h"
42//
43// template<typename Lhs,typename Rhs,typename Res>
44// void mkl_multiply(const Lhs& lhs, const Rhs& rhs, Res& res)
45// {
46//   char n = 'N';
47//   float alpha = 1;
48//   char matdescra[6];
49//   matdescra[0] = 'G';
50//   matdescra[1] = 0;
51//   matdescra[2] = 0;
52//   matdescra[3] = 'C';
53//   mkl_scscmm(&n, lhs.rows(), rhs.cols(), lhs.cols(), &alpha, matdescra,
54//              lhs._valuePtr(), lhs._innerIndexPtr(), lhs.outerIndexPtr(),
55//              pntre, b, &ldb, &beta, c, &ldc);
56// //   mkl_somatcopy('C', 'T', lhs.rows(), lhs.cols(), 1,
57// //                 lhs._valuePtr(), lhs.rows(), DST, dst_stride);
58// }
59//
60// #endif
61
62
63#ifdef CSPARSE
64cs* cs_sorted_multiply(const cs* a, const cs* b)
65{
66//   return cs_multiply(a,b);
67
68  cs* A = cs_transpose(a, 1);
69  cs* B = cs_transpose(b, 1);
70  cs* D = cs_multiply(B,A);   /* D = B'*A' */
71  cs_spfree (A) ;
72  cs_spfree (B) ;
73  cs_dropzeros (D) ;      /* drop zeros from D */
74  cs* C = cs_transpose (D, 1) ;   /* C = D', so that C is sorted */
75  cs_spfree (D) ;
76  return C;
77
78//   cs* A = cs_transpose(a, 1);
79//   cs* C = cs_transpose(A, 1);
80//   return C;
81}
82
83cs* cs_sorted_multiply2(const cs* a, const cs* b)
84{
85  cs* D = cs_multiply(a,b);
86  cs* E = cs_transpose(D,1);
87  cs_spfree(D);
88  cs* C = cs_transpose(E,1);
89  cs_spfree(E);
90  return C;
91}
92#endif
93
94void bench_sort();
95
96int main(int argc, char *argv[])
97{
98//   bench_sort();
99
100  int rows = SIZE;
101  int cols = SIZE;
102  float density = DENSITY;
103
104  EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols);
105
106  BenchTimer timer;
107  for (int nnzPerCol = NNZPERCOL; nnzPerCol>1; nnzPerCol/=1.1)
108  {
109    sm1.setZero();
110    sm2.setZero();
111    fillMatrix2(nnzPerCol, rows, cols, sm1);
112    fillMatrix2(nnzPerCol, rows, cols, sm2);
113//     std::cerr << "filling OK\n";
114
115    // dense matrices
116    #ifdef DENSEMATRIX
117    {
118      std::cout << "Eigen Dense\t" << nnzPerCol << "%\n";
119      DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
120      eiToDense(sm1, m1);
121      eiToDense(sm2, m2);
122
123      timer.reset();
124      timer.start();
125      for (int k=0; k<REPEAT; ++k)
126        m3 = m1 * m2;
127      timer.stop();
128      std::cout << "   a * b:\t" << timer.value() << endl;
129
130      timer.reset();
131      timer.start();
132      for (int k=0; k<REPEAT; ++k)
133        m3 = m1.transpose() * m2;
134      timer.stop();
135      std::cout << "   a' * b:\t" << timer.value() << endl;
136
137      timer.reset();
138      timer.start();
139      for (int k=0; k<REPEAT; ++k)
140        m3 = m1.transpose() * m2.transpose();
141      timer.stop();
142      std::cout << "   a' * b':\t" << timer.value() << endl;
143
144      timer.reset();
145      timer.start();
146      for (int k=0; k<REPEAT; ++k)
147        m3 = m1 * m2.transpose();
148      timer.stop();
149      std::cout << "   a * b':\t" << timer.value() << endl;
150    }
151    #endif
152
153    // eigen sparse matrices
154    {
155      std::cout << "Eigen sparse\t" << sm1.nonZeros()/(float(sm1.rows())*float(sm1.cols()))*100 << "% * "
156                << sm2.nonZeros()/(float(sm2.rows())*float(sm2.cols()))*100 << "%\n";
157
158      BENCH(sm3 = sm1 * sm2; )
159      std::cout << "   a * b:\t" << timer.value() << endl;
160
161//       BENCH(sm3 = sm1.transpose() * sm2; )
162//       std::cout << "   a' * b:\t" << timer.value() << endl;
163// //
164//       BENCH(sm3 = sm1.transpose() * sm2.transpose(); )
165//       std::cout << "   a' * b':\t" << timer.value() << endl;
166// //
167//       BENCH(sm3 = sm1 * sm2.transpose(); )
168//       std::cout << "   a * b' :\t" << timer.value() << endl;
169
170
171//       std::cout << "\n";
172//
173//       BENCH( sm3._experimentalNewProduct(sm1, sm2); )
174//       std::cout << "   a * b:\t" << timer.value() << endl;
175//
176//       BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2); )
177//       std::cout << "   a' * b:\t" << timer.value() << endl;
178// //
179//       BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2.transpose()); )
180//       std::cout << "   a' * b':\t" << timer.value() << endl;
181// //
182//       BENCH(sm3._experimentalNewProduct(sm1, sm2.transpose());)
183//       std::cout << "   a * b' :\t" << timer.value() << endl;
184    }
185
186    // eigen dyn-sparse matrices
187    /*{
188      DynamicSparseMatrix<Scalar> m1(sm1), m2(sm2), m3(sm3);
189      std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/(float(m1.rows())*float(m1.cols()))*100 << "% * "
190                << m2.nonZeros()/(float(m2.rows())*float(m2.cols()))*100 << "%\n";
191
192//       timer.reset();
193//       timer.start();
194      BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1 * m2;)
195//       timer.stop();
196      std::cout << "   a * b:\t" << timer.value() << endl;
197//       std::cout << sm3 << "\n";
198
199      timer.reset();
200      timer.start();
201//       std::cerr << "transpose...\n";
202//       EigenSparseMatrix sm4 = sm1.transpose();
203//       std::cout << sm4.nonZeros() << " == " << sm1.nonZeros() << "\n";
204//       exit(1);
205//       std::cerr << "transpose OK\n";
206//       std::cout << sm1 << "\n\n" << sm1.transpose() << "\n\n" << sm4.transpose() << "\n\n";
207      BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2;)
208//       timer.stop();
209      std::cout << "   a' * b:\t" << timer.value() << endl;
210
211//       timer.reset();
212//       timer.start();
213      BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2.transpose(); )
214//       timer.stop();
215      std::cout << "   a' * b':\t" << timer.value() << endl;
216
217//       timer.reset();
218//       timer.start();
219      BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1 * m2.transpose(); )
220//       timer.stop();
221      std::cout << "   a * b' :\t" << timer.value() << endl;
222    }*/
223
224    // CSparse
225    #ifdef CSPARSE
226    {
227      std::cout << "CSparse \t" << nnzPerCol << "%\n";
228      cs *m1, *m2, *m3;
229      eiToCSparse(sm1, m1);
230      eiToCSparse(sm2, m2);
231
232      BENCH(
233      {
234        m3 = cs_sorted_multiply(m1, m2);
235        if (!m3)
236        {
237          std::cerr << "cs_multiply failed\n";
238        }
239//         cs_print(m3, 0);
240        cs_spfree(m3);
241      }
242      );
243//       timer.stop();
244      std::cout << "   a * b:\t" << timer.value() << endl;
245
246//       BENCH( { m3 = cs_sorted_multiply2(m1, m2); cs_spfree(m3); } );
247//       std::cout << "   a * b:\t" << timer.value() << endl;
248    }
249    #endif
250
251    #ifndef NOUBLAS
252    {
253      std::cout << "ublas\t" << nnzPerCol << "%\n";
254      UBlasSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
255      eiToUblas(sm1, m1);
256      eiToUblas(sm2, m2);
257
258      BENCH(boost::numeric::ublas::prod(m1, m2, m3););
259      std::cout << "   a * b:\t" << timer.value() << endl;
260    }
261    #endif
262
263    // GMM++
264    #ifndef NOGMM
265    {
266      std::cout << "GMM++ sparse\t" << nnzPerCol << "%\n";
267      GmmDynSparse  gmmT3(rows,cols);
268      GmmSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
269      eiToGmm(sm1, m1);
270      eiToGmm(sm2, m2);
271
272      BENCH(gmm::mult(m1, m2, gmmT3););
273      std::cout << "   a * b:\t" << timer.value() << endl;
274
275//       BENCH(gmm::mult(gmm::transposed(m1), m2, gmmT3););
276//       std::cout << "   a' * b:\t" << timer.value() << endl;
277//
278//       if (rows<500)
279//       {
280//         BENCH(gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3););
281//         std::cout << "   a' * b':\t" << timer.value() << endl;
282//
283//         BENCH(gmm::mult(m1, gmm::transposed(m2), gmmT3););
284//         std::cout << "   a * b':\t" << timer.value() << endl;
285//       }
286//       else
287//       {
288//         std::cout << "   a' * b':\t" << "forever" << endl;
289//         std::cout << "   a * b':\t" << "forever" << endl;
290//       }
291    }
292    #endif
293
294    // MTL4
295    #ifndef NOMTL
296    {
297      std::cout << "MTL4\t" << nnzPerCol << "%\n";
298      MtlSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
299      eiToMtl(sm1, m1);
300      eiToMtl(sm2, m2);
301
302      BENCH(m3 = m1 * m2;);
303      std::cout << "   a * b:\t" << timer.value() << endl;
304
305//       BENCH(m3 = trans(m1) * m2;);
306//       std::cout << "   a' * b:\t" << timer.value() << endl;
307//
308//       BENCH(m3 = trans(m1) * trans(m2););
309//       std::cout << "  a' * b':\t" << timer.value() << endl;
310//
311//       BENCH(m3 = m1 * trans(m2););
312//       std::cout << "   a * b' :\t" << timer.value() << endl;
313    }
314    #endif
315
316    std::cout << "\n\n";
317  }
318
319  return 0;
320}
321
322
323
324