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