Visual Servoing Platform version 3.6.0
Loading...
Searching...
No Matches
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
53namespace
54{
55
56bool runBenchmark = false;
57bool runBenchmarkAll = false;
58
59double getRandomValues(double min, double max) { return (max - min) * ((double)rand() / (double)RAND_MAX) + min; }
60
61vpMatrix 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
74vpColVector 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
86vpMatrix 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
119vpMatrix 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()
149vpMatrix 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
174vpColVector 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
195bool 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
214TEST_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
298TEST_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
391TEST_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
487TEST_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
575TEST_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
652TEST_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
730TEST_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
820TEST_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
910int 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
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
955int main() { return EXIT_SUCCESS; }
956#endif
unsigned int getCols() const
Definition vpArray2D.h:280
Type * data
Address of the first element of the data array.
Definition vpArray2D.h:144
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition vpArray2D.h:305
unsigned int getRows() const
Definition vpArray2D.h:290
Implementation of column vector and the associated operations.
void resize(unsigned int i, bool flagNullify=true)
error that can be emitted by ViSP classes.
Definition vpException.h:59
@ dimensionError
Bad dimension.
Definition vpException.h:83
Implementation of an homogeneous matrix and operations on such kind of matrices.
static bool equal(double x, double y, double threshold=0.001)
Definition vpMath.h:369
static double deg(double rad)
Definition vpMath.h:106
Implementation of a matrix and operations on matrices.
Definition vpMatrix.h:152
static void setLapackMatrixMinSize(unsigned int min_size)
Definition vpMatrix.h:259
static unsigned int getLapackMatrixMinSize()
Definition vpMatrix.h:245
vpMatrix AtA() const
Definition vpMatrix.cpp:625
vpMatrix AAt() const
Definition vpMatrix.cpp:501
static void mult2Matrices(const vpMatrix &A, const vpMatrix &B, vpMatrix &C)
Implementation of a rotation matrix and operations on such kind of matrices.
Implementation of a rotation vector as axis-angle minimal representation.
Class that consider the case of a translation vector.