50 #include <visp/vpTime.h>
52 #include <visp/vpMatrix.h>
53 #include <visp/vpColVector.h>
54 #include <visp/vpParseArgv.h>
62 #define GETOPTARGS "n:i:pf:r:c:vh"
72 void usage(
const char *name,
const char *badparam)
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\
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);
86 -n <number of matrices> \n\
87 Number of matrices inverted during each test loop.\n\
89 -i <number of iterations> \n\
90 Number of iterations of the test.\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\
98 -r <number of rows>\n\
99 Number of rows of the automatically generated matrices \n\
102 -c <number of columns>\n\
103 Number of colums of the automatically generated matrices \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\
112 Print the help.\n\n");
115 fprintf(stderr,
"ERROR: \n" );
116 fprintf(stderr,
"\nBad parameter [%s]\n", badparam);
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)
137 case 'h': usage(argv[0], NULL);
return false;
break;
139 nb_matrices = (
unsigned int)atoi(optarg);
142 nb_iterations = (
unsigned int)atoi(optarg);
146 use_plot_file =
true;
149 use_plot_file =
true;
152 nbrows = (
unsigned int)atoi(optarg);
155 nbcols = (
unsigned int)atoi(optarg);
161 usage(argv[0], optarg);
166 if ((c == 1) || (c == -1)) {
168 usage(argv[0], NULL);
169 std::cerr <<
"ERROR: " << std::endl;
170 std::cerr <<
" Bad argument " << optarg << std::endl << std::endl;
177 void writeTime(
const char *name,
double time){
178 std::cout << name <<
" time=" << time <<
" ms." << std::endl;
183 vpMatrix makeRandomMatrix(
unsigned int nbrows,
unsigned int nbcols){
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;
195 main(
int argc,
const char ** argv)
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;
207 double t, qr_time, lu_time,pi_time,chol_time;
209 if (getOptions(argc, argv, nb_matrices,nb_iterations,use_plot_file,plotfile,nb_rows,nb_cols,verbose) ==
false) {
214 of.open(plotfile.c_str());
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;
223 std::cout <<
"********* generating matrices for iteration " << iter <<
"." << std::endl;
224 for(
unsigned int i=0;i<nb_matrices;i++){
228 for(cur=makeRandomMatrix(nb_rows,nb_cols);std::abs(det=cur.
AtA().
det())<.01;cur = makeRandomMatrix(nb_rows,nb_cols))
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;
233 benchCholesky.push_back(cur);
234 benchQR.push_back(cur);
235 benchLU.push_back(cur);
236 benchPseudoInverse.push_back(cur);
240 std::cout <<
"\t Inverting " << benchCholesky[0].
AtA().
getRows() <<
"x" << benchCholesky[0].AtA().getCols() <<
" matrix using cholesky decomposition." << std::endl;
242 for(
unsigned int i=0;i<nb_matrices;i++){
243 benchCholesky[i]=benchCholesky[i].AtA().inverseByCholesky()*benchCholesky[i].transpose();
248 std::cout <<
"\t Inverting " << benchLU[0].AtA().getRows() <<
"x" << benchLU[0].AtA().getCols() <<
" matrix using LU decomposition." << std::endl;
250 for(
unsigned int i=0;i<nb_matrices;i++)
251 benchLU[i] = benchLU[i].AtA().inverseByLU()*benchLU[i].transpose();
255 std::cout <<
"\t Inverting " << benchQR[0].AtA().getRows() <<
"x" << benchQR[0].AtA().getCols() <<
" matrix using QR decomposition." << std::endl;
257 for(
unsigned int i=0;i<nb_matrices;i++){
258 benchQR[i]=benchQR[i].AtA().inverseByQR()*benchQR[i].transpose();
263 std::cout <<
"\t Inverting " << benchPseudoInverse[0].AtA().getRows() <<
"x" << benchPseudoInverse[0].AtA().getCols() <<
" matrix while computing pseudo-inverse." << std::endl;
265 for(
unsigned int i=0;i<nb_matrices;i++){
266 benchPseudoInverse[i]=benchPseudoInverse[i].pseudoInverse();
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.;
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();
286 avg_err_lu_qr/=nb_matrices;
287 avg_err_lu_pi/=nb_matrices;
288 avg_err_qr_pi/=nb_matrices;
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;
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);
303 std::cout <<
"You don't have lapack installed" << std::endl;
Definition of the vpMatrix class.
void resize(const unsigned int nrows, const unsigned int ncols, const bool nullify=true)
static double measureTimeMs()
static bool parse(int *argcPtr, const char **argv, vpArgvInfo *argTable, int flags)
double det(vpDetMethod method=LU_DECOMPOSITION) const
unsigned int getCols() const
Return the number of columns of the matrix.
unsigned int getRows() const
Return the number of rows of the matrix.