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