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