36 #include <visp3/core/vpConfig.h> 38 #ifdef VISP_HAVE_CATCH2 39 #define CATCH_CONFIG_ENABLE_BENCHMARKING 40 #define CATCH_CONFIG_RUNNER 43 #include <visp3/core/vpMatrix.h> 45 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000) 46 #include <opencv2/core.hpp> 49 #ifdef VISP_HAVE_EIGEN3 50 #include <Eigen/Dense> 56 bool runBenchmark =
false;
57 bool runBenchmarkAll =
false;
59 double getRandomValues(
double min,
double max)
61 return (max - min) * ((double)rand() / (double)RAND_MAX) + min;
64 vpMatrix generateRandomMatrix(
unsigned int rows,
unsigned int cols,
double min=-1,
double max=1)
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);
77 vpColVector generateRandomVector(
unsigned int rows,
double min=-1,
double max=1)
81 for (
unsigned int i = 0; i < v.getRows(); i++) {
82 v[i] = getRandomValues(min, max);
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];
109 for (j = 0; j < BcolNum; j++) {
111 for (k = 0; k < BrowNum; k++) {
112 s += rowptri[k] * B[k][j];
127 for (
unsigned int i = 0; i < A.
getCols(); i++) {
129 for (
unsigned int j = 0; j < i; j++) {
130 double *ptr = A.
data;
132 for (
unsigned int k = 0; k < A.
getRows(); k++) {
133 s += (*(ptr + i)) * (*(ptr + j));
139 double *ptr = A.
data;
141 for (
unsigned int k = 0; k < A.
getRows(); k++) {
142 s += (*(ptr + i)) * (*(ptr + i));
158 for (
unsigned int i = 0; i < A.
getRows(); i++) {
159 for (
unsigned int j = i; j < A.
getRows(); j++) {
165 for (
unsigned int k = 0; k < A.
getCols(); k++)
166 ssum += *(pi++) * *(pj++);
188 for (
unsigned int j = 0; j < A.
getCols(); j++) {
190 for (
unsigned int i = 0; i < A.
getRows(); i++) {
191 w[i] += A[i][j] * vj;
204 for (
unsigned int i = 0; i < A.
getRows(); i++) {
205 for (
unsigned int j = 0; j < A.
getCols(); j++) {
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} };
221 for (
auto sz : sizes) {
222 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
223 vpMatrix B = generateRandomMatrix(sz.second, sz.first);
225 std::ostringstream oss;
228 BENCHMARK(oss.str().c_str()) {
229 C_true = dgemm_regular(A, B);
235 BENCHMARK(oss.str().c_str()) {
239 REQUIRE(equalMatrix(C, C_true));
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);
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];
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;
261 #ifdef VISP_HAVE_EIGEN3 262 Eigen::MatrixXd eigenA(sz.first, sz.second);
263 Eigen::MatrixXd eigenB(sz.second, sz.first);
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];
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;
284 const unsigned int rows = 47, cols = 63;
285 vpMatrix A = generateRandomMatrix(rows, cols);
286 vpMatrix B = generateRandomMatrix(cols, rows);
288 vpMatrix C_true = dgemm_regular(A, B);
290 REQUIRE(equalMatrix(C, C_true));
294 TEST_CASE(
"Benchmark matrix-rotation matrix multiplication",
"[benchmark]") {
295 if (runBenchmark || runBenchmarkAll) {
296 std::vector<std::pair<int, int>> sizes = { {3, 3} };
298 for (
auto sz : sizes) {
299 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
302 std::ostringstream oss;
305 BENCHMARK(oss.str().c_str()) {
306 AB_true = dgemm_regular(A, B);
312 BENCHMARK(oss.str().c_str()) {
316 REQUIRE(equalMatrix(AB, AB_true));
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);
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];
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];
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;
342 #ifdef VISP_HAVE_EIGEN3 343 Eigen::MatrixXd eigenA(sz.first, sz.second);
344 Eigen::MatrixXd eigenB(3, 3);
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];
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];
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;
369 const unsigned int rows = 3, cols = 3;
370 vpMatrix A = generateRandomMatrix(rows, cols);
373 vpMatrix AB_true = dgemm_regular(A, B);
375 REQUIRE(equalMatrix(AB, AB_true));
379 TEST_CASE(
"Benchmark matrix-homogeneous matrix multiplication",
"[benchmark]") {
380 if (runBenchmark || runBenchmarkAll) {
381 std::vector<std::pair<int, int>> sizes = { {4, 4} };
383 for (
auto sz : sizes) {
384 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
388 std::ostringstream oss;
391 BENCHMARK(oss.str().c_str()) {
392 AB_true = dgemm_regular(A, B);
398 BENCHMARK(oss.str().c_str()) {
402 REQUIRE(equalMatrix(AB, AB_true));
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);
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];
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];
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;
428 #ifdef VISP_HAVE_EIGEN3 429 Eigen::MatrixXd eigenA(sz.first, sz.second);
430 Eigen::MatrixXd eigenB(4, 4);
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];
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];
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;
455 const unsigned int rows = 4, cols = 4;
456 vpMatrix A = generateRandomMatrix(rows, cols);
460 vpMatrix AB_true = dgemm_regular(A, B);
463 REQUIRE(equalMatrix(AB, AB_true));
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} };
471 for (
auto sz : sizes) {
472 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
475 std::ostringstream oss;
478 BENCHMARK(oss.str().c_str()) {
479 C_true = dgemv_regular(A, B);
485 BENCHMARK(oss.str().c_str()) {
489 REQUIRE(equalMatrix(C, C_true));
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);
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];
500 matB.at<
double>(j, 0) = B[j];
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;
513 #ifdef VISP_HAVE_EIGEN3 514 Eigen::MatrixXd eigenA(sz.first, sz.second);
515 Eigen::MatrixXd eigenB(sz.second, 1);
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];
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;
538 const unsigned int rows = 47, cols = 63;
539 vpMatrix A = generateRandomMatrix(rows, cols);
544 REQUIRE(equalMatrix(C, C_true));
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} };
552 for (
auto sz : sizes) {
553 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
555 std::ostringstream oss;
556 oss <<
"(" << A.
getRows() <<
"x" << A.
getCols() <<
") - Naive code";
558 BENCHMARK(oss.str().c_str()) {
559 AtA_true = AtA_regular(A);
565 BENCHMARK(oss.str().c_str()) {
569 REQUIRE(equalMatrix(AtA, AtA_true));
571 if(runBenchmarkAll) {
572 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000) 573 cv::Mat matA(sz.first, sz.second, CV_64FC1);
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];
582 oss <<
"(" << matA.rows <<
"x" << matA.cols <<
") - OpenCV";
583 BENCHMARK(oss.str().c_str()) {
584 cv::Mat matAtA = matA.
t() * matA;
589 #ifdef VISP_HAVE_EIGEN3 590 Eigen::MatrixXd eigenA(sz.first, sz.second);
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];
599 oss <<
"(" << eigenA.rows() <<
"x" << eigenA.cols() <<
") - Eigen";
600 BENCHMARK(oss.str().c_str()) {
601 Eigen::MatrixXd eigenAtA = eigenA.
transpose() * eigenA;
610 const unsigned int rows = 47, cols = 63;
611 vpMatrix A = generateRandomMatrix(rows, cols);
615 REQUIRE(equalMatrix(AtA, AtA_true));
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} };
623 for (
auto sz : sizes) {
624 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
626 std::ostringstream oss;
627 oss <<
"(" << A.
getRows() <<
"x" << A.
getCols() <<
") - Naive code";
629 BENCHMARK(oss.str().c_str()) {
630 AAt_true = AAt_regular(A);
636 BENCHMARK(oss.str().c_str()) {
640 REQUIRE(equalMatrix(AAt, AAt_true));
642 if(runBenchmarkAll) {
643 #if (VISP_HAVE_OPENCV_VERSION >= 0x030000) 644 cv::Mat matA(sz.first, sz.second, CV_64FC1);
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];
653 oss <<
"(" << matA.rows <<
"x" << matA.cols <<
") - OpenCV";
654 BENCHMARK(oss.str().c_str()) {
655 cv::Mat matAAt = matA * matA.
t();
660 #ifdef VISP_HAVE_EIGEN3 661 Eigen::MatrixXd eigenA(sz.first, sz.second);
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];
670 oss <<
"(" << eigenA.rows() <<
"x" << eigenA.cols() <<
") - Eigen";
671 BENCHMARK(oss.str().c_str()) {
672 Eigen::MatrixXd eigenAAt = eigenA * eigenA.
transpose();
681 const unsigned int rows = 47, cols = 63;
682 vpMatrix A = generateRandomMatrix(rows, cols);
686 REQUIRE(equalMatrix(AAt, AAt_true));
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} };
694 for (
auto sz : sizes) {
695 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
698 std::ostringstream oss;
699 oss <<
"(" << A.
getRows() <<
"x" << A.
getCols() <<
")x(6x6) - Naive code";
701 BENCHMARK(oss.str().c_str()) {
702 AV_true = dgemm_regular(A, V);
707 oss <<
"(" << A.
getRows() <<
"x" << A.
getCols() <<
")x(6x6) - ViSP";
708 BENCHMARK(oss.str().c_str()) {
712 REQUIRE(equalMatrix(AV, AV_true));
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);
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];
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];
731 oss <<
"(" << matA.rows <<
"x" << matA.cols <<
")x(6x6) - OpenCV";
732 BENCHMARK(oss.str().c_str()) {
733 cv::Mat matAV = matA * matV;
738 #ifdef VISP_HAVE_EIGEN3 739 Eigen::MatrixXd eigenA(sz.first, sz.second);
740 Eigen::MatrixXd eigenV(6, 6);
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];
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];
754 oss <<
"(" << eigenA.rows() <<
"x" << eigenA.cols() <<
")x(6x6) - Eigen";
755 BENCHMARK(oss.str().c_str()) {
756 Eigen::MatrixXd eigenAV = eigenA * eigenV;
765 const unsigned int rows = 47, cols = 6;
766 vpMatrix A = generateRandomMatrix(rows, cols);
769 vpMatrix AV_true = dgemm_regular(A, V);
771 REQUIRE(equalMatrix(AV, AV_true));
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} };
779 for (
auto sz : sizes) {
780 vpMatrix A = generateRandomMatrix(sz.first, sz.second);
783 std::ostringstream oss;
784 oss <<
"(" << A.
getRows() <<
"x" << A.
getCols() <<
")x(6x6) - Naive code";
786 BENCHMARK(oss.str().c_str()) {
787 AV_true = dgemm_regular(A, V);
792 oss <<
"(" << A.
getRows() <<
"x" << A.
getCols() <<
")x(6x6) - ViSP";
793 BENCHMARK(oss.str().c_str()) {
797 REQUIRE(equalMatrix(AV, AV_true));
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);
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];
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];
816 oss <<
"(" << matA.rows <<
"x" << matA.cols <<
")x(6x6) - OpenCV";
817 BENCHMARK(oss.str().c_str()) {
818 cv::Mat matAV = matA * matV;
823 #ifdef VISP_HAVE_EIGEN3 824 Eigen::MatrixXd eigenA(sz.first, sz.second);
825 Eigen::MatrixXd eigenV(6, 6);
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];
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];
839 oss <<
"(" << eigenA.rows() <<
"x" << eigenA.cols() <<
")x(6x6) - Eigen";
840 BENCHMARK(oss.str().c_str()) {
841 Eigen::MatrixXd eigenAV = eigenA * eigenV;
850 const unsigned int rows = 47, cols = 6;
851 vpMatrix A = generateRandomMatrix(rows, cols);
854 vpMatrix AV_true = dgemm_regular(A, V);
856 REQUIRE(equalMatrix(AV, AV_true));
860 int main(
int argc,
char *argv[])
867 Catch::Session session;
870 std::cout <<
"Default matrix/vector min size to enable Blas/Lapack optimization: " 871 << lapackMinSize << std::endl;
873 using namespace Catch::clara;
874 auto cli = session.cli()
877 (
"run benchmark comparing naive code with ViSP implementation")
878 | Opt(runBenchmarkAll)
880 (
"run benchmark comparing naive code with ViSP, OpenCV, Eigen implementation")
881 | Opt(lapackMinSize,
"min size")
882 [
"--lapack-min-size"]
883 (
"matrix/vector min size to enable blas/lapack usage");
889 session.applyCommandLine(argc, argv);
892 std::cout <<
"Used matrix/vector min size to enable Blas/Lapack optimization: " 895 int numFailed = session.run();
Implementation of a matrix and operations on matrices.
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Implementation of an homogeneous matrix and operations on such kind of matrices.
static bool equal(double x, double y, double s=0.001)
error that can be emited by ViSP classes.
Type * data
Address of the first element of the data array.
unsigned int getCols() const
Implementation of a rotation matrix and operations on such kind of matrices.
static unsigned int getLapackMatrixMinSize()
static void mult2Matrices(const vpMatrix &A, const vpMatrix &B, vpMatrix &C)
vpMatrix transpose() const
unsigned int getRows() const
static void setLapackMatrixMinSize(unsigned int min_size)
vp_deprecated vpColVector rows(unsigned int first_row, unsigned int last_row) const
void resize(unsigned int i, bool flagNullify=true)
static double deg(double rad)
Implementation of column vector and the associated operations.
Class that consider the case of a translation vector.
Implementation of a rotation vector as axis-angle minimal representation.