48 #include <visp3/core/vpColVector.h>
49 #include <visp3/core/vpMatrix.h>
50 #include <visp3/core/vpTime.h>
51 #include <visp3/io/vpParseArgv.h>
54 #define GETOPTARGS "cdn:i:pf:R:C:vh"
64 void usage(
const char *name,
const char *badparam)
67 Test matrix inversions\n\
68 using LU, QR and Cholesky methods as well as Pseudo-inverse.\n\
69 Outputs a comparison of these methods.\n\
72 %s [-n <number of matrices>] [-f <plot filename>]\n\
73 [-R <number of rows>] [-C <number of columns>]\n\
74 [-i <number of iterations>] [-p] [-h]\n",
79 -n <number of matrices> \n\
80 Number of matrices inverted during each test loop.\n\
82 -i <number of iterations> \n\
83 Number of iterations of the test.\n\
85 -f <plot filename> \n\
86 Set output path for plot output.\n\
87 The plot logs the times of \n\
88 the different inversion methods: \n\
89 QR,LU,Cholesky and Pseudo-inverse.\n\
91 -R <number of rows>\n\
92 Number of rows of the automatically generated matrices \n\
95 -C <number of columns>\n\
96 Number of colums of the automatically generated matrices \n\
100 Plot into filename in the gnuplot format. \n\
101 If this option is used, tests results will be logged \n\
102 into a filename specified with -f.\n\
105 Print the help.\n\n");
108 fprintf(stderr,
"ERROR: \n");
109 fprintf(stderr,
"\nBad parameter [%s]\n", badparam);
120 bool getOptions(
int argc,
const char **argv,
unsigned int &nb_matrices,
unsigned int &nb_iterations,
121 bool &use_plot_file, std::string &plotfile,
unsigned int &nbrows,
unsigned int &nbcols,
bool &verbose)
129 usage(argv[0],
nullptr);
133 nb_matrices = (
unsigned int)atoi(optarg_);
136 nb_iterations = (
unsigned int)atoi(optarg_);
140 use_plot_file =
true;
143 use_plot_file =
true;
146 nbrows = (
unsigned int)atoi(optarg_);
149 nbcols = (
unsigned int)atoi(optarg_);
160 usage(argv[0], optarg_);
166 if ((c == 1) || (c == -1)) {
168 usage(argv[0],
nullptr);
169 std::cerr <<
"ERROR: " << std::endl;
170 std::cerr <<
" Bad argument " << optarg_ << std::endl << std::endl;
177 vpMatrix make_random_matrix(
unsigned int nbrows,
unsigned int nbcols)
182 for (
unsigned int i = 0; i < A.
getRows(); i++)
183 for (
unsigned int j = 0; j < A.
getCols(); j++)
184 A[i][j] = (
double)rand() / (double)RAND_MAX;
188 vpMatrix make_random_symmetric_positive_matrix(
unsigned int n)
195 for (
unsigned int i = 0; i < A.
getRows(); i++)
196 for (
unsigned int j = 0; j < A.
getCols(); j++)
197 A[i][j] = (
double)rand() / (double)RAND_MAX;
199 A = 0.5 * (A + A.
t());
204 vpMatrix make_random_triangular_matrix(
unsigned int nbrows)
209 for (
unsigned int i = 0; i < A.
getRows(); i++) {
210 for (
unsigned int j = i; j < A.
getCols(); j++) {
211 A[i][j] =
static_cast<double>(rand()) /
static_cast<double>(RAND_MAX);
221 void create_bench_random_matrix(
unsigned int nb_matrices,
unsigned int nb_rows,
unsigned int nb_cols,
bool verbose,
222 std::vector<vpMatrix> &bench)
225 std::cout <<
"Create a bench of " << nb_matrices <<
" " << nb_rows <<
" by " << nb_cols <<
" matrices" << std::endl;
227 for (
unsigned int i = 0; i < nb_matrices; i++) {
229 #if defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_OPENCV)
232 for (M = make_random_matrix(nb_rows, nb_cols); std::fabs(det = M.
AtA().
det()) < .01;
233 M = make_random_matrix(nb_rows, nb_cols)) {
235 std::cout <<
" Generated random matrix AtA=" << std::endl << M.
AtA() << std::endl;
236 std::cout <<
" Generated random matrix not invertible: det=" << det <<
". Retrying..." << std::endl;
240 M = make_random_matrix(nb_rows, nb_cols);
246 void create_bench_symmetric_positive_matrix(
unsigned int nb_matrices,
unsigned int n,
bool verbose,
247 std::vector<vpMatrix> &bench)
250 std::cout <<
"Create a bench of " << nb_matrices <<
" " << n <<
" by " << n <<
" symmetric positive matrices"
253 for (
unsigned int i = 0; i < nb_matrices; i++) {
255 #if defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_OPENCV)
258 for (M = make_random_symmetric_positive_matrix(n); std::fabs(det = M.
det()) < .01;
259 M = make_random_symmetric_positive_matrix(n)) {
261 std::cout <<
" Generated random symmetric positive matrix A=" << std::endl << M << std::endl;
262 std::cout <<
" Generated random symmetric positive matrix not "
264 << det <<
". Retrying..." << std::endl;
268 M = make_random_symmetric_positive_matrix(n);
274 void create_bench_random_triangular_matrix(
unsigned int nb_matrices,
unsigned int n,
bool verbose,
275 std::vector<vpMatrix> &bench)
278 std::cout <<
"Create a bench of " << nb_matrices <<
" " << n <<
" by " << n <<
" triangular matrices" << std::endl;
280 for (
unsigned int i = 0; i < nb_matrices; i++) {
282 #if defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_OPENCV)
285 for (M = make_random_triangular_matrix(n); std::fabs(det = M.
det()) < .01; M = make_random_triangular_matrix(n)) {
287 std::cout <<
" Generated random symmetric positive matrix A=" << std::endl << M << std::endl;
288 std::cout <<
" Generated random symmetric positive matrix not "
290 << det <<
". Retrying..." << std::endl;
294 M = make_random_triangular_matrix(n);
300 int test_inverse(
const std::vector<vpMatrix> &bench,
const std::vector<vpMatrix> &result)
302 double epsilon = 1e-10;
303 for (
unsigned int i = 0; i < bench.size(); i++) {
305 if (std::fabs(I.frobeniusNorm() - sqrt(
static_cast<double>(bench[0].AtA().getRows()))) > epsilon) {
306 std::cout <<
"Bad inverse[" << i <<
"]: " << I.frobeniusNorm() <<
" " << sqrt((
double)bench[0].AtA().getRows())
314 int test_inverse_lu_small(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
317 std::cout <<
"Test inverse by LU on small matrices" << std::endl;
320 std::cout <<
" Inverting " << bench[0].getRows() <<
"x" << bench[0].getCols() <<
" small matrix." << std::endl;
321 std::vector<vpMatrix> result(bench.size());
323 for (
unsigned int i = 0; i < bench.size(); i++) {
324 result[i] = bench[i].inverseByLU();
329 return test_inverse(bench, result);
332 #if defined(VISP_HAVE_EIGEN3)
333 int test_inverse_lu_eigen3(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
336 std::cout <<
"Test inverse by LU using Eigen3 3rd party" << std::endl;
339 std::cout <<
" Inverting " << bench[0].AtA().getRows() <<
"x" << bench[0].AtA().getCols()
340 <<
" matrix using LU decomposition (Eigen3)." << std::endl;
341 std::vector<vpMatrix> result(bench.size());
343 for (
unsigned int i = 0; i < bench.size(); i++) {
344 result[i] = bench[i].AtA().inverseByLUEigen3() * bench[i].transpose();
349 return test_inverse(bench, result);
353 #if defined(VISP_HAVE_LAPACK)
354 int test_inverse_lu_lapack(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
357 std::cout <<
"Test inverse by LU using Lapack 3rd party" << std::endl;
360 std::cout <<
" Inverting " << bench[0].AtA().getRows() <<
"x" << bench[0].AtA().getCols()
361 <<
" matrix using LU decomposition (Lapack)." << std::endl;
362 std::vector<vpMatrix> result(bench.size());
364 for (
unsigned int i = 0; i < bench.size(); i++) {
365 result[i] = bench[i].AtA().inverseByLULapack() * bench[i].transpose();
370 return test_inverse(bench, result);
373 int test_inverse_cholesky_lapack(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
376 std::cout <<
"Test inverse by Cholesky using Lapack 3rd party" << std::endl;
379 std::cout <<
" Inverting " << bench[0].AtA().getRows() <<
"x" << bench[0].AtA().getCols()
380 <<
" matrix using cholesky decomposition (Lapack)." << std::endl;
381 std::vector<vpMatrix> result(bench.size());
383 for (
unsigned int i = 0; i < bench.size(); i++) {
384 result[i] = bench[i].AtA().inverseByCholeskyLapack() * bench[i].transpose();
389 return test_inverse(bench, result);
392 int test_inverse_qr_lapack(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
395 std::cout <<
"Test inverse by QR using Lapack 3rd party" << std::endl;
398 std::cout <<
" Inverting " << bench[0].AtA().getRows() <<
"x" << bench[0].AtA().getCols()
399 <<
" matrix using QR decomposition (Lapack)" << std::endl;
400 std::vector<vpMatrix> result(bench.size());
402 for (
unsigned int i = 0; i < bench.size(); i++) {
403 result[i] = bench[i].AtA().inverseByQRLapack() * bench[i].transpose();
408 return test_inverse(bench, result);
412 #if defined(VISP_HAVE_OPENCV)
413 int test_inverse_lu_opencv(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
416 std::cout <<
"Test inverse by LU using OpenCV 3rd party" << std::endl;
419 std::cout <<
" Inverting " << bench[0].AtA().getRows() <<
"x" << bench[0].AtA().getCols()
420 <<
" matrix using LU decomposition (OpenCV)" << std::endl;
421 std::vector<vpMatrix> result(bench.size());
423 for (
unsigned int i = 0; i < bench.size(); i++) {
424 result[i] = bench[i].AtA().inverseByLUOpenCV() * bench[i].transpose();
429 return test_inverse(bench, result);
432 int test_inverse_cholesky_opencv(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
435 std::cout <<
"Test inverse by Cholesky using OpenCV 3rd party" << std::endl;
438 std::cout <<
" Inverting " << bench[0].AtA().getRows() <<
"x" << bench[0].AtA().getCols()
439 <<
" matrix using Cholesky decomposition (OpenCV)" << std::endl;
440 std::vector<vpMatrix> result(bench.size());
442 for (
unsigned int i = 0; i < bench.size(); i++) {
443 result[i] = bench[i].AtA().inverseByCholeskyOpenCV() * bench[i].transpose();
448 return test_inverse(bench, result);
452 #if defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_OPENCV)
454 int test_pseudo_inverse(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
457 std::cout <<
"Test pseudo inverse using either Lapack, Eigen3 or OpenCV 3rd party" << std::endl;
460 std::cout <<
" Pseudo inverting " << bench[0].AtA().getRows() <<
"x" << bench[0].AtA().getCols() <<
" matrix"
462 std::vector<vpMatrix> result(bench.size());
464 for (
unsigned int i = 0; i < bench.size(); i++) {
465 result[i] = bench[i].AtA().pseudoInverse() * bench[i].transpose();
470 return test_inverse(bench, result);
473 int test_inverse_triangular(
bool verbose,
const std::vector<vpMatrix> &bench,
double &time)
476 std::cout <<
"Test inverse triangular using Lapack" << std::endl;
479 std::cout <<
" Triangular inverse " << bench[0].getRows() <<
"x" << bench[0].getCols() <<
" matrix" << std::endl;
480 std::vector<vpMatrix> result(bench.size());
482 for (
unsigned int i = 0; i < bench.size(); i++) {
483 result[i] = bench[i].inverseTriangular(
true);
488 return test_inverse(bench, result);
492 void save_time(
const std::string &method,
bool verbose,
bool use_plot_file, std::ofstream &of,
double time)
496 if (verbose || !use_plot_file) {
497 std::cout << method << time << std::endl;
501 int main(
int argc,
const char *argv[])
504 unsigned int nb_matrices = 1000;
505 unsigned int nb_iterations = 10;
506 unsigned int nb_rows = 6;
507 unsigned int nb_cols = 6;
508 bool verbose =
false;
509 std::string plotfile(
"plot-inv.csv");
510 bool use_plot_file =
false;
514 if (getOptions(argc, argv, nb_matrices, nb_iterations, use_plot_file, plotfile, nb_rows, nb_cols, verbose) ==
520 of.open(plotfile.c_str());
524 #if defined(VISP_HAVE_LAPACK)
525 of <<
"\"LU Lapack\""
528 #if defined(VISP_HAVE_EIGEN3)
529 of <<
"\"LU Eigen3\""
532 #if defined(VISP_HAVE_OPENCV)
533 of <<
"\"LU OpenCV\""
537 #if defined(VISP_HAVE_LAPACK)
538 of <<
"\"Cholesky Lapack\""
542 #if defined(VISP_HAVE_OPENCV)
543 of <<
"\"Cholesky OpenCV\""
547 #if defined(VISP_HAVE_LAPACK)
548 of <<
"\"QR Lapack\""
552 #if defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_OPENCV)
553 of <<
"\"Pseudo inverse (Lapack, Eigen3 or OpenCV)\""
559 int ret = EXIT_SUCCESS;
560 for (
unsigned int iter = 0; iter < nb_iterations; iter++) {
561 std::vector<vpMatrix> bench_random_matrices_11;
562 create_bench_random_matrix(nb_matrices, 1, 1, verbose, bench_random_matrices_11);
563 std::vector<vpMatrix> bench_random_matrices_22;
564 create_bench_random_matrix(nb_matrices, 2, 2, verbose, bench_random_matrices_22);
565 std::vector<vpMatrix> bench_random_matrices_33;
566 create_bench_random_matrix(nb_matrices, 3, 3, verbose, bench_random_matrices_33);
567 #if defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_OPENCV)
568 std::vector<vpMatrix> bench_random_matrices;
569 create_bench_random_matrix(nb_matrices, nb_rows, nb_cols, verbose, bench_random_matrices);
570 std::vector<vpMatrix> bench_symmetric_positive_matrices;
571 create_bench_symmetric_positive_matrix(nb_matrices, nb_rows, verbose, bench_symmetric_positive_matrices);
572 std::vector<vpMatrix> bench_triangular_matrices;
573 create_bench_random_triangular_matrix(nb_matrices, nb_rows, verbose, bench_triangular_matrices);
582 ret += test_inverse_lu_small(verbose, bench_random_matrices_11, time);
583 save_time(
"Inverse by LU 1x1: ", verbose, use_plot_file, of, time);
585 ret += test_inverse_lu_small(verbose, bench_random_matrices_22, time);
586 save_time(
"Inverse by LU 2x2: ", verbose, use_plot_file, of, time);
588 ret += test_inverse_lu_small(verbose, bench_random_matrices_33, time);
589 save_time(
"Inverse by LU 3x3: ", verbose, use_plot_file, of, time);
592 #if defined(VISP_HAVE_LAPACK)
593 ret += test_inverse_lu_lapack(verbose, bench_random_matrices, time);
594 save_time(
"Inverse by LU (Lapack): ", verbose, use_plot_file, of, time);
597 #if defined(VISP_HAVE_EIGEN3)
598 ret += test_inverse_lu_eigen3(verbose, bench_random_matrices, time);
599 save_time(
"Inverse by LU (Eigen3): ", verbose, use_plot_file, of, time);
602 #if defined(VISP_HAVE_OPENCV)
603 ret += test_inverse_lu_opencv(verbose, bench_random_matrices, time);
604 save_time(
"Inverse by LU (OpenCV): ", verbose, use_plot_file, of, time);
608 #if defined(VISP_HAVE_LAPACK)
609 ret += test_inverse_cholesky_lapack(verbose, bench_symmetric_positive_matrices, time);
610 save_time(
"Inverse by Cholesly (Lapack): ", verbose, use_plot_file, of, time);
613 #if defined(VISP_HAVE_OPENCV)
614 ret += test_inverse_cholesky_opencv(verbose, bench_symmetric_positive_matrices, time);
615 save_time(
"Inverse by Cholesky (OpenCV): ", verbose, use_plot_file, of, time);
619 #if defined(VISP_HAVE_LAPACK)
620 ret += test_inverse_qr_lapack(verbose, bench_random_matrices, time);
621 save_time(
"Inverse by QR (Lapack): ", verbose, use_plot_file, of, time);
625 #if defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_OPENCV)
626 ret += test_pseudo_inverse(verbose, bench_random_matrices, time);
627 save_time(
"Pseudo inverse (Lapack, Eigen3, OpenCV): ", verbose, use_plot_file, of, time);
631 #if defined(VISP_HAVE_LAPACK)
632 ret += test_inverse_triangular(verbose, bench_triangular_matrices, time);
633 save_time(
"Triangular inverse (Lapack): ", verbose, use_plot_file, of, time);
641 std::cout <<
"Result saved in " << plotfile << std::endl;
644 if (ret == EXIT_SUCCESS) {
645 std::cout <<
"Test succeed" << std::endl;
647 std::cout <<
"Test failed" << std::endl;
unsigned int getCols() const
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
unsigned int getRows() const
error that can be emitted by ViSP classes.
const std::string & getStringMessage() const
Implementation of a matrix and operations on matrices.
double det(vpDetMethod method=LU_DECOMPOSITION) const
static bool parse(int *argcPtr, const char **argv, vpArgvInfo *argTable, int flags)
VISP_EXPORT double measureTimeMs()