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