• 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
5 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
6 
7 #ifndef SIZE
8 #define SIZE 10000
9 #endif
10 
11 #ifndef DENSITY
12 #define DENSITY 0.01
13 #endif
14 
15 #ifndef REPEAT
16 #define REPEAT 1
17 #endif
18 
19 #include "BenchSparseUtil.h"
20 
21 #ifndef MINDENSITY
22 #define MINDENSITY 0.0004
23 #endif
24 
25 #ifndef NBTRIES
26 #define NBTRIES 10
27 #endif
28 
29 #define BENCH(X) \
30   timer.reset(); \
31   for (int _j=0; _j<NBTRIES; ++_j) { \
32     timer.start(); \
33     for (int _k=0; _k<REPEAT; ++_k) { \
34         X  \
35   } timer.stop(); }
36 
37 typedef SparseMatrix<Scalar,UpperTriangular> EigenSparseTriMatrix;
38 typedef SparseMatrix<Scalar,RowMajorBit|UpperTriangular> EigenSparseTriMatrixRow;
39 
fillMatrix(float density,int rows,int cols,EigenSparseTriMatrix & dst)40 void fillMatrix(float density, int rows, int cols,  EigenSparseTriMatrix& dst)
41 {
42   dst.startFill(rows*cols*density);
43   for(int j = 0; j < cols; j++)
44   {
45     for(int i = 0; i < j; i++)
46     {
47       Scalar v = (internal::random<float>(0,1) < density) ? internal::random<Scalar>() : 0;
48       if (v!=0)
49         dst.fill(i,j) = v;
50     }
51     dst.fill(j,j) = internal::random<Scalar>();
52   }
53   dst.endFill();
54 }
55 
main(int argc,char * argv[])56 int main(int argc, char *argv[])
57 {
58   int rows = SIZE;
59   int cols = SIZE;
60   float density = DENSITY;
61   BenchTimer timer;
62   #if 1
63   EigenSparseTriMatrix sm1(rows,cols);
64   typedef Matrix<Scalar,Dynamic,1> DenseVector;
65   DenseVector b = DenseVector::Random(cols);
66   DenseVector x = DenseVector::Random(cols);
67 
68   bool densedone = false;
69 
70   for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
71   {
72     EigenSparseTriMatrix sm1(rows, cols);
73     fillMatrix(density, rows, cols, sm1);
74 
75     // dense matrices
76     #ifdef DENSEMATRIX
77     if (!densedone)
78     {
79       densedone = true;
80       std::cout << "Eigen Dense\t" << density*100 << "%\n";
81       DenseMatrix m1(rows,cols);
82       Matrix<Scalar,Dynamic,Dynamic,Dynamic,Dynamic,RowMajorBit> m2(rows,cols);
83       eiToDense(sm1, m1);
84       m2 = m1;
85 
86       BENCH(x = m1.marked<UpperTriangular>().solveTriangular(b);)
87       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
88 //       std::cerr << x.transpose() << "\n";
89 
90       BENCH(x = m2.marked<UpperTriangular>().solveTriangular(b);)
91       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
92 //       std::cerr << x.transpose() << "\n";
93     }
94     #endif
95 
96     // eigen sparse matrices
97     {
98       std::cout << "Eigen sparse\t" << density*100 << "%\n";
99       EigenSparseTriMatrixRow sm2 = sm1;
100 
101       BENCH(x = sm1.solveTriangular(b);)
102       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
103 //       std::cerr << x.transpose() << "\n";
104 
105       BENCH(x = sm2.solveTriangular(b);)
106       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
107 //       std::cerr << x.transpose() << "\n";
108 
109 //       x = b;
110 //       BENCH(sm1.inverseProductInPlace(x);)
111 //       std::cout << "   colmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl;
112 //       std::cerr << x.transpose() << "\n";
113 //
114 //       x = b;
115 //       BENCH(sm2.inverseProductInPlace(x);)
116 //       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl;
117 //       std::cerr << x.transpose() << "\n";
118     }
119 
120 
121 
122     // CSparse
123     #ifdef CSPARSE
124     {
125       std::cout << "CSparse \t" << density*100 << "%\n";
126       cs *m1;
127       eiToCSparse(sm1, m1);
128 
129       BENCH(x = b; if (!cs_lsolve (m1, x.data())){std::cerr << "cs_lsolve failed\n"; break;}; )
130       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
131     }
132     #endif
133 
134     // GMM++
135     #ifndef NOGMM
136     {
137       std::cout << "GMM++ sparse\t" << density*100 << "%\n";
138       GmmSparse m1(rows,cols);
139       gmm::csr_matrix<Scalar> m2;
140       eiToGmm(sm1, m1);
141       gmm::copy(m1,m2);
142       std::vector<Scalar> gmmX(cols), gmmB(cols);
143       Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols) = x;
144       Map<Matrix<Scalar,Dynamic,1> >(&gmmB[0], cols) = b;
145 
146       gmmX = gmmB;
147       BENCH(gmm::upper_tri_solve(m1, gmmX, false);)
148       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
149 //       std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n";
150 
151       gmmX = gmmB;
152       BENCH(gmm::upper_tri_solve(m2, gmmX, false);)
153       timer.stop();
154       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
155 //       std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n";
156     }
157     #endif
158 
159     // MTL4
160     #ifndef NOMTL
161     {
162       std::cout << "MTL4\t" << density*100 << "%\n";
163       MtlSparse m1(rows,cols);
164       MtlSparseRowMajor m2(rows,cols);
165       eiToMtl(sm1, m1);
166       m2 = m1;
167       mtl::dense_vector<Scalar> x(rows, 1.0);
168       mtl::dense_vector<Scalar> b(rows, 1.0);
169 
170       BENCH(x = mtl::upper_trisolve(m1,b);)
171       std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
172 //       std::cerr << x << "\n";
173 
174       BENCH(x = mtl::upper_trisolve(m2,b);)
175       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
176 //       std::cerr << x << "\n";
177     }
178     #endif
179 
180 
181     std::cout << "\n\n";
182   }
183   #endif
184 
185   #if 0
186     // bench small matrices (in-place versus return bye value)
187     {
188       timer.reset();
189       for (int _j=0; _j<10; ++_j) {
190         Matrix4f m = Matrix4f::Random();
191         Vector4f b = Vector4f::Random();
192         Vector4f x = Vector4f::Random();
193         timer.start();
194         for (int _k=0; _k<1000000; ++_k) {
195           b = m.inverseProduct(b);
196         }
197         timer.stop();
198       }
199       std::cout << "4x4 :\t" << timer.value() << endl;
200     }
201 
202     {
203       timer.reset();
204       for (int _j=0; _j<10; ++_j) {
205         Matrix4f m = Matrix4f::Random();
206         Vector4f b = Vector4f::Random();
207         Vector4f x = Vector4f::Random();
208         timer.start();
209         for (int _k=0; _k<1000000; ++_k) {
210           m.inverseProductInPlace(x);
211         }
212         timer.stop();
213       }
214       std::cout << "4x4 IP :\t" << timer.value() << endl;
215     }
216   #endif
217 
218   return 0;
219 }
220 
221