Visual Servoing Platform version 3.7.0
Loading...
Searching...
No Matches
perfMatrixMultiplication.cpp
1/*
2 * ViSP, open source Visual Servoing Platform software.
3 * Copyright (C) 2005 - 2025 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
37
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 defined(VISP_HAVE_OPENCV) && (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
55using namespace VISP_NAMESPACE_NAME;
56#endif
57
58namespace
59{
60
61bool runBenchmark = false;
62bool runBenchmarkAll = false;
63
64double getRandomValues(double min, double max) { return (max - min) * (static_cast<double>(rand()) / static_cast<double>(RAND_MAX)) + min; }
65
66vpMatrix 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
79vpColVector 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
91vpMatrix 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
124vpMatrix 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()
154vpMatrix 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
179vpColVector 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
200bool 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
219TEST_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 defined(VISP_HAVE_OPENCV) && (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
303TEST_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 defined(VISP_HAVE_OPENCV) && (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
396TEST_CASE("Benchmark rotation matrix-matrix multiplication", "[benchmark]")
397{
398 if (runBenchmark || runBenchmarkAll) {
399 std::vector<std::pair<int, int> > sizes = { {3, 3} };
400
401 for (auto sz : sizes) {
402 vpRotationMatrix A(vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)),
403 vpMath::deg(getRandomValues(0, 360)));
404 vpMatrix B = generateRandomMatrix(sz.first, sz.second);
405
406 std::ostringstream oss;
407 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - Naive code";
408 vpMatrix AB, AB_true;
409 BENCHMARK(oss.str().c_str())
410 {
411 AB_true = dgemm_regular(static_cast<vpMatrix>(A), B);
412 return AB_true;
413 };
414
415 oss.str("");
416 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - ViSP";
417 BENCHMARK(oss.str().c_str())
418 {
419 AB = A * B;
420 return AB;
421 };
422 REQUIRE(equalMatrix(AB, AB_true));
423
424 if (runBenchmarkAll) {
425#if defined(VISP_HAVE_OPENCV) && (VISP_HAVE_OPENCV_VERSION >= 0x030000)
426 cv::Mat matA(3, 3, CV_64FC1);
427 cv::Mat matB(sz.first, sz.second, CV_64FC1);
428
429 for (unsigned int i = 0; i < A.getRows(); i++) {
430 for (unsigned int j = 0; j < A.getCols(); j++) {
431 matA.at<double>(i, j) = A[i][j];
432 }
433 }
434 for (unsigned int i = 0; i < B.getRows(); i++) {
435 for (unsigned int j = 0; j < B.getCols(); j++) {
436 matB.at<double>(j, i) = B[j][i];
437 }
438 }
439
440 oss.str("");
441 oss << "(" << matA.rows << "x" << matA.cols << ")x(" << matB.rows << "x" << matB.cols << ") - OpenCV";
442 BENCHMARK(oss.str().c_str())
443 {
444 cv::Mat matC = matA * matB;
445 return matC;
446 };
447#endif
448
449#ifdef VISP_HAVE_EIGEN3
450 Eigen::MatrixXd eigenA(3, 3);
451 Eigen::MatrixXd eigenB(sz.first, sz.second);
452
453 for (unsigned int i = 0; i < A.getRows(); i++) {
454 for (unsigned int j = 0; j < A.getCols(); j++) {
455 eigenA(i, j) = A[i][j];
456 }
457 }
458 for (unsigned int i = 0; i < B.getRows(); i++) {
459 for (unsigned int j = 0; j < B.getCols(); j++) {
460 eigenB(j, i) = B[j][i];
461 }
462 }
463
464 oss.str("");
465 oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(" << eigenB.rows() << "x" << eigenB.cols()
466 << ") - Eigen";
467 BENCHMARK(oss.str().c_str())
468 {
469 Eigen::MatrixXd eigenC = eigenA * eigenB;
470 return eigenC;
471 };
472#endif
473 }
474 }
475 }
476
477 {
478 const unsigned int rows = 3, cols = 3;
479 vpRotationMatrix A(vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)),
480 vpMath::deg(getRandomValues(0, 360)));
481 vpMatrix B = generateRandomMatrix(rows, cols);
482
483 vpMatrix AB_true = dgemm_regular(static_cast<vpMatrix>(A), B);
484 vpMatrix AB = A * B;
485 REQUIRE(equalMatrix(AB, AB_true));
486 }
487}
488
489TEST_CASE("Benchmark matrix-homogeneous matrix multiplication", "[benchmark]")
490{
491 if (runBenchmark || runBenchmarkAll) {
492 std::vector<std::pair<int, int> > sizes = { {4, 4} };
493
494 for (auto sz : sizes) {
495 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
496 vpHomogeneousMatrix B(getRandomValues(0, 1), getRandomValues(0, 1), getRandomValues(0, 1),
497 vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)),
498 vpMath::deg(getRandomValues(0, 360)));
499
500 std::ostringstream oss;
501 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - Naive code";
502 vpMatrix AB, AB_true;
503 BENCHMARK(oss.str().c_str())
504 {
505 AB_true = dgemm_regular(A, static_cast<vpMatrix>(B));
506 return AB_true;
507 };
508
509 oss.str("");
510 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - ViSP";
511 BENCHMARK(oss.str().c_str())
512 {
513 AB = A * B;
514 return AB;
515 };
516 REQUIRE(equalMatrix(AB, AB_true));
517
518 if (runBenchmarkAll) {
519#if defined(VISP_HAVE_OPENCV) && (VISP_HAVE_OPENCV_VERSION >= 0x030000)
520 cv::Mat matA(sz.first, sz.second, CV_64FC1);
521 cv::Mat matB(4, 4, CV_64FC1);
522
523 for (unsigned int i = 0; i < A.getRows(); i++) {
524 for (unsigned int j = 0; j < A.getCols(); j++) {
525 matA.at<double>(i, j) = A[i][j];
526 }
527 }
528 for (unsigned int i = 0; i < B.getRows(); i++) {
529 for (unsigned int j = 0; j < B.getCols(); j++) {
530 matB.at<double>(j, i) = B[j][i];
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(4, 4);
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 }
551 }
552 for (unsigned int i = 0; i < B.getRows(); i++) {
553 for (unsigned int j = 0; j < B.getCols(); j++) {
554 eigenB(j, i) = B[j][i];
555 }
556 }
557
558 oss.str("");
559 oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(" << eigenB.rows() << "x" << eigenB.cols()
560 << ") - Eigen";
561 BENCHMARK(oss.str().c_str())
562 {
563 Eigen::MatrixXd eigenC = eigenA * eigenB;
564 return eigenC;
565 };
566#endif
567 }
568 }
569 }
570
571 {
572 const unsigned int rows = 4, cols = 4;
573 vpMatrix A = generateRandomMatrix(rows, cols);
574 vpHomogeneousMatrix B(getRandomValues(0, 1), getRandomValues(0, 1), getRandomValues(0, 1),
575 vpMath::deg(getRandomValues(0, 360)), vpMath::deg(getRandomValues(0, 360)),
576 vpMath::deg(getRandomValues(0, 360)));
577
578 vpMatrix AB_true = dgemm_regular(A, static_cast<vpMatrix>(B));
579 vpMatrix AB;
580 vpMatrix::mult2Matrices(A, static_cast<vpMatrix>(B), AB);
581 REQUIRE(equalMatrix(AB, AB_true));
582 }
583}
584
585TEST_CASE("Benchmark matrix-vector multiplication", "[benchmark]")
586{
587 if (runBenchmark || runBenchmarkAll) {
588 std::vector<std::pair<int, int> > sizes = { {3, 3}, {6, 6}, {8, 8}, {10, 10}, {20, 20}, {6, 200},
589 {200, 6}, {207, 119}, {83, 201}, {600, 400}, {400, 600} };
590
591 for (auto sz : sizes) {
592 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
593 vpColVector B = generateRandomVector(sz.second);
594
595 std::ostringstream oss;
596 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - Naive code";
597 vpColVector C, C_true;
598 BENCHMARK(oss.str().c_str())
599 {
600 C_true = dgemv_regular(A, static_cast<vpColVector>(B));
601 return C_true;
602 };
603
604 oss.str("");
605 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(" << B.getRows() << "x" << B.getCols() << ") - ViSP";
606 BENCHMARK(oss.str().c_str())
607 {
608 C = A * B;
609 return C;
610 };
611 REQUIRE(equalMatrix(static_cast<vpMatrix>(C), static_cast<vpMatrix>(C_true)));
612
613 if (runBenchmarkAll) {
614#if defined(VISP_HAVE_OPENCV) && (VISP_HAVE_OPENCV_VERSION >= 0x030000)
615 cv::Mat matA(sz.first, sz.second, CV_64FC1);
616 cv::Mat matB(sz.second, 1, CV_64FC1);
617
618 for (unsigned int i = 0; i < A.getRows(); i++) {
619 for (unsigned int j = 0; j < A.getCols(); j++) {
620 matA.at<double>(i, j) = A[i][j];
621 if (i == 0) {
622 matB.at<double>(j, 0) = B[j];
623 }
624 }
625 }
626
627 oss.str("");
628 oss << "(" << matA.rows << "x" << matA.cols << ")x(" << matB.rows << "x" << matB.cols << ") - OpenCV";
629 BENCHMARK(oss.str().c_str())
630 {
631 cv::Mat matC = matA * matB;
632 return matC;
633 };
634#endif
635
636#ifdef VISP_HAVE_EIGEN3
637 Eigen::MatrixXd eigenA(sz.first, sz.second);
638 Eigen::MatrixXd eigenB(sz.second, 1);
639
640 for (unsigned int i = 0; i < A.getRows(); i++) {
641 for (unsigned int j = 0; j < A.getCols(); j++) {
642 eigenA(i, j) = A[i][j];
643 if (i == 0) {
644 eigenB(j, 0) = B[j];
645 }
646 }
647 }
648
649 oss.str("");
650 oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(" << eigenB.rows() << "x" << eigenB.cols()
651 << ") - Eigen";
652 BENCHMARK(oss.str().c_str())
653 {
654 Eigen::MatrixXd eigenC = eigenA * eigenB;
655 return eigenC;
656 };
657#endif
658 }
659 }
660 }
661
662 {
663 const unsigned int rows = 47, cols = 63;
664 vpMatrix A = generateRandomMatrix(rows, cols);
665 vpColVector B = generateRandomVector(cols);
666
667 vpColVector C_true = dgemv_regular(A, B);
668 vpColVector C = A * B;
669 REQUIRE(equalMatrix(static_cast<vpMatrix>(C), static_cast<vpMatrix>(C_true)));
670 }
671}
672
673TEST_CASE("Benchmark AtA", "[benchmark]")
674{
675 if (runBenchmark || runBenchmarkAll) {
676 std::vector<std::pair<int, int> > sizes = { {3, 3}, {6, 6}, {8, 8}, {10, 10}, {20, 20}, {6, 200},
677 {200, 6}, {207, 119}, {83, 201}, {600, 400}, {400, 600} };
678
679 for (auto sz : sizes) {
680 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
681
682 std::ostringstream oss;
683 oss << "(" << A.getRows() << "x" << A.getCols() << ") - Naive code";
684 vpMatrix AtA, AtA_true;
685 BENCHMARK(oss.str().c_str())
686 {
687 AtA_true = AtA_regular(A);
688 return AtA_true;
689 };
690
691 oss.str("");
692 oss << "(" << A.getRows() << "x" << A.getCols() << ") - ViSP";
693 BENCHMARK(oss.str().c_str())
694 {
695 AtA = A.AtA();
696 return AtA;
697 };
698 REQUIRE(equalMatrix(AtA, AtA_true));
699
700 if (runBenchmarkAll) {
701#if defined(VISP_HAVE_OPENCV) && (VISP_HAVE_OPENCV_VERSION >= 0x030000)
702 cv::Mat matA(sz.first, sz.second, CV_64FC1);
703
704 for (unsigned int i = 0; i < A.getRows(); i++) {
705 for (unsigned int j = 0; j < A.getCols(); j++) {
706 matA.at<double>(i, j) = A[i][j];
707 }
708 }
709
710 oss.str("");
711 oss << "(" << matA.rows << "x" << matA.cols << ") - OpenCV";
712 BENCHMARK(oss.str().c_str())
713 {
714 cv::Mat matAtA = matA.t() * matA;
715 return matAtA;
716 };
717#endif
718
719#ifdef VISP_HAVE_EIGEN3
720 Eigen::MatrixXd eigenA(sz.first, sz.second);
721
722 for (unsigned int i = 0; i < A.getRows(); i++) {
723 for (unsigned int j = 0; j < A.getCols(); j++) {
724 eigenA(i, j) = A[i][j];
725 }
726 }
727
728 oss.str("");
729 oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ") - Eigen";
730 BENCHMARK(oss.str().c_str())
731 {
732 Eigen::MatrixXd eigenAtA = eigenA.transpose() * eigenA;
733 return eigenAtA;
734 };
735#endif
736 }
737 }
738 }
739
740 {
741 const unsigned int rows = 47, cols = 63;
742 vpMatrix A = generateRandomMatrix(rows, cols);
743
744 vpMatrix AtA_true = AtA_regular(A);
745 vpMatrix AtA = A.AtA();
746 REQUIRE(equalMatrix(AtA, AtA_true));
747 }
748}
749
750TEST_CASE("Benchmark AAt", "[benchmark]")
751{
752 if (runBenchmark || runBenchmarkAll) {
753 std::vector<std::pair<int, int> > sizes = {
754 {3, 3}, {6, 6}, {8, 8}, {10, 10},
755 {20, 20}, {6, 200}, {200, 6} }; //, {207, 119}, {83, 201}, {600, 400}, {400, 600} };
756
757 for (auto sz : sizes) {
758 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
759
760 std::ostringstream oss;
761 oss << "(" << A.getRows() << "x" << A.getCols() << ") - Naive code";
762 vpMatrix AAt_true, AAt;
763 BENCHMARK(oss.str().c_str())
764 {
765 AAt_true = AAt_regular(A);
766 return AAt_true;
767 };
768
769 oss.str("");
770 oss << "(" << A.getRows() << "x" << A.getCols() << ") - ViSP";
771 BENCHMARK(oss.str().c_str())
772 {
773 AAt = A.AAt();
774 return AAt;
775 };
776 REQUIRE(equalMatrix(AAt, AAt_true));
777
778 if (runBenchmarkAll) {
779#if defined(VISP_HAVE_OPENCV) && (VISP_HAVE_OPENCV_VERSION >= 0x030000)
780 cv::Mat matA(sz.first, sz.second, CV_64FC1);
781
782 for (unsigned int i = 0; i < A.getRows(); i++) {
783 for (unsigned int j = 0; j < A.getCols(); j++) {
784 matA.at<double>(i, j) = A[i][j];
785 }
786 }
787
788 oss.str("");
789 oss << "(" << matA.rows << "x" << matA.cols << ") - OpenCV";
790 BENCHMARK(oss.str().c_str())
791 {
792 cv::Mat matAAt = matA * matA.t();
793 return matAAt;
794 };
795#endif
796
797#ifdef VISP_HAVE_EIGEN3
798 Eigen::MatrixXd eigenA(sz.first, sz.second);
799
800 for (unsigned int i = 0; i < A.getRows(); i++) {
801 for (unsigned int j = 0; j < A.getCols(); j++) {
802 eigenA(i, j) = A[i][j];
803 }
804 }
805
806 oss.str("");
807 oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ") - Eigen";
808 BENCHMARK(oss.str().c_str())
809 {
810 Eigen::MatrixXd eigenAAt = eigenA * eigenA.transpose();
811 return eigenAAt;
812 };
813#endif
814 }
815 }
816 }
817
818 {
819 const unsigned int rows = 47, cols = 63;
820 vpMatrix A = generateRandomMatrix(rows, cols);
821
822 vpMatrix AAt_true = AAt_regular(A);
823 vpMatrix AAt = A.AAt();
824 REQUIRE(equalMatrix(AAt, AAt_true));
825 }
826}
827
828TEST_CASE("Benchmark matrix-velocity twist multiplication", "[benchmark]")
829{
830 if (runBenchmark || runBenchmarkAll) {
831 std::vector<std::pair<int, int> > sizes = { {6, 6}, {20, 6}, {207, 6}, {600, 6}, {1201, 6} };
832
833 for (auto sz : sizes) {
834 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
835 vpVelocityTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
836
837 std::ostringstream oss;
838 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - Naive code";
839 vpMatrix AV, AV_true;
840 BENCHMARK(oss.str().c_str())
841 {
842 AV_true = dgemm_regular(A, static_cast<vpMatrix>(V));
843 return AV_true;
844 };
845
846 oss.str("");
847 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - ViSP";
848 BENCHMARK(oss.str().c_str())
849 {
850 AV = A * V;
851 return AV;
852 };
853 REQUIRE(equalMatrix(AV, AV_true));
854
855 if (runBenchmarkAll) {
856#if defined(VISP_HAVE_OPENCV) && (VISP_HAVE_OPENCV_VERSION >= 0x030000)
857 cv::Mat matA(sz.first, sz.second, CV_64FC1);
858 cv::Mat matV(6, 6, CV_64FC1);
859
860 for (unsigned int i = 0; i < A.getRows(); i++) {
861 for (unsigned int j = 0; j < A.getCols(); j++) {
862 matA.at<double>(i, j) = A[i][j];
863 }
864 }
865 for (unsigned int i = 0; i < V.getRows(); i++) {
866 for (unsigned int j = 0; j < V.getCols(); j++) {
867 matV.at<double>(i, j) = V[i][j];
868 }
869 }
870
871 oss.str("");
872 oss << "(" << matA.rows << "x" << matA.cols << ")x(6x6) - OpenCV";
873 BENCHMARK(oss.str().c_str())
874 {
875 cv::Mat matAV = matA * matV;
876 return matAV;
877 };
878#endif
879
880#ifdef VISP_HAVE_EIGEN3
881 Eigen::MatrixXd eigenA(sz.first, sz.second);
882 Eigen::MatrixXd eigenV(6, 6);
883
884 for (unsigned int i = 0; i < A.getRows(); i++) {
885 for (unsigned int j = 0; j < A.getCols(); j++) {
886 eigenA(i, j) = A[i][j];
887 }
888 }
889 for (unsigned int i = 0; i < V.getRows(); i++) {
890 for (unsigned int j = 0; j < V.getCols(); j++) {
891 eigenV(i, j) = V[i][j];
892 }
893 }
894
895 oss.str("");
896 oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(6x6) - Eigen";
897 BENCHMARK(oss.str().c_str())
898 {
899 Eigen::MatrixXd eigenAV = eigenA * eigenV;
900 return eigenAV;
901 };
902#endif
903 }
904 }
905 }
906
907 {
908 const unsigned int rows = 47, cols = 6;
909 vpMatrix A = generateRandomMatrix(rows, cols);
910 vpVelocityTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
911
912 vpMatrix AV_true = dgemm_regular(A, static_cast<vpMatrix>(V));
913 vpMatrix AV = A * V;
914 REQUIRE(equalMatrix(AV, AV_true));
915 }
916}
917
918TEST_CASE("Benchmark matrix-force twist multiplication", "[benchmark]")
919{
920 if (runBenchmark || runBenchmarkAll) {
921 std::vector<std::pair<int, int> > sizes = { {6, 6}, {20, 6}, {207, 6}, {600, 6}, {1201, 6} };
922
923 for (auto sz : sizes) {
924 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
925 vpForceTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
926
927 std::ostringstream oss;
928 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - Naive code";
929 vpMatrix AV, AV_true;
930 BENCHMARK(oss.str().c_str())
931 {
932 AV_true = dgemm_regular(A, static_cast<vpMatrix>(V));
933 return AV_true;
934 };
935
936 oss.str("");
937 oss << "(" << A.getRows() << "x" << A.getCols() << ")x(6x6) - ViSP";
938 BENCHMARK(oss.str().c_str())
939 {
940 AV = A * V;
941 return AV;
942 };
943 REQUIRE(equalMatrix(AV, AV_true));
944
945 if (runBenchmarkAll) {
946#if defined(VISP_HAVE_OPENCV) && (VISP_HAVE_OPENCV_VERSION >= 0x030000)
947 cv::Mat matA(sz.first, sz.second, CV_64FC1);
948 cv::Mat matV(6, 6, CV_64FC1);
949
950 for (unsigned int i = 0; i < A.getRows(); i++) {
951 for (unsigned int j = 0; j < A.getCols(); j++) {
952 matA.at<double>(i, j) = A[i][j];
953 }
954 }
955 for (unsigned int i = 0; i < V.getRows(); i++) {
956 for (unsigned int j = 0; j < V.getCols(); j++) {
957 matV.at<double>(i, j) = V[i][j];
958 }
959 }
960
961 oss.str("");
962 oss << "(" << matA.rows << "x" << matA.cols << ")x(6x6) - OpenCV";
963 BENCHMARK(oss.str().c_str())
964 {
965 cv::Mat matAV = matA * matV;
966 return matAV;
967 };
968#endif
969
970#ifdef VISP_HAVE_EIGEN3
971 Eigen::MatrixXd eigenA(sz.first, sz.second);
972 Eigen::MatrixXd eigenV(6, 6);
973
974 for (unsigned int i = 0; i < A.getRows(); i++) {
975 for (unsigned int j = 0; j < A.getCols(); j++) {
976 eigenA(i, j) = A[i][j];
977 }
978 }
979 for (unsigned int i = 0; i < V.getRows(); i++) {
980 for (unsigned int j = 0; j < V.getCols(); j++) {
981 eigenV(i, j) = V[i][j];
982 }
983 }
984
985 oss.str("");
986 oss << "(" << eigenA.rows() << "x" << eigenA.cols() << ")x(6x6) - Eigen";
987 BENCHMARK(oss.str().c_str())
988 {
989 Eigen::MatrixXd eigenAV = eigenA * eigenV;
990 return eigenAV;
991 };
992#endif
993 }
994 }
995 }
996
997 {
998 const unsigned int rows = 47, cols = 6;
999 vpMatrix A = generateRandomMatrix(rows, cols);
1000 vpForceTwistMatrix V(vpTranslationVector(0.1, -0.4, 1.5), vpThetaUVector(0.4, -0.1, 0.7));
1001
1002 vpMatrix AV_true = dgemm_regular(A, static_cast<vpMatrix>(V));
1003 vpMatrix AV = A * V;
1004 REQUIRE(equalMatrix(AV, AV_true));
1005 }
1006}
1007
1008int main(int argc, char *argv[])
1009{
1010 // Set random seed explicitly to avoid confusion
1011 // See: https://en.cppreference.com/w/cpp/numeric/random/srand
1012 // If rand() is used before any calls to srand(), rand() behaves as if it was seeded with srand(1).
1013 srand(1);
1014
1015 Catch::Session session;
1016 unsigned int lapackMinSize = vpMatrix::getLapackMatrixMinSize();
1017
1018 std::cout << "Default matrix/vector min size to enable Blas/Lapack optimization: " << lapackMinSize << std::endl;
1019
1020 auto cli = session.cli()
1021 | Catch::Clara::Opt(runBenchmark)["--benchmark"]("run benchmark comparing naive code with ViSP implementation")
1022 | Catch::Clara::Opt(runBenchmarkAll)["--benchmark-all"]("run benchmark comparing naive code with ViSP, OpenCV, Eigen implementation")
1023 | Catch::Clara::Opt(lapackMinSize, "min size")["--lapack-min-size"]("matrix/vector min size to enable blas/lapack usage");
1024
1025 session.cli(cli);
1026 session.applyCommandLine(argc, argv);
1027
1028 vpMatrix::setLapackMatrixMinSize(lapackMinSize);
1029 std::cout << "Used matrix/vector min size to enable Blas/Lapack optimization: " << vpMatrix::getLapackMatrixMinSize()
1030 << std::endl;
1031
1032 int numFailed = session.run();
1033
1034 return numFailed;
1035}
1036#else
1037#include <iostream>
1038
1039int main() { return EXIT_SUCCESS; }
1040#endif
vpArray2D::getCols
unsigned int getCols() const
Definition vpArray2D.h:423
vpArray2D::data
Type * data
Address of the first element of the data array.
Definition vpArray2D.h:149
vpArray2D::resize
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition vpArray2D.h:448
vpArray2D::getRows
unsigned int getRows() const
Definition vpArray2D.h:433
vpColVector
Implementation of column vector and the associated operations.
Definition vpColVector.h:191
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:470
vpMath::deg
static double deg(double rad)
Definition vpMath.h:119
vpMatrix
Implementation of a matrix and operations on matrices.
Definition vpMatrix.h:175
vpMatrix::setLapackMatrixMinSize
static void setLapackMatrixMinSize(unsigned int min_size)
Definition vpMatrix.h:290
vpMatrix::getLapackMatrixMinSize
static unsigned int getLapackMatrixMinSize()
Definition vpMatrix.h:276
vpMatrix::AtA
vpMatrix AtA() const
Definition vpMatrix_operations.cpp:785
vpMatrix::AAt
vpMatrix AAt() const
Definition vpMatrix_operations.cpp:644
vpMatrix::mult2Matrices
static void mult2Matrices(const vpMatrix &A, const vpMatrix &B, vpMatrix &C)
Definition vpMatrix_operations.cpp:158
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