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