1/* ----------------------------------------------------------------------------
2
3 * GTSAM Copyright 2010, Georgia Tech Research Corporation,
4 * Atlanta, Georgia 30332-0415
5 * All Rights Reserved
6 * Authors: Frank Dellaert, et al. (see THANKS for the full author list)
7
8 * See LICENSE for the license information
9
10 * -------------------------------------------------------------------------- */
11
12/**
13 * @file timeublas.cpp
14 * @brief Tests to help determine which way of accomplishing something with Eigen is faster
15 * @author Richard Roberts
16 * @date Sep 18, 2010
17 */
18
19#include <gtsam/base/timing.h>
20#include <gtsam/base/Matrix.h>
21
22#include <iostream>
23#include <random>
24#include <vector>
25#include <utility>
26
27using namespace std;
28//namespace ublas = boost::numeric::ublas;
29//using namespace Eigen;
30
31static std::mt19937 rng;
32static std::uniform_real_distribution<> uniform(-1.0, 0.0);
33//typedef ublas::matrix<double> matrix;
34//typedef ublas::matrix_range<matrix> matrix_range;
35//typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix;
36//typedef Eigen::Block<matrix> matrix_block;
37
38//using ublas::range;
39//using ublas::triangular_matrix;
40
41int main(int argc, char* argv[]) {
42
43 if(true) {
44 cout << "\nTiming matrix_block:" << endl;
45
46 // We use volatile here to make these appear to the optimizing compiler as
47 // if their values are only known at run-time.
48 volatile size_t m=500;
49 volatile size_t n=300;
50 volatile size_t nReps = 1000;
51 assert(m > n);
52 std::uniform_int_distribution<size_t> uniform_i(0,m-1);
53 std::uniform_int_distribution<size_t> uniform_j(0,n-1);
54 gtsam::Matrix mat((int)m,(int)n);
55 gtsam::SubMatrix full = mat.block(startRow: 0, startCol: 0, blockRows: m, blockCols: n);
56 gtsam::SubMatrix top = mat.block(startRow: 0, startCol: 0, blockRows: n, blockCols: n);
57 gtsam::SubMatrix block = mat.block(startRow: m/4, startCol: n/4, blockRows: m-m/2, blockCols: n-n/2);
58
59 cout << " Basic: " << (int)m << "x" << (int)n << endl;
60 cout << " Full: mat(" << 0 << ":" << (int)m << ", " << 0 << ":" << (int)n << ")" << endl;
61 cout << " Top: mat(" << 0 << ":" << (int)n << ", " << 0 << ":" << (int)n << ")" << endl;
62 cout << " Block: mat(" << size_t(m/4) << ":" << size_t(m-m/4) << ", " << size_t(n/4) << ":" << size_t(n-n/4) << ")" << endl;
63 cout << endl;
64
65 {
66 double basicTime, fullTime, topTime, blockTime;
67
68 cout << "Row-major matrix, row-major assignment:" << endl;
69
70 // Do a few initial assignments to let any cache effects stabilize
71 for(size_t rep=0; rep<1000; ++rep)
72 for(size_t i=0; i<(size_t)mat.rows(); ++i)
73 for(size_t j=0; j<(size_t)mat.cols(); ++j)
74 mat(i,j) = uniform(rng);
75
76 gttic_(basicTime);
77 for(size_t rep=0; rep<nReps; ++rep)
78 for(size_t i=0; i<(size_t)mat.rows(); ++i)
79 for(size_t j=0; j<(size_t)mat.cols(); ++j)
80 mat(i,j) = uniform(rng);
81 gttoc_(basicTime);
82 tictoc_getNode(basicTimeNode, basicTime);
83 basicTime = basicTimeNode->secs();
84 gtsam::tictoc_reset_();
85 cout << " Basic: " << double(1000000 * basicTime / double(mat.rows()*mat.cols()*nReps)) << " μs/element" << endl;
86
87 gttic_(fullTime);
88 for(size_t rep=0; rep<nReps; ++rep)
89 for(size_t i=0; i<(size_t)full.rows(); ++i)
90 for(size_t j=0; j<(size_t)full.cols(); ++j)
91 full(i,j) = uniform(rng);
92 gttoc_(fullTime);
93 tictoc_getNode(fullTimeNode, fullTime);
94 fullTime = fullTimeNode->secs();
95 gtsam::tictoc_reset_();
96 cout << " Full: " << double(1000000 * fullTime / double(full.rows()*full.cols()*nReps)) << " μs/element" << endl;
97
98 gttic_(topTime);
99 for(size_t rep=0; rep<nReps; ++rep)
100 for(size_t i=0; i<(size_t)top.rows(); ++i)
101 for(size_t j=0; j<(size_t)top.cols(); ++j)
102 top(i,j) = uniform(rng);
103 gttoc_(topTime);
104 tictoc_getNode(topTimeNode, topTime);
105 topTime = topTimeNode->secs();
106 gtsam::tictoc_reset_();
107 cout << " Top: " << double(1000000 * topTime / double(top.rows()*top.cols()*nReps)) << " μs/element" << endl;
108
109 gttic_(blockTime);
110 for(size_t rep=0; rep<nReps; ++rep)
111 for(size_t i=0; i<(size_t)block.rows(); ++i)
112 for(size_t j=0; j<(size_t)block.cols(); ++j)
113 block(i,j) = uniform(rng);
114 gttoc_(blockTime);
115 tictoc_getNode(blockTimeNode, blockTime);
116 blockTime = blockTimeNode->secs();
117 gtsam::tictoc_reset_();
118 cout << " Block: " << double(1000000 * blockTime / double(block.rows()*block.cols()*nReps)) << " μs/element" << endl;
119
120 cout << endl;
121 }
122
123 {
124 double basicTime, fullTime, topTime, blockTime;
125
126 cout << "Row-major matrix, column-major assignment:" << endl;
127
128 // Do a few initial assignments to let any cache effects stabilize
129 for(size_t rep=0; rep<1000; ++rep)
130 for(size_t j=0; j<(size_t)mat.cols(); ++j)
131 for(size_t i=0; i<(size_t)mat.rows(); ++i)
132 mat(i,j) = uniform(rng);
133
134 gttic_(basicTime);
135 for(size_t rep=0; rep<nReps; ++rep)
136 for(size_t j=0; j<(size_t)mat.cols(); ++j)
137 for(size_t i=0; i<(size_t)mat.rows(); ++i)
138 mat(i,j) = uniform(rng);
139 gttoc_(basicTime);
140 tictoc_getNode(basicTimeNode, basicTime);
141 basicTime = basicTimeNode->secs();
142 gtsam::tictoc_reset_();
143 cout << " Basic: " << double(1000000 * basicTime / double(mat.rows()*mat.cols()*nReps)) << " μs/element" << endl;
144
145 gttic_(fullTime);
146 for(size_t rep=0; rep<nReps; ++rep)
147 for(size_t j=0; j<(size_t)full.cols(); ++j)
148 for(size_t i=0; i<(size_t)full.rows(); ++i)
149 full(i,j) = uniform(rng);
150 gttoc_(fullTime);
151 tictoc_getNode(fullTimeNode, fullTime);
152 fullTime = fullTimeNode->secs();
153 gtsam::tictoc_reset_();
154 cout << " Full: " << double(1000000 * fullTime / double(full.rows()*full.cols()*nReps)) << " μs/element" << endl;
155
156 gttic_(topTime);
157 for(size_t rep=0; rep<nReps; ++rep)
158 for(size_t j=0; j<(size_t)top.cols(); ++j)
159 for(size_t i=0; i<(size_t)top.rows(); ++i)
160 top(i,j) = uniform(rng);
161 gttoc_(topTime);
162 tictoc_getNode(topTimeNode, topTime);
163 topTime = topTimeNode->secs();
164 gtsam::tictoc_reset_();
165 cout << " Top: " << double(1000000 * topTime / double(top.rows()*top.cols()*nReps)) << " μs/element" << endl;
166
167 gttic_(blockTime);
168 for(size_t rep=0; rep<nReps; ++rep)
169 for(size_t j=0; j<(size_t)block.cols(); ++j)
170 for(size_t i=0; i<(size_t)block.rows(); ++i)
171 block(i,j) = uniform(rng);
172 gttoc_(blockTime);
173 tictoc_getNode(blockTimeNode, blockTime);
174 blockTime = blockTimeNode->secs();
175 gtsam::tictoc_reset_();
176 cout << " Block: " << double(1000000 * blockTime / double(block.rows()*block.cols()*nReps)) << " μs/element" << endl;
177
178 cout << endl;
179 }
180
181 {
182 double basicTime, fullTime, topTime, blockTime;
183 typedef std::pair<size_t,size_t> ij_t;
184 std::vector<ij_t> ijs(100000);
185
186 cout << "Row-major matrix, random assignment:" << endl;
187
188 // Do a few initial assignments to let any cache effects stabilize
189 for (auto& ij : ijs) ij = {uniform_i(rng), uniform_j(rng)};
190 for(size_t rep=0; rep<1000; ++rep)
191 for(const auto& [i, j]: ijs) { mat(i, j) = uniform(rng); }
192
193 gttic_(basicTime);
194 for (auto& ij : ijs) ij = {uniform_i(rng), uniform_j(rng)};
195 for(size_t rep=0; rep<1000; ++rep)
196 for(const auto& [i, j]: ijs) { mat(i, j) = uniform(rng); }
197 gttoc_(basicTime);
198 tictoc_getNode(basicTimeNode, basicTime);
199 basicTime = basicTimeNode->secs();
200 gtsam::tictoc_reset_();
201 cout << " Basic: " << double(1000000 * basicTime / double(ijs.size()*nReps)) << " μs/element" << endl;
202
203 gttic_(fullTime);
204 for (auto& ij : ijs) ij = {uniform_i(rng), uniform_j(rng)};
205 for(size_t rep=0; rep<1000; ++rep)
206 for(const auto& [i, j]: ijs) { full(i, j) = uniform(rng); }
207 gttoc_(fullTime);
208 tictoc_getNode(fullTimeNode, fullTime);
209 fullTime = fullTimeNode->secs();
210 gtsam::tictoc_reset_();
211 cout << " Full: " << double(1000000 * fullTime / double(ijs.size()*nReps)) << " μs/element" << endl;
212
213 gttic_(topTime);
214 for (auto& ij : ijs) ij = {uniform_i(rng) % top.rows(), uniform_j(rng)};
215 for(size_t rep=0; rep<1000; ++rep)
216 for(const auto& [i, j]: ijs) { top(i, j) = uniform(rng); }
217 gttoc_(topTime);
218 tictoc_getNode(topTimeNode, topTime);
219 topTime = topTimeNode->secs();
220 gtsam::tictoc_reset_();
221 cout << " Top: " << double(1000000 * topTime / double(ijs.size()*nReps)) << " μs/element" << endl;
222
223 gttic_(blockTime);
224 for (auto& ij : ijs)
225 ij = {uniform_i(rng) % block.rows(), uniform_j(rng) % block.cols()};
226 for(size_t rep=0; rep<1000; ++rep)
227 for(const auto& [i, j]: ijs) { block(i, j) = uniform(rng); }
228 gttoc_(blockTime);
229 tictoc_getNode(blockTimeNode, blockTime);
230 blockTime = blockTimeNode->secs();
231 gtsam::tictoc_reset_();
232 cout << " Block: " << double(1000000 * blockTime / double(ijs.size()*nReps)) << " μs/element" << endl;
233
234 cout << endl;
235 }
236 }
237
238// if(true) {
239// cout << "\nTesting square triangular matrices:" << endl;
240//
241//// typedef triangular_matrix<double, ublas::upper, ublas::column_major> triangular;
242//// typedef ublas::matrix<double, ublas::column_major> matrix;
243// typedef MatrixXd matrix; // default col major
244//
245//// triangular tri(5,5);
246//
247// matrix mat(5,5);
248// for(size_t j=0; j<(size_t)mat.cols(); ++j)
249// for(size_t i=0; i<(size_t)mat.rows(); ++i)
250// mat(i,j) = uniform(rng);
251//
252// tri = ublas::triangular_adaptor<matrix, ublas::upper>(mat);
253// cout << " Assigned from triangular adapter: " << tri << endl;
254//
255// cout << " Triangular adapter of mat: " << ublas::triangular_adaptor<matrix, ublas::upper>(mat) << endl;
256//
257// for(size_t j=0; j<(size_t)mat.cols(); ++j)
258// for(size_t i=0; i<(size_t)mat.rows(); ++i)
259// mat(i,j) = uniform(rng);
260// mat = tri;
261// cout << " Assign matrix from triangular: " << mat << endl;
262//
263// for(size_t j=0; j<(size_t)mat.cols(); ++j)
264// for(size_t i=0; i<(size_t)mat.rows(); ++i)
265// mat(i,j) = uniform(rng);
266// (ublas::triangular_adaptor<matrix, ublas::upper>(mat)) = tri;
267// cout << " Assign triangular adaptor from triangular: " << mat << endl;
268// }
269
270// {
271// cout << "\nTesting wide triangular matrices:" << endl;
272//
273// typedef triangular_matrix<double, ublas::upper, ublas::column_major> triangular;
274// typedef ublas::matrix<double, ublas::column_major> matrix;
275//
276// triangular tri(5,7);
277//
278// matrix mat(5,7);
279// for(size_t j=0; j<(size_t)mat.cols(); ++j)
280// for(size_t i=0; i<(size_t)mat.rows(); ++i)
281// mat(i,j) = uniform(rng);
282//
283// tri = ublas::triangular_adaptor<matrix, ublas::upper>(mat);
284// cout << " Assigned from triangular adapter: " << tri << endl;
285//
286// cout << " Triangular adapter of mat: " << ublas::triangular_adaptor<matrix, ublas::upper>(mat) << endl;
287//
288// for(size_t j=0; j<(size_t)mat.cols(); ++j)
289// for(size_t i=0; i<(size_t)mat.rows(); ++i)
290// mat(i,j) = uniform(rng);
291// mat = tri;
292// cout << " Assign matrix from triangular: " << mat << endl;
293//
294// for(size_t j=0; j<(size_t)mat.cols(); ++j)
295// for(size_t i=0; i<(size_t)mat.rows(); ++i)
296// mat(i,j) = uniform(rng);
297// mat = ublas::triangular_adaptor<matrix, ublas::upper>(mat);
298// cout << " Assign matrix from triangular adaptor of self: " << mat << endl;
299// }
300
301// {
302// cout << "\nTesting subvectors:" << endl;
303//
304// typedef MatrixXd matrix;
305// matrix mat(4,4);
306//
307// for(size_t j=0; j<(size_t)mat.cols(); ++j)
308// for(size_t i=0; i<(size_t)mat.rows(); ++i)
309// mat(i,j) = i*mat.rows() + j;
310// cout << " mat = " << mat;
311//
312// cout << " vec(1:4, 2:2) = " << mat.block(1,2, ), ublas::range(1,4), ublas::range(2,2));
313//
314// }
315
316 return 0;
317
318}
319