Visual Servoing Platform  version 3.4.0
perfMatrixMultiplication.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2019 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See http://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31  * Description:
32  * Benchmark matrix multiplication.
33  *
34  *****************************************************************************/
35 
36 #include <visp3/core/vpConfig.h>
37 
38 #ifdef VISP_HAVE_CATCH2
39 #define CATCH_CONFIG_ENABLE_BENCHMARKING
40 #define CATCH_CONFIG_RUNNER
41 #include <catch.hpp>
42 
43 #include <visp3/core/vpMatrix.h>
44 
45 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
46 #include <opencv2/core.hpp>
47 #endif
48 
49 #ifdef VISP_HAVE_EIGEN3
50 #include <Eigen/Dense>
51 #endif
52 
53 namespace
54 {
55 
56 bool runBenchmark = false;
57 bool runBenchmarkAll = false;
58 
59 double getRandomValues(double min, double max)
60 {
61  return (max - min) * ((double)rand() / (double)RAND_MAX) + min;
62 }
63 
64 vpMatrix generateRandomMatrix(unsigned int rows, unsigned int cols, double min=-1, double max=1)
65 {
66  vpMatrix M(rows, cols);
67 
68  for (unsigned int i = 0; i < M.getRows(); i++) {
69  for (unsigned int j = 0; j < M.getCols(); j++) {
70  M[i][j] = getRandomValues(min, max);
71  }
72  }
73 
74  return M;
75 }
76 
77 vpColVector generateRandomVector(unsigned int rows, double min=-1, double max=1)
78 {
79  vpColVector v(rows);
80 
81  for (unsigned int i = 0; i < v.getRows(); i++) {
82  v[i] = getRandomValues(min, max);
83  }
84 
85  return v;
86 }
87 
88 // Copy of vpMatrix::mult2Matrices
89 vpMatrix dgemm_regular(const vpMatrix& A, const vpMatrix& B)
90 {
91  vpMatrix C;
92 
93  if ((A.getRows() != C.getRows()) || (B.getCols() != C.getCols())) {
94  C.resize(A.getRows(), B.getCols(), false);
95  }
96 
97  if (A.getCols() != B.getRows()) {
98  throw(vpException(vpException::dimensionError, "Cannot multiply (%dx%d) matrix by (%dx%d) matrix", A.getRows(),
99  A.getCols(), B.getRows(), B.getCols()));
100  }
101 
102  // 5/12/06 some "very" simple optimization to avoid indexation
103  unsigned int BcolNum = B.getCols();
104  unsigned int BrowNum = B.getRows();
105  unsigned int i, j, k;
106  for (i = 0; i < A.getRows(); i++) {
107  double *rowptri = A[i];
108  double *ci = C[i];
109  for (j = 0; j < BcolNum; j++) {
110  double s = 0;
111  for (k = 0; k < BrowNum; k++) {
112  s += rowptri[k] * B[k][j];
113  }
114  ci[j] = s;
115  }
116  }
117 
118  return C;
119 }
120 
121 // Copy of vpMatrix::AtA
122 vpMatrix AtA_regular(const vpMatrix& A)
123 {
124  vpMatrix B;
125  B.resize(A.getCols(), A.getCols(), false);
126 
127  for (unsigned int i = 0; i < A.getCols(); i++) {
128  double *Bi = B[i];
129  for (unsigned int j = 0; j < i; j++) {
130  double *ptr = A.data;
131  double s = 0;
132  for (unsigned int k = 0; k < A.getRows(); k++) {
133  s += (*(ptr + i)) * (*(ptr + j));
134  ptr += A.getCols();
135  }
136  *Bi++ = s;
137  B[j][i] = s;
138  }
139  double *ptr = A.data;
140  double s = 0;
141  for (unsigned int k = 0; k < A.getRows(); k++) {
142  s += (*(ptr + i)) * (*(ptr + i));
143  ptr += A.getCols();
144  }
145  *Bi = s;
146  }
147 
148  return B;
149 }
150 
151 // Copy of vpMatrix::AAt()
152 vpMatrix AAt_regular(const vpMatrix& A)
153 {
154  vpMatrix B;
155  B.resize(A.getRows(), A.getRows(), false);
156 
157  // compute A*A^T
158  for (unsigned int i = 0; i < A.getRows(); i++) {
159  for (unsigned int j = i; j < A.getRows(); j++) {
160  double *pi = A[i]; // row i
161  double *pj = A[j]; // row j
162 
163  // sum (row i .* row j)
164  double ssum = 0;
165  for (unsigned int k = 0; k < A.getCols(); k++)
166  ssum += *(pi++) * *(pj++);
167 
168  B[i][j] = ssum; // upper triangle
169  if (i != j)
170  B[j][i] = ssum; // lower triangle
171  }
172  }
173  return B;
174 }
175 
176 // Copy of vpMatrix::multMatrixVector
177 vpColVector dgemv_regular(const vpMatrix& A, const vpColVector& v)
178 {
179  vpColVector w;
180 
181  if (A.getCols() != v.getRows()) {
182  throw(vpException(vpException::dimensionError, "Cannot multiply a (%dx%d) matrix by a (%d) column vector",
183  A.getRows(), A.getCols(), v.getRows()));
184  }
185 
186  w.resize(A.getRows(), true);
187 
188  for (unsigned int j = 0; j < A.getCols(); j++) {
189  double vj = v[j]; // optimization em 5/12/2006
190  for (unsigned int i = 0; i < A.getRows(); i++) {
191  w[i] += A[i][j] * vj;
192  }
193  }
194 
195  return w;
196 }
197 
198 bool equalMatrix(const vpMatrix& A, const vpMatrix& B, double tol=1e-9)
199 {
200  if (A.getRows() != B.getRows() || A.getCols() != B.getCols()) {
201  return false;
202  }
203 
204  for (unsigned int i = 0; i < A.getRows(); i++) {
205  for (unsigned int j = 0; j < A.getCols(); j++) {
206  if (!vpMath::equal(A[i][j], B[i][j], tol)) {
207  return false;
208  }
209  }
210  }
211 
212  return true;
213 }
214 
215 }
216 
217 TEST_CASE("Benchmark matrix-matrix multiplication", "[benchmark]") {
218  if (runBenchmark || runBenchmarkAll) {
219  std::vector<std::pair<int, int>> sizes = { {3, 3}, {6, 6}, {8, 8}, {10, 10}, {20, 20}, {6, 200}, {200, 6}, {207, 119}, {83, 201}, {600, 400}, {400, 600} };
220 
221  for (auto sz : sizes) {
222  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
223  vpMatrix B = generateRandomMatrix(sz.second, sz.first);
224 
225  std::ostringstream oss;
226  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - Naive code";
227  vpMatrix C, C_true;
228  BENCHMARK(oss.str().c_str()) {
229  C_true = dgemm_regular(A, B);
230  return C_true;
231  };
232 
233  oss.str("");
234  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - ViSP";
235  BENCHMARK(oss.str().c_str()) {
236  C = A * B;
237  return C;
238  };
239  REQUIRE(equalMatrix(C, C_true));
240 
241  if(runBenchmarkAll) {
242 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
243  cv::Mat matA(sz.first, sz.second, CV_64FC1);
244  cv::Mat matB(sz.second, sz.first, CV_64FC1);
245 
246  for (unsigned int i = 0; i < A.getRows(); i++) {
247  for (unsigned int j = 0; j < A.getCols(); j++) {
248  matA.at<double>(i, j) = A[i][j];
249  matB.at<double>(j, i) = B[j][i];
250  }
251  }
252 
253  oss.str("");
254  oss << "(" << matA.rows << "x" << matA.cols << ")x(" << matB.rows << "x" << matB.cols << ") - OpenCV";
255  BENCHMARK(oss.str().c_str()) {
256  cv::Mat matC = matA * matB;
257  return matC;
258  };
259 #endif
260 
261 #ifdef VISP_HAVE_EIGEN3
262  Eigen::MatrixXd eigenA(sz.first, sz.second);
263  Eigen::MatrixXd eigenB(sz.second, sz.first);
264 
265  for (unsigned int i = 0; i < A.getRows(); i++) {
266  for (unsigned int j = 0; j < A.getCols(); j++) {
267  eigenA(i, j) = A[i][j];
268  eigenB(j, i) = B[j][i];
269  }
270  }
271 
272  oss.str("");
273  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(" << eigenB.rows() << "x" << eigenB.cols() << ") - Eigen";
274  BENCHMARK(oss.str().c_str()) {
275  Eigen::MatrixXd eigenC = eigenA * eigenB;
276  return eigenC;
277  };
278 #endif
279  }
280  }
281  }
282 
283  {
284  const unsigned int rows = 47, cols = 63;
285  vpMatrix A = generateRandomMatrix(rows, cols);
286  vpMatrix B = generateRandomMatrix(cols, rows);
287 
288  vpMatrix C_true = dgemm_regular(A, B);
289  vpMatrix C = A * B;
290  REQUIRE(equalMatrix(C, C_true));
291  }
292 }
293 
294 TEST_CASE("Benchmark matrix-rotation matrix multiplication", "[benchmark]") {
295  if (runBenchmark || runBenchmarkAll) {
296  std::vector<std::pair<int, int>> sizes = { {3, 3} };
297 
298  for (auto sz : sizes) {
299  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
300  vpRotationMatrix B(vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)));
301 
302  std::ostringstream oss;
303  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - Naive code";
304  vpMatrix AB, AB_true;
305  BENCHMARK(oss.str().c_str()) {
306  AB_true = dgemm_regular(A, B);
307  return AB_true;
308  };
309 
310  oss.str("");
311  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - ViSP";
312  BENCHMARK(oss.str().c_str()) {
313  AB = A * B;
314  return AB;
315  };
316  REQUIRE(equalMatrix(AB, AB_true));
317 
318  if(runBenchmarkAll) {
319 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
320  cv::Mat matA(sz.first, sz.second, CV_64FC1);
321  cv::Mat matB(3, 3, CV_64FC1);
322 
323  for (unsigned int i = 0; i < A.getRows(); i++) {
324  for (unsigned int j = 0; j < A.getCols(); j++) {
325  matA.at<double>(i, j) = A[i][j];
326  }
327  }
328  for (unsigned int i = 0; i < B.getRows(); i++) {
329  for (unsigned int j = 0; j < B.getCols(); j++) {
330  matB.at<double>(j, i) = B[j][i];
331  }
332  }
333 
334  oss.str("");
335  oss << "(" << matA.rows << "x" << matA.cols << ")x(" << matB.rows << "x" << matB.cols << ") - OpenCV";
336  BENCHMARK(oss.str().c_str()) {
337  cv::Mat matC = matA * matB;
338  return matC;
339  };
340 #endif
341 
342 #ifdef VISP_HAVE_EIGEN3
343  Eigen::MatrixXd eigenA(sz.first, sz.second);
344  Eigen::MatrixXd eigenB(3, 3);
345 
346  for (unsigned int i = 0; i < A.getRows(); i++) {
347  for (unsigned int j = 0; j < A.getCols(); j++) {
348  eigenA(i, j) = A[i][j];
349  }
350  }
351  for (unsigned int i = 0; i < B.getRows(); i++) {
352  for (unsigned int j = 0; j < B.getCols(); j++) {
353  eigenB(j, i) = B[j][i];
354  }
355  }
356 
357  oss.str("");
358  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(" << eigenB.rows() << "x" << eigenB.cols() << ") - Eigen";
359  BENCHMARK(oss.str().c_str()) {
360  Eigen::MatrixXd eigenC = eigenA * eigenB;
361  return eigenC;
362  };
363 #endif
364  }
365  }
366  }
367 
368  {
369  const unsigned int rows = 3, cols = 3;
370  vpMatrix A = generateRandomMatrix(rows, cols);
371  vpRotationMatrix B(vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)));
372 
373  vpMatrix AB_true = dgemm_regular(A, B);
374  vpMatrix AB = A * B;
375  REQUIRE(equalMatrix(AB, AB_true));
376  }
377 }
378 
379 TEST_CASE("Benchmark matrix-homogeneous matrix multiplication", "[benchmark]") {
380  if (runBenchmark || runBenchmarkAll) {
381  std::vector<std::pair<int, int>> sizes = { {4, 4} };
382 
383  for (auto sz : sizes) {
384  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
385  vpHomogeneousMatrix B(getRandomValues(0, 1), getRandomValues(0, 1), getRandomValues(0, 1),
386  vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)));
387 
388  std::ostringstream oss;
389  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - Naive code";
390  vpMatrix AB, AB_true;
391  BENCHMARK(oss.str().c_str()) {
392  AB_true = dgemm_regular(A, B);
393  return AB_true;
394  };
395 
396  oss.str("");
397  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - ViSP";
398  BENCHMARK(oss.str().c_str()) {
399  AB = A * B;
400  return AB;
401  };
402  REQUIRE(equalMatrix(AB, AB_true));
403 
404  if(runBenchmarkAll) {
405 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
406  cv::Mat matA(sz.first, sz.second, CV_64FC1);
407  cv::Mat matB(4, 4, CV_64FC1);
408 
409  for (unsigned int i = 0; i < A.getRows(); i++) {
410  for (unsigned int j = 0; j < A.getCols(); j++) {
411  matA.at<double>(i, j) = A[i][j];
412  }
413  }
414  for (unsigned int i = 0; i < B.getRows(); i++) {
415  for (unsigned int j = 0; j < B.getCols(); j++) {
416  matB.at<double>(j, i) = B[j][i];
417  }
418  }
419 
420  oss.str("");
421  oss << "(" << matA.rows << "x" << matA.cols << ")x(" << matB.rows << "x" << matB.cols << ") - OpenCV";
422  BENCHMARK(oss.str().c_str()) {
423  cv::Mat matC = matA * matB;
424  return matC;
425  };
426 #endif
427 
428 #ifdef VISP_HAVE_EIGEN3
429  Eigen::MatrixXd eigenA(sz.first, sz.second);
430  Eigen::MatrixXd eigenB(4, 4);
431 
432  for (unsigned int i = 0; i < A.getRows(); i++) {
433  for (unsigned int j = 0; j < A.getCols(); j++) {
434  eigenA(i, j) = A[i][j];
435  }
436  }
437  for (unsigned int i = 0; i < B.getRows(); i++) {
438  for (unsigned int j = 0; j < B.getCols(); j++) {
439  eigenB(j, i) = B[j][i];
440  }
441  }
442 
443  oss.str("");
444  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(" << eigenB.rows() << "x" << eigenB.cols() << ") - Eigen";
445  BENCHMARK(oss.str().c_str()) {
446  Eigen::MatrixXd eigenC = eigenA * eigenB;
447  return eigenC;
448  };
449 #endif
450  }
451  }
452  }
453 
454  {
455  const unsigned int rows = 4, cols = 4;
456  vpMatrix A = generateRandomMatrix(rows, cols);
457  vpHomogeneousMatrix B(getRandomValues(0, 1), getRandomValues(0, 1), getRandomValues(0, 1),
458  vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)));
459 
460  vpMatrix AB_true = dgemm_regular(A, B);
461  vpMatrix AB;
462  vpMatrix::mult2Matrices(A, B, AB);
463  REQUIRE(equalMatrix(AB, AB_true));
464  }
465 }
466 
467 TEST_CASE("Benchmark matrix-vector multiplication", "[benchmark]") {
468  if (runBenchmark || runBenchmarkAll) {
469  std::vector<std::pair<int, int>> sizes = { {3, 3}, {6, 6}, {8, 8}, {10, 10}, {20, 20}, {6, 200}, {200, 6}, {207, 119}, {83, 201}, {600, 400}, {400, 600} };
470 
471  for (auto sz : sizes) {
472  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
473  vpColVector B = generateRandomVector(sz.second);
474 
475  std::ostringstream oss;
476  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - Naive code";
477  vpColVector C, C_true;
478  BENCHMARK(oss.str().c_str()) {
479  C_true = dgemv_regular(A, B);
480  return C_true;
481  };
482 
483  oss.str("");
484  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - ViSP";
485  BENCHMARK(oss.str().c_str()) {
486  C = A * B;
487  return C;
488  };
489  REQUIRE(equalMatrix(C, C_true));
490 
491  if(runBenchmarkAll) {
492 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
493  cv::Mat matA(sz.first, sz.second, CV_64FC1);
494  cv::Mat matB(sz.second, 1, CV_64FC1);
495 
496  for (unsigned int i = 0; i < A.getRows(); i++) {
497  for (unsigned int j = 0; j < A.getCols(); j++) {
498  matA.at<double>(i, j) = A[i][j];
499  if (i == 0) {
500  matB.at<double>(j, 0) = B[j];
501  }
502  }
503  }
504 
505  oss.str("");
506  oss << "(" << matA.rows << "x" << matA.cols << ")x(" << matB.rows << "x" << matB.cols << ") - OpenCV";
507  BENCHMARK(oss.str().c_str()) {
508  cv::Mat matC = matA * matB;
509  return matC;
510  };
511 #endif
512 
513 #ifdef VISP_HAVE_EIGEN3
514  Eigen::MatrixXd eigenA(sz.first, sz.second);
515  Eigen::MatrixXd eigenB(sz.second, 1);
516 
517  for (unsigned int i = 0; i < A.getRows(); i++) {
518  for (unsigned int j = 0; j < A.getCols(); j++) {
519  eigenA(i, j) = A[i][j];
520  if (i == 0) {
521  eigenB(j, 0) = B[j];
522  }
523  }
524  }
525 
526  oss.str("");
527  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(" << eigenB.rows() << "x" << eigenB.cols() << ") - Eigen";
528  BENCHMARK(oss.str().c_str()) {
529  Eigen::MatrixXd eigenC = eigenA * eigenB;
530  return eigenC;
531  };
532 #endif
533  }
534  }
535  }
536 
537  {
538  const unsigned int rows = 47, cols = 63;
539  vpMatrix A = generateRandomMatrix(rows, cols);
540  vpColVector B = generateRandomVector(cols);
541 
542  vpColVector C_true = dgemv_regular(A, B);
543  vpColVector C = A * B;
544  REQUIRE(equalMatrix(C, C_true));
545  }
546 }
547 
548 TEST_CASE("Benchmark AtA", "[benchmark]") {
549  if (runBenchmark || runBenchmarkAll) {
550  std::vector<std::pair<int, int>> sizes = { {3, 3}, {6, 6}, {8, 8}, {10, 10}, {20, 20}, {6, 200}, {200, 6}, {207, 119}, {83, 201}, {600, 400}, {400, 600} };
551 
552  for (auto sz : sizes) {
553  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
554 
555  std::ostringstream oss;
556  oss << "(" << A.getRows() << "x" << A.getCols() << ") - Naive code";
557  vpMatrix AtA, AtA_true;
558  BENCHMARK(oss.str().c_str()) {
559  AtA_true = AtA_regular(A);
560  return AtA_true;
561  };
562 
563  oss.str("");
564  oss << "(" << A.getRows() << "x" << A.getCols() << ") - ViSP";
565  BENCHMARK(oss.str().c_str()) {
566  AtA = A.AtA();
567  return AtA;
568  };
569  REQUIRE(equalMatrix(AtA, AtA_true));
570 
571  if(runBenchmarkAll) {
572 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
573  cv::Mat matA(sz.first, sz.second, CV_64FC1);
574 
575  for (unsigned int i = 0; i < A.getRows(); i++) {
576  for (unsigned int j = 0; j < A.getCols(); j++) {
577  matA.at<double>(i, j) = A[i][j];
578  }
579  }
580 
581  oss.str("");
582  oss << "(" << matA.rows << "x" << matA.cols << ") - OpenCV";
583  BENCHMARK(oss.str().c_str()) {
584  cv::Mat matAtA = matA.t() * matA;
585  return matAtA;
586  };
587 #endif
588 
589 #ifdef VISP_HAVE_EIGEN3
590  Eigen::MatrixXd eigenA(sz.first, sz.second);
591 
592  for (unsigned int i = 0; i < A.getRows(); i++) {
593  for (unsigned int j = 0; j < A.getCols(); j++) {
594  eigenA(i, j) = A[i][j];
595  }
596  }
597 
598  oss.str("");
599  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ") - Eigen";
600  BENCHMARK(oss.str().c_str()) {
601  Eigen::MatrixXd eigenAtA = eigenA.transpose() * eigenA;
602  return eigenAtA;
603  };
604 #endif
605  }
606  }
607  }
608 
609  {
610  const unsigned int rows = 47, cols = 63;
611  vpMatrix A = generateRandomMatrix(rows, cols);
612 
613  vpMatrix AtA_true = AtA_regular(A);
614  vpMatrix AtA = A.AtA();
615  REQUIRE(equalMatrix(AtA, AtA_true));
616  }
617 }
618 
619 TEST_CASE("Benchmark AAt", "[benchmark]") {
620  if (runBenchmark || runBenchmarkAll) {
621  std::vector<std::pair<int, int>> sizes = { {3, 3}, {6, 6}, {8, 8}, {10, 10}, {20, 20}, {6, 200}, {200, 6} };//, {207, 119}, {83, 201}, {600, 400}, {400, 600} };
622 
623  for (auto sz : sizes) {
624  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
625 
626  std::ostringstream oss;
627  oss << "(" << A.getRows() << "x" << A.getCols() << ") - Naive code";
628  vpMatrix AAt_true, AAt;
629  BENCHMARK(oss.str().c_str()) {
630  AAt_true = AAt_regular(A);
631  return AAt_true;
632  };
633 
634  oss.str("");
635  oss << "(" << A.getRows() << "x" << A.getCols() << ") - ViSP";
636  BENCHMARK(oss.str().c_str()) {
637  AAt = A.AAt();
638  return AAt;
639  };
640  REQUIRE(equalMatrix(AAt, AAt_true));
641 
642  if(runBenchmarkAll) {
643 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
644  cv::Mat matA(sz.first, sz.second, CV_64FC1);
645 
646  for (unsigned int i = 0; i < A.getRows(); i++) {
647  for (unsigned int j = 0; j < A.getCols(); j++) {
648  matA.at<double>(i, j) = A[i][j];
649  }
650  }
651 
652  oss.str("");
653  oss << "(" << matA.rows << "x" << matA.cols << ") - OpenCV";
654  BENCHMARK(oss.str().c_str()) {
655  cv::Mat matAAt = matA * matA.t();
656  return matAAt;
657  };
658 #endif
659 
660 #ifdef VISP_HAVE_EIGEN3
661  Eigen::MatrixXd eigenA(sz.first, sz.second);
662 
663  for (unsigned int i = 0; i < A.getRows(); i++) {
664  for (unsigned int j = 0; j < A.getCols(); j++) {
665  eigenA(i, j) = A[i][j];
666  }
667  }
668 
669  oss.str("");
670  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ") - Eigen";
671  BENCHMARK(oss.str().c_str()) {
672  Eigen::MatrixXd eigenAAt = eigenA * eigenA.transpose();
673  return eigenAAt;
674  };
675 #endif
676  }
677  }
678  }
679 
680  {
681  const unsigned int rows = 47, cols = 63;
682  vpMatrix A = generateRandomMatrix(rows, cols);
683 
684  vpMatrix AAt_true = AAt_regular(A);
685  vpMatrix AAt = A.AAt();
686  REQUIRE(equalMatrix(AAt, AAt_true));
687  }
688 }
689 
690 TEST_CASE("Benchmark matrix-velocity twist multiplication", "[benchmark]") {
691  if (runBenchmark || runBenchmarkAll) {
692  std::vector<std::pair<int, int>> sizes = { {6, 6}, {20, 6}, {207, 6}, {600, 6}, {1201, 6} };
693 
694  for (auto sz : sizes) {
695  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
696  vpVelocityTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
697 
698  std::ostringstream oss;
699  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - Naive code";
700  vpMatrix AV, AV_true;
701  BENCHMARK(oss.str().c_str()) {
702  AV_true = dgemm_regular(A, V);
703  return AV_true;
704  };
705 
706  oss.str("");
707  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - ViSP";
708  BENCHMARK(oss.str().c_str()) {
709  AV = A * V;
710  return AV;
711  };
712  REQUIRE(equalMatrix(AV, AV_true));
713 
714  if(runBenchmarkAll) {
715 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
716  cv::Mat matA(sz.first, sz.second, CV_64FC1);
717  cv::Mat matV(6, 6, CV_64FC1);
718 
719  for (unsigned int i = 0; i < A.getRows(); i++) {
720  for (unsigned int j = 0; j < A.getCols(); j++) {
721  matA.at<double>(i, j) = A[i][j];
722  }
723  }
724  for (unsigned int i = 0; i < V.getRows(); i++) {
725  for (unsigned int j = 0; j < V.getCols(); j++) {
726  matV.at<double>(i, j) = V[i][j];
727  }
728  }
729 
730  oss.str("");
731  oss << "(" << matA.rows << "x" << matA.cols << ")x(6x6) - OpenCV";
732  BENCHMARK(oss.str().c_str()) {
733  cv::Mat matAV = matA * matV;
734  return matAV;
735  };
736 #endif
737 
738 #ifdef VISP_HAVE_EIGEN3
739  Eigen::MatrixXd eigenA(sz.first, sz.second);
740  Eigen::MatrixXd eigenV(6, 6);
741 
742  for (unsigned int i = 0; i < A.getRows(); i++) {
743  for (unsigned int j = 0; j < A.getCols(); j++) {
744  eigenA(i, j) = A[i][j];
745  }
746  }
747  for (unsigned int i = 0; i < V.getRows(); i++) {
748  for (unsigned int j = 0; j < V.getCols(); j++) {
749  eigenV(i, j) = V[i][j];
750  }
751  }
752 
753  oss.str("");
754  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(6x6) - Eigen";
755  BENCHMARK(oss.str().c_str()) {
756  Eigen::MatrixXd eigenAV = eigenA * eigenV;
757  return eigenAV;
758  };
759 #endif
760  }
761  }
762  }
763 
764  {
765  const unsigned int rows = 47, cols = 6;
766  vpMatrix A = generateRandomMatrix(rows, cols);
767  vpVelocityTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
768 
769  vpMatrix AV_true = dgemm_regular(A, V);
770  vpMatrix AV = A * V;
771  REQUIRE(equalMatrix(AV, AV_true));
772  }
773 }
774 
775 TEST_CASE("Benchmark matrix-force twist multiplication", "[benchmark]") {
776  if (runBenchmark || runBenchmarkAll) {
777  std::vector<std::pair<int, int>> sizes = { {6, 6}, {20, 6}, {207, 6}, {600, 6}, {1201, 6} };
778 
779  for (auto sz : sizes) {
780  vpMatrix A = generateRandomMatrix(sz.first, sz.second);
781  vpForceTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
782 
783  std::ostringstream oss;
784  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - Naive code";
785  vpMatrix AV, AV_true;
786  BENCHMARK(oss.str().c_str()) {
787  AV_true = dgemm_regular(A, V);
788  return AV_true;
789  };
790 
791  oss.str("");
792  oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - ViSP";
793  BENCHMARK(oss.str().c_str()) {
794  AV = A * V;
795  return AV;
796  };
797  REQUIRE(equalMatrix(AV, AV_true));
798 
799  if(runBenchmarkAll) {
800 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000)
801  cv::Mat matA(sz.first, sz.second, CV_64FC1);
802  cv::Mat matV(6, 6, CV_64FC1);
803 
804  for (unsigned int i = 0; i < A.getRows(); i++) {
805  for (unsigned int j = 0; j < A.getCols(); j++) {
806  matA.at<double>(i, j) = A[i][j];
807  }
808  }
809  for (unsigned int i = 0; i < V.getRows(); i++) {
810  for (unsigned int j = 0; j < V.getCols(); j++) {
811  matV.at<double>(i, j) = V[i][j];
812  }
813  }
814 
815  oss.str("");
816  oss << "(" << matA.rows << "x" << matA.cols << ")x(6x6) - OpenCV";
817  BENCHMARK(oss.str().c_str()) {
818  cv::Mat matAV = matA * matV;
819  return matAV;
820  };
821 #endif
822 
823 #ifdef VISP_HAVE_EIGEN3
824  Eigen::MatrixXd eigenA(sz.first, sz.second);
825  Eigen::MatrixXd eigenV(6, 6);
826 
827  for (unsigned int i = 0; i < A.getRows(); i++) {
828  for (unsigned int j = 0; j < A.getCols(); j++) {
829  eigenA(i, j) = A[i][j];
830  }
831  }
832  for (unsigned int i = 0; i < V.getRows(); i++) {
833  for (unsigned int j = 0; j < V.getCols(); j++) {
834  eigenV(i, j) = V[i][j];
835  }
836  }
837 
838  oss.str("");
839  oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(6x6) - Eigen";
840  BENCHMARK(oss.str().c_str()) {
841  Eigen::MatrixXd eigenAV = eigenA * eigenV;
842  return eigenAV;
843  };
844 #endif
845  }
846  }
847  }
848 
849  {
850  const unsigned int rows = 47, cols = 6;
851  vpMatrix A = generateRandomMatrix(rows, cols);
852  vpForceTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
853 
854  vpMatrix AV_true = dgemm_regular(A, V);
855  vpMatrix AV = A * V;
856  REQUIRE(equalMatrix(AV, AV_true));
857  }
858 }
859 
860 int main(int argc, char *argv[])
861 {
862  // Set random seed explicitly to avoid confusion
863  // See: https://en.cppreference.com/w/cpp/numeric/random/srand
864  // If rand() is used before any calls to srand(), rand() behaves as if it was seeded with srand(1).
865  srand(1);
866 
867  Catch::Session session; // There must be exactly one instance
868  unsigned int lapackMinSize = vpMatrix::getLapackMatrixMinSize();
869 
870  std::cout << "Default matrix/vector min size to enable Blas/Lapack optimization: "
871  << lapackMinSize << std::endl;
872  // Build a new parser on top of Catch's
873  using namespace Catch::clara;
874  auto cli = session.cli() // Get Catch's composite command line parser
875  | Opt(runBenchmark) // bind variable to a new option, with a hint string
876  ["--benchmark"] // the option names it will respond to
877  ("run benchmark comparing naive code with ViSP implementation") // description string for the help output
878  | Opt(runBenchmarkAll) // bind variable to a new option, with a hint string
879  ["--benchmark-all"] // the option names it will respond to
880  ("run benchmark comparing naive code with ViSP, OpenCV, Eigen implementation") // description string for the help output
881  | Opt(lapackMinSize, "min size") // bind variable to a new option, with a hint string
882  ["--lapack-min-size"] // the option names it will respond to
883  ("matrix/vector min size to enable blas/lapack usage"); // description string for the help output
884 
885  // Now pass the new composite back to Catch so it uses that
886  session.cli(cli);
887 
888  // Let Catch (using Clara) parse the command line
889  session.applyCommandLine(argc, argv);
890 
891  vpMatrix::setLapackMatrixMinSize(lapackMinSize);
892  std::cout << "Used matrix/vector min size to enable Blas/Lapack optimization: "
893  << vpMatrix::getLapackMatrixMinSize() << std::endl;
894 
895  int numFailed = session.run();
896 
897  // numFailed is clamped to 255 as some unices only use the lower 8 bits.
898  // This clamping has already been applied, so just return it here
899  // You can also do any post run clean-up here
900  return numFailed;
901 }
902 #else
903 #include <iostream>
904 
905 int main()
906 {
907  return 0;
908 }
909 #endif
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:153
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:304
Implementation of an homogeneous matrix and operations on such kind of matrices.
static bool equal(double x, double y, double s=0.001)
Definition: vpMath.h:293
error that can be emited by ViSP classes.
Definition: vpException.h:71
Type * data
Address of the first element of the data array.
Definition: vpArray2D.h:145
unsigned int getCols() const
Definition: vpArray2D.h:279
Implementation of a rotation matrix and operations on such kind of matrices.
static unsigned int getLapackMatrixMinSize()
Definition: vpMatrix.h:246
vpMatrix AtA() const
Definition: vpMatrix.cpp:629
vpMatrix AAt() const
Definition: vpMatrix.cpp:507
static void mult2Matrices(const vpMatrix &A, const vpMatrix &B, vpMatrix &C)
Definition: vpMatrix.cpp:1009
vpMatrix transpose() const
Definition: vpMatrix.cpp:474
unsigned int getRows() const
Definition: vpArray2D.h:289
static void setLapackMatrixMinSize(unsigned int min_size)
Definition: vpMatrix.h:262
vp_deprecated vpColVector rows(unsigned int first_row, unsigned int last_row) const
Definition: vpColVector.h:381
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:310
vpMatrix t() const
Definition: vpMatrix.cpp:464
static double deg(double rad)
Definition: vpMath.h:103
Implementation of column vector and the associated operations.
Definition: vpColVector.h:130
Class that consider the case of a translation vector.
Implementation of a rotation vector as axis-angle minimal representation.