timeMatrixOps.cpp
Go to the documentation of this file.
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 
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 
27 using namespace std;
28 //namespace ublas = boost::numeric::ublas;
29 //using namespace Eigen;
30 
31 static std::mt19937 rng;
32 static 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 
41 int 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(0, 0, m, n);
56  gtsam::SubMatrix top = mat.block(0, 0, n, n);
57  gtsam::SubMatrix block = mat.block(m/4, n/4, m-m/2, 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();
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();
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();
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();
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();
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();
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();
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();
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();
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();
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();
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();
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 //
243 // typedef MatrixXd matrix; // default col major
244 //
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 }
timing.h
Timing utilities.
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
main
int main(int argc, char *argv[])
Definition: timeMatrixOps.cpp:41
Eigen::Block
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
gtsam::tictoc_reset_
void tictoc_reset_()
Definition: timing.h:282
Matrix.h
typedef and functions to augment Eigen's MatrixXd
gtsam::Matrix
Eigen::MatrixXd Matrix
Definition: base/Matrix.h:39
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
block
m m block(1, 0, 2, 2)<< 4
n
int n
Definition: BiCGSTAB_simple.cpp:1
rng
static std::mt19937 rng
Definition: timeMatrixOps.cpp:31
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
gttoc_
#define gttoc_(label)
Definition: timing.h:250
gttic_
#define gttic_(label)
Definition: timing.h:245
tictoc_getNode
#define tictoc_getNode(variable, label)
Definition: timing.h:276
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
size_t
std::size_t size_t
Definition: wrap/pybind11/include/pybind11/detail/common.h:490
std
Definition: BFloat16.h:88
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
uniform
static std::uniform_real_distribution uniform(-1.0, 0.0)


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:42:09