Visual Servoing Platform version 3.5.0
testMatrixDeterminant.cpp
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 * Test various determinant computation methods.
33 *
34 * Authors:
35 * Fabien Spindler
36 *
37 *****************************************************************************/
38
44#include <visp3/core/vpMatrix.h>
45#include <visp3/core/vpTime.h>
46#include <visp3/io/vpParseArgv.h>
47
48// List of allowed command line options
49#define GETOPTARGS "cdn:i:pf:R:C:vh"
50
59void usage(const char *name, const char *badparam)
60{
61 fprintf(stdout, "\n\
62Test matrix inversions\n\
63using LU, QR and Cholesky methods as well as Pseudo-inverse.\n\
64Outputs a comparison of these methods.\n\
65\n\
66SYNOPSIS\n\
67 %s [-n <number of matrices>] [-f <plot filename>]\n\
68 [-R <number of rows>] [-C <number of columns>]\n\
69 [-i <number of iterations>] [-p] [-h]\n", name);
70
71 fprintf(stdout, "\n\
72OPTIONS: Default\n\
73 -n <number of matrices> \n\
74 Number of matrices inverted during each test loop.\n\
75\n\
76 -i <number of iterations> \n\
77 Number of iterations of the test.\n\
78\n\
79 -f <plot filename> \n\
80 Set output path for plot output.\n\
81 The plot logs the times of \n\
82 the different inversion methods: \n\
83 QR,LU,Cholesky and Pseudo-inverse.\n\
84\n\
85 -R <number of rows>\n\
86 Number of rows of the automatically generated matrices \n\
87 we test on.\n\
88\n\
89 -C <number of columns>\n\
90 Number of colums of the automatically generated matrices \n\
91 we test on.\n\
92\n\
93 -p \n\
94 Plot into filename in the gnuplot format. \n\
95 If this option is used, tests results will be logged \n\
96 into a filename specified with -f.\n\
97\n\
98 -h\n\
99 Print the help.\n\n");
100
101 if (badparam) {
102 fprintf(stderr, "ERROR: \n");
103 fprintf(stderr, "\nBad parameter [%s]\n", badparam);
104 }
105}
106
114bool getOptions(int argc, const char **argv, unsigned int &nb_matrices, unsigned int &nb_iterations,
115 bool &use_plot_file, std::string &plotfile, unsigned int &nbrows, unsigned int &nbcols, bool &verbose)
116{
117 const char *optarg_;
118 int c;
119 while ((c = vpParseArgv::parse(argc, argv, GETOPTARGS, &optarg_)) > 1) {
120
121 switch (c) {
122 case 'h':
123 usage(argv[0], NULL);
124 return false;
125 break;
126 case 'n':
127 nb_matrices = (unsigned int)atoi(optarg_);
128 break;
129 case 'i':
130 nb_iterations = (unsigned int)atoi(optarg_);
131 break;
132 case 'f':
133 plotfile = optarg_;
134 use_plot_file = true;
135 break;
136 case 'p':
137 use_plot_file = true;
138 break;
139 case 'R':
140 nbrows = (unsigned int)atoi(optarg_);
141 break;
142 case 'C':
143 nbcols = (unsigned int)atoi(optarg_);
144 break;
145 case 'v':
146 verbose = true;
147 break;
148 // add default options -c -d
149 case 'c':
150 break;
151 case 'd':
152 break;
153 default:
154 usage(argv[0], optarg_);
155 return false;
156 break;
157 }
158 }
159
160 if ((c == 1) || (c == -1)) {
161 // standalone param or error
162 usage(argv[0], NULL);
163 std::cerr << "ERROR: " << std::endl;
164 std::cerr << " Bad argument " << optarg_ << std::endl << std::endl;
165 return false;
166 }
167
168 return true;
169}
170
171vpMatrix make_random_matrix(unsigned int nbrows, unsigned int nbcols)
172{
173 vpMatrix A;
174 A.resize(nbrows, nbcols);
175
176 for (unsigned int i = 0; i < A.getRows(); i++)
177 for (unsigned int j = 0; j < A.getCols(); j++)
178 A[i][j] = (double)rand() / (double)RAND_MAX;
179 return A;
180}
181
182void create_bench(unsigned int nb_matrices, unsigned int nb_rows, unsigned int nb_cols, bool verbose,
183 std::vector<vpMatrix> &bench)
184{
185 if (verbose)
186 std::cout << "Create a bench of " << nb_matrices << " " << nb_rows << " by " << nb_cols << " matrices" << std::endl;
187 bench.clear();
188 for (unsigned int i = 0; i < nb_matrices; i++) {
189 vpMatrix M = make_random_matrix(nb_rows, nb_cols);
190 bench.push_back(M);
191 }
192}
193
194void test_det_default(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
195{
196 if (verbose)
197 std::cout << "Test determinant using default method" << std::endl;
198 // Compute inverse
199 if (verbose)
200 std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
201
202 result.resize(bench.size());
203 double t = vpTime::measureTimeMs();
204 for (unsigned int i = 0; i < bench.size(); i++) {
205 result[i] = bench[i].AtA().det();
206 }
207 time = vpTime::measureTimeMs() - t;
208}
209
210#if defined(VISP_HAVE_EIGEN3)
211void test_det_eigen3(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
212{
213 if (verbose)
214 std::cout << "Test determinant using Eigen3 3rd party" << std::endl;
215 // Compute inverse
216 if (verbose)
217 std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
218
219 result.resize(bench.size());
220 double t = vpTime::measureTimeMs();
221 for (unsigned int i = 0; i < bench.size(); i++) {
222 result[i] = bench[i].AtA().detByLUEigen3();
223 }
224 time = vpTime::measureTimeMs() - t;
225}
226#endif
227
228#if defined(VISP_HAVE_LAPACK)
229void test_det_lapack(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
230{
231 if (verbose)
232 std::cout << "Test determinant using Lapack 3rd party" << std::endl;
233 // Compute inverse
234 if (verbose)
235 std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
236
237 result.resize(bench.size());
238 double t = vpTime::measureTimeMs();
239 for (unsigned int i = 0; i < bench.size(); i++) {
240 result[i] = bench[i].AtA().detByLULapack();
241 }
242 time = vpTime::measureTimeMs() - t;
243}
244#endif
245
246#if (VISP_HAVE_OPENCV_VERSION >= 0x020101)
247void test_det_opencv(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
248{
249 if (verbose)
250 std::cout << "Test determinant using OpenCV 3rd party" << std::endl;
251 // Compute inverse
252 if (verbose)
253 std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
254
255 result.resize(bench.size());
256 double t = vpTime::measureTimeMs();
257 for (unsigned int i = 0; i < bench.size(); i++) {
258 result[i] = bench[i].AtA().detByLUOpenCV();
259 }
260 time = vpTime::measureTimeMs() - t;
261}
262#endif
263
264void save_time(const std::string &method, bool verbose, bool use_plot_file, std::ofstream &of, double time)
265{
266 if (use_plot_file)
267 of << time << "\t";
268 if (verbose || !use_plot_file) {
269 std::cout << method << time << std::endl;
270 }
271}
272
273int main(int argc, const char *argv[])
274{
275 try {
276#if defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_LAPACK) || (VISP_HAVE_OPENCV_VERSION >= 0x020101)
277 unsigned int nb_matrices = 1000;
278 unsigned int nb_iterations = 10;
279 unsigned int nb_rows = 6;
280 unsigned int nb_cols = 6;
281 bool verbose = false;
282 std::string plotfile("plot-det.csv");
283 bool use_plot_file = false;
284 std::ofstream of;
285
286 // Read the command line options
287 if (getOptions(argc, argv, nb_matrices, nb_iterations, use_plot_file, plotfile, nb_rows, nb_cols, verbose) ==
288 false) {
289 exit(-1);
290 }
291
292 if (use_plot_file) {
293 of.open(plotfile.c_str());
294 of << "iter"
295 << "\t";
296
297 of << "\"Determinant default\""
298 << "\t";
299
300#if defined(VISP_HAVE_LAPACK)
301 of << "\"Determinant Lapack\""
302 << "\t";
303#endif
304#if defined(VISP_HAVE_EIGEN3)
305 of << "\"Determinant Eigen3\""
306 << "\t";
307#endif
308#if (VISP_HAVE_OPENCV_VERSION >= 0x020101)
309 of << "\"Determinant OpenCV\""
310 << "\t";
311#endif
312 of << std::endl;
313 }
314
315 int ret = EXIT_SUCCESS;
316 for (unsigned int iter = 0; iter < nb_iterations; iter++) {
317 std::vector<vpMatrix> bench;
318 create_bench(nb_matrices, nb_rows, nb_cols, verbose, bench);
319
320 if (use_plot_file)
321 of << iter << "\t";
322
323 double time;
324
325 std::vector<double> result_default;
326 test_det_default(verbose, bench, time, result_default);
327 save_time("Determinant default: ", verbose, use_plot_file, of, time);
328
329#if defined(VISP_HAVE_LAPACK)
330 std::vector<double> result_lapack;
331 test_det_lapack(verbose, bench, time, result_lapack);
332 save_time("Determinant by Lapack: ", verbose, use_plot_file, of, time);
333#endif
334
335#if defined(VISP_HAVE_EIGEN3)
336 std::vector<double> result_eigen3;
337 test_det_eigen3(verbose, bench, time, result_eigen3);
338 save_time("Determinant by Eigen3: ", verbose, use_plot_file, of, time);
339#endif
340
341#if (VISP_HAVE_OPENCV_VERSION >= 0x020101)
342 std::vector<double> result_opencv;
343 test_det_opencv(verbose, bench, time, result_opencv);
344 save_time("Determinant by OpenCV: ", verbose, use_plot_file, of, time);
345#endif
346
347 if (use_plot_file)
348 of << std::endl;
349
350#if defined(VISP_HAVE_LAPACK) && (VISP_HAVE_OPENCV_VERSION >= 0x020101)
351 // Compare results
352 for (unsigned int i = 0; i < bench.size(); i++) {
353 if (std::fabs(result_lapack[i] - result_opencv[i]) > 1e-6) {
354 std::cout << "Determinant differ between Lapack and OpenCV: " << result_lapack[i] << " " << result_opencv[i]
355 << std::endl;
356 ret = EXIT_FAILURE;
357 }
358 }
359#endif
360#if defined(VISP_HAVE_EIGEN3) && (VISP_HAVE_OPENCV_VERSION >= 0x020101)
361 // Compare results
362 for (unsigned int i = 0; i < bench.size(); i++) {
363 if (std::fabs(result_eigen3[i] - result_opencv[i]) > 1e-6) {
364 std::cout << "Determinant differ between Eigen3 and OpenCV: " << result_eigen3[i] << " " << result_opencv[i]
365 << std::endl;
366 ret = EXIT_FAILURE;
367 }
368 }
369#endif
370#if defined(VISP_HAVE_EIGEN3) && defined(VISP_HAVE_LAPACK)
371 // Compare results
372 for (unsigned int i = 0; i < bench.size(); i++) {
373 if (std::fabs(result_eigen3[i] - result_lapack[i]) > 1e-6) {
374 std::cout << "Determinant differ between Eigen3 and Lapack: " << result_eigen3[i] << " " << result_lapack[i]
375 << std::endl;
376 ret = EXIT_FAILURE;
377 }
378 }
379#endif
380 }
381 if (use_plot_file) {
382 of.close();
383 std::cout << "Result saved in " << plotfile << std::endl;
384 }
385
386 if (ret == EXIT_SUCCESS) {
387 std::cout << "Test succeed" << std::endl;
388 } else {
389 std::cout << "Test failed" << std::endl;
390 }
391
392 return ret;
393#else
394 (void)argc;
395 (void)argv;
396 std::cout << "Test does nothing since you dont't have Lapack, Eigen3 or OpenCV 3rd party"
397 << std::endl;
398 return EXIT_SUCCESS;
399#endif
400 } catch (const vpException &e) {
401 std::cout << "Catch an exception: " << e.getStringMessage() << std::endl;
402 return EXIT_FAILURE;
403 }
404}
unsigned int getCols() const
Definition: vpArray2D.h:279
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:304
unsigned int getRows() const
Definition: vpArray2D.h:289
error that can be emited by ViSP classes.
Definition: vpException.h:72
const std::string & getStringMessage() const
Send a reference (constant) related the error message (can be empty).
Definition: vpException.cpp:92
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:154
static bool parse(int *argcPtr, const char **argv, vpArgvInfo *argTable, int flags)
Definition: vpParseArgv.cpp:69
VISP_EXPORT double measureTimeMs()