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