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