ViSP  3.0.0
testMatrixInverse.cpp
1 /****************************************************************************
2  *
3  * This file is part of the ViSP software.
4  * Copyright (C) 2005 - 2015 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * ("GPL") version 2 as published by the Free Software Foundation.
9  * See the file LICENSE.txt at the root directory of this source
10  * distribution for additional information about the GNU GPL.
11  *
12  * For using ViSP with software that can not be combined with the GNU
13  * GPL, please contact Inria about acquiring a ViSP Professional
14  * Edition License.
15  *
16  * See http://visp.inria.fr for more information.
17  *
18  * This software was developed at:
19  * Inria Rennes - Bretagne Atlantique
20  * Campus Universitaire de Beaulieu
21  * 35042 Rennes Cedex
22  * France
23  *
24  * If you have questions regarding the use of this file, please contact
25  * Inria at visp@inria.fr
26  *
27  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
28  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
29  *
30  * Description:
31  * Test various inversions.
32  *
33  * Authors:
34  * Filip Novotny
35  *
36  *****************************************************************************/
37 
38 
46 #include <visp3/core/vpTime.h>
47 
48 #include <visp3/core/vpMatrix.h>
49 #include <visp3/core/vpColVector.h>
50 #include <visp3/io/vpParseArgv.h>
51 #include <vector>
52 #include <stdlib.h>
53 #include <stdio.h>
54 #include <fstream>
55 #include <iostream>
56 #include <cmath>
57 // List of allowed command line options
58 #define GETOPTARGS "cdn:i:pf:R:C:vh"
59 
60 void usage(const char *name, const char *badparam);
61 bool getOptions(int argc, const char **argv,
62  unsigned int& nb_matrices, unsigned int& nb_iterations,
63  bool& use_plot_file, std::string& plotfile,
64  unsigned int& nbrows, unsigned int& nbcols, bool& verbose);
65 void writeTime(const char *name, double time);
66 vpMatrix makeRandomMatrix(unsigned int nbrows, unsigned int nbcols);
67 
76 void usage(const char *name, const char *badparam)
77 {
78  fprintf(stdout, "\n\
79 Test matrix inversions\n\
80 using LU, QR and Cholesky methods as well as Pseudo-inverse.\n\
81 Outputs a comparison of these methods.\n\
82 \n\
83 SYNOPSIS\n\
84  %s [-n <number of matrices>] [-f <plot filename>]\n\
85  [-R <number of rows>] [-C <number of columns>]\n\
86  [-i <number of iterations>] [-p] [-h]\n", name);
87 
88  fprintf(stdout, "\n\
89 OPTIONS: Default\n\
90  -n <number of matrices> \n\
91  Number of matrices inverted during each test loop.\n\
92 \n\
93  -i <number of iterations> \n\
94  Number of iterations of the test.\n\
95 \n\
96  -f <plot filename> \n\
97  Set output path for plot output.\n\
98  The plot logs the times of \n\
99  the different inversion methods: \n\
100  QR,LU,Cholesky and Pseudo-inverse.\n\
101 \n\
102  -R <number of rows>\n\
103  Number of rows of the automatically generated matrices \n\
104  we test on.\n\
105 \n\
106  -C <number of columns>\n\
107  Number of colums of the automatically generated matrices \n\
108  we test on.\n\
109 \n\
110  -p \n\
111  Plot into filename in the gnuplot format. \n\
112  If this option is used, tests results will be logged \n\
113  into a filename specified with -f.\n\
114 \n\
115  -h\n\
116  Print the help.\n\n");
117 
118  if (badparam) {
119  fprintf(stderr, "ERROR: \n" );
120  fprintf(stderr, "\nBad parameter [%s]\n", badparam);
121  }
122 
123 }
131 bool getOptions(int argc, const char **argv,
132  unsigned int& nb_matrices, unsigned int& nb_iterations,
133  bool& use_plot_file, std::string& plotfile,
134  unsigned int& nbrows, unsigned int& nbcols, bool& verbose)
135 {
136  const char *optarg_;
137  int c;
138  while ((c = vpParseArgv::parse(argc, argv, GETOPTARGS, &optarg_)) > 1) {
139 
140  switch (c) {
141  case 'h': usage(argv[0], NULL); return false; break;
142  case 'n':
143  nb_matrices = (unsigned int)atoi(optarg_);
144  break;
145  case 'i':
146  nb_iterations = (unsigned int)atoi(optarg_);
147  break;
148  case 'f':
149  plotfile = optarg_;
150  use_plot_file = true;
151  break;
152  case 'p':
153  use_plot_file = true;
154  break;
155  case 'R':
156  nbrows = (unsigned int)atoi(optarg_);
157  break;
158  case 'C':
159  nbcols = (unsigned int)atoi(optarg_);
160  break;
161  case 'v':
162  verbose = true;
163  break;
164  // add default options -c -d
165  case 'c':
166  break;
167  case 'd':
168  break;
169  default:
170  usage(argv[0], optarg_);
171  return false; break;
172  }
173  }
174 
175  if ((c == 1) || (c == -1)) {
176  // standalone param or error
177  usage(argv[0], NULL);
178  std::cerr << "ERROR: " << std::endl;
179  std::cerr << " Bad argument " << optarg_ << std::endl << std::endl;
180  return false;
181  }
182 
183  return true;
184 }
185 
186 void writeTime(const char *name, double time)
187 {
188  std::cout << name << " time=" << time << " ms." << std::endl;
189 }
190 
191 vpMatrix makeRandomMatrix(unsigned int nbrows, unsigned int nbcols)
192 {
193  vpMatrix A;
194  A.resize(nbrows,nbcols);
195 
196  for (unsigned int i=0 ; i < A.getRows() ; i++)
197  for (unsigned int j=0 ; j < A.getCols() ; j++)
198  A[i][j] = (double)rand()/(double)RAND_MAX;
199  return A;
200 }
201 
202 
203 int
204 main(int argc, const char ** argv)
205 {
206 #ifdef VISP_HAVE_LAPACK_C
207  try {
208  unsigned int nb_matrices=1000;
209  unsigned int nb_iterations=10;
210  unsigned int nb_rows = 6;
211  unsigned int nb_cols = 6;
212  bool verbose = false;
213  std::string plotfile("plot.txt");
214  bool use_plot_file=false;
215  std::ofstream of;
216 
217  double t, qr_time, lu_time,pi_time,chol_time;
218  // Read the command line options
219  if (getOptions(argc, argv, nb_matrices,nb_iterations,use_plot_file,plotfile,nb_rows,nb_cols,verbose) == false) {
220  exit (-1);
221  }
222 
223  if(use_plot_file){
224  of.open(plotfile.c_str());
225  }
226 
227  for(unsigned int iter=0;iter<nb_iterations;iter++){
228  std::vector<vpMatrix> benchQR;
229  std::vector<vpMatrix> benchLU;
230  std::vector<vpMatrix> benchCholesky;
231  std::vector<vpMatrix> benchPseudoInverse;
232  if(verbose)
233  std::cout << "********* generating matrices for iteration " << iter << "." << std::endl;
234  for(unsigned int i=0;i<nb_matrices;i++){
235  vpMatrix cur;
236  double det=0.;
237  //don't put singular matrices in the benchmark
238  for(cur=makeRandomMatrix(nb_rows,nb_cols);std::abs(det=cur.AtA().det())<.01;cur = makeRandomMatrix(nb_rows,nb_cols))
239  if(verbose){
240  std::cout << "Generated random matrix A*tA=" << std::endl << cur.AtA() << std::endl;
241  std::cout << "generated random matrix not invertibleL: det="<<det<< ". Retrying..." << std::endl;
242  }
243  benchCholesky.push_back(cur);
244  benchQR.push_back(cur);
245  benchLU.push_back(cur);
246  benchPseudoInverse.push_back(cur);
247  }
248 
249  if(verbose)
250  std::cout << "\t Inverting " << benchCholesky[0].AtA().getRows() << "x" << benchCholesky[0].AtA().getCols() << " matrix using cholesky decomposition." << std::endl;
251  t = vpTime::measureTimeMs() ;
252  for(unsigned int i=0;i<nb_matrices;i++){
253  benchCholesky[i]=benchCholesky[i].AtA().inverseByCholesky()*benchCholesky[i].transpose();
254  }
255  chol_time = vpTime::measureTimeMs() - t ;
256 
257  if(verbose)
258  std::cout << "\t Inverting " << benchLU[0].AtA().getRows() << "x" << benchLU[0].AtA().getCols() << " matrix using LU decomposition." << std::endl;
259  t = vpTime::measureTimeMs() ;
260  for(unsigned int i=0;i<nb_matrices;i++)
261  benchLU[i] = benchLU[i].AtA().inverseByLU()*benchLU[i].transpose();
262  lu_time = vpTime::measureTimeMs() -t ;
263 
264  if(verbose)
265  std::cout << "\t Inverting " << benchQR[0].AtA().getRows() << "x" << benchQR[0].AtA().getCols() << " matrix using QR decomposition." << std::endl;
266  t = vpTime::measureTimeMs() ;
267  for(unsigned int i=0;i<nb_matrices;i++){
268  benchQR[i]=benchQR[i].AtA().inverseByQR()*benchQR[i].transpose();
269  }
270  qr_time = vpTime::measureTimeMs() - t ;
271 
272  if(verbose)
273  std::cout << "\t Inverting " << benchPseudoInverse[0].AtA().getRows() << "x" << benchPseudoInverse[0].AtA().getCols() << " matrix while computing pseudo-inverse." << std::endl;
274  t = vpTime::measureTimeMs() ;
275  for(unsigned int i=0;i<nb_matrices;i++){
276  benchPseudoInverse[i]=benchPseudoInverse[i].pseudoInverse();
277  }
278  pi_time = vpTime::measureTimeMs() - t ;
279 
280  double avg_err_lu_qr=0.;
281  double avg_err_lu_pi=0.;
282  double avg_err_lu_chol=0.;
283  double avg_err_qr_pi=0.;
284  double avg_err_qr_chol=0.;
285  double avg_err_pi_chol=0.;
286 
287  for(unsigned int i=0;i<nb_matrices;i++){
288  avg_err_lu_qr+= (benchQR[i]-benchLU[i]).euclideanNorm();
289  avg_err_lu_pi+= (benchPseudoInverse[i]-benchLU[i]).euclideanNorm();
290  avg_err_qr_pi+= (benchPseudoInverse[i]-benchQR[i]).euclideanNorm();
291  avg_err_qr_chol+= (benchCholesky[i]-benchQR[i]).euclideanNorm();
292  avg_err_lu_chol+= (benchCholesky[i]-benchLU[i]).euclideanNorm();
293  avg_err_pi_chol+= (benchCholesky[i]-benchPseudoInverse[i]).euclideanNorm();
294  }
295 
296  avg_err_lu_qr/=nb_matrices;
297  avg_err_lu_pi/=nb_matrices;
298  avg_err_qr_pi/=nb_matrices;
299 
300  if(use_plot_file){
301  of << iter << "\t" << lu_time << "\t" << qr_time << "\t" << pi_time << "\t" << chol_time << "\t" << avg_err_lu_qr << "\t" << avg_err_qr_pi << "\t" << avg_err_lu_pi << "\t" << avg_err_qr_chol << "\t" << avg_err_lu_chol << "\t" << avg_err_pi_chol << std::endl;
302  }
303  if(verbose || !use_plot_file){
304  writeTime("LU",lu_time);
305  writeTime("QR",qr_time);
306  writeTime("Pseudo-inverse",pi_time);
307  writeTime("Cholesky",chol_time);
308  }
309  }
310  return 0;
311  }
312  catch(vpException e) {
313  std::cout << "Catch an exception: " << e << std::endl;
314  return 1;
315  }
316 
317 #else
318  (void)argc;
319  (void)argv;
320  std::cout << "You don't have lapack installed" << std::endl;
321  return 0;
322 #endif
323 }
324 
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:92
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true)
Definition: vpArray2D.h:167
error that can be emited by ViSP classes.
Definition: vpException.h:73
unsigned int getCols() const
Return the number of columns of the 2D array.
Definition: vpArray2D.h:154
static bool parse(int *argcPtr, const char **argv, vpArgvInfo *argTable, int flags)
Definition: vpParseArgv.cpp:76
vpMatrix AtA() const
Definition: vpMatrix.cpp:391
unsigned int getRows() const
Return the number of rows of the 2D array.
Definition: vpArray2D.h:152
double det(vpDetMethod method=LU_DECOMPOSITION) const
Definition: vpMatrix.cpp:3297
VISP_EXPORT double measureTimeMs()
Definition: vpTime.cpp:93