matrixwrapper_test.cpp
Go to the documentation of this file.
1 // Copyright (C) 2007 Wim Meeussen <wim.meeussen@mech.kuleuven.be>
2 // Tinne De Laet<first DOT last AT mech.kuleuven.be>
3 //
4 // This program is free software; you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation; either version 2 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
17 //
18 
19 
20 #include "matrixwrapper_test.hpp"
21 #include "approxEqual.hpp"
22 
23 
24 // Registers the fixture into the 'registry'
26 
27 using namespace MatrixWrapper;
28 
29 
30 void
32 {
33 }
34 
35 
36 void
38 {
39 }
40 
41 void
43 {
44  double epsilon = 0.00001;
45 
46  unsigned int one = 1;
47  unsigned int r = 4;
48  unsigned int c = 3;
49 
50  // reference numbers
51  vector< vector<double> > REF;
52  for (unsigned int i=0; i<r+c; i++){
53  vector<double> row;
54  for (unsigned int j=0; j<r+c; j++)
55  row.push_back( (r+c)*i+j );
56  REF.push_back(row);
57  }
58 
59  // TEST DIMENSIONS
60  //
61  // TEST DIMENSIONS MATRIX
62  Matrix Am(r,c);
63  // rows()
64  CPPUNIT_ASSERT_EQUAL(Am.rows(), r);
65  // columns()
66  CPPUNIT_ASSERT_EQUAL(Am.columns(), c);
67  // TEST DIMENSIONS SYMMETRICMATRIX
68  SymmetricMatrix As(r);
69  // rows()
70  CPPUNIT_ASSERT_EQUAL(As.rows(), r);
71  // columns()
72  CPPUNIT_ASSERT_EQUAL(As.columns(), r);
73  // TEST DIMENSIONS COLUMNVECTOR
74  ColumnVector Ac(r);
75  // rows()
76  CPPUNIT_ASSERT_EQUAL(Ac.rows(), r);
77  // columns()
78  CPPUNIT_ASSERT_EQUAL(Ac.columns(), one);
79  // TEST DIMENSIONS ROWVECTOR
80  RowVector Ar(c);
81  // rows()
82  CPPUNIT_ASSERT_EQUAL(Ar.rows(), one);
83  // columns()
84  CPPUNIT_ASSERT_EQUAL(Ar.columns(), c);
85 
86  // test operator = double
87  double v = 3.5;
88  Am = v; As = v; Ac = v; Ar = v;
89  // MATRIX
90  for (unsigned int i=0; i<r; i++)
91  for (unsigned int j=0; j<c; j++)
92  CPPUNIT_ASSERT_EQUAL(Am(i+1,j+1), v);
93  // SYMMETRICMATRIX
94  for (unsigned int i=0; i<r; i++)
95  for (unsigned int j=0; j<r; j++)
96  CPPUNIT_ASSERT_EQUAL(As(i+1,j+1), v);
97  // COLUMNVECTOR
98  for (unsigned int i=0; i<r; i++)
99  CPPUNIT_ASSERT_EQUAL(Ac(i+1), v);
100  // ROWVECTOR
101  for (unsigned int i=0; i<c; i++)
102  CPPUNIT_ASSERT_EQUAL(Ar(i+1), v);
103 
104  // test operator ()
105  // MATRIX
106  Matrix Bm(r,c);
107  for (unsigned int i=0; i<r; i++){
108  for (unsigned int j=0; j<c; j++){
109  Bm(i+1,j+1) = REF[i][j];
110  CPPUNIT_ASSERT_EQUAL(Bm(i+1,j+1), REF[i][j]);
111  }
112  }
113 
114  // SYMMETRICMATRIX
115  SymmetricMatrix Bs(r);
116  for (unsigned int i=0; i<r; i++){ // fill in upper triangle
117  for (unsigned int j=0; j<=i; j++){
118  Bs(j+1,i+1) = REF[i][j];
119  CPPUNIT_ASSERT_EQUAL(Bs(i+1,j+1), REF[i][j]);
120  CPPUNIT_ASSERT_EQUAL(Bs(j+1,i+1), REF[i][j]);
121  }
122  }
123  for (unsigned int i=0; i<r; i++){ // fill in lower triangle
124  for (unsigned int j=0; j<=i; j++){
125  Bs(i+1,j+1) = REF[i][j];
126  CPPUNIT_ASSERT_EQUAL(Bs(i+1,j+1), REF[i][j]);
127  CPPUNIT_ASSERT_EQUAL(Bs(j+1,i+1), REF[i][j]);
128  }
129  }
130  // COLUMNVECTOR
131  ColumnVector Bc(r);
132  for (unsigned int i=0; i<r; i++){
133  Bc(i+1) = REF[0][i];
134  CPPUNIT_ASSERT_EQUAL(Bc(i+1), REF[0][i]);
135  }
136  // ROWVECTOR
137  RowVector Br(c);
138  for (unsigned int i=0; i<c; i++){
139  Br(i+1) = REF[0][i];
140  CPPUNIT_ASSERT_EQUAL(Br(i+1), REF[0][i]);
141  }
142 
143  // test operator =
144  // MATRIX
145  Am = Bm;
146  for (unsigned int i=0; i<r; i++){
147  for (unsigned int j=0; j<c; j++){
148  CPPUNIT_ASSERT_EQUAL(Am(i+1,j+1), REF[i][j]);
149  CPPUNIT_ASSERT_EQUAL(Bm(i+1,j+1), REF[i][j]);
150  }
151  }
152  // SYMMETRICMATRIX
153  As = Bs;
154  for (unsigned int i=0; i<r; i++){
155  for (unsigned int j=0; j<=i; j++){
156  CPPUNIT_ASSERT_EQUAL(As(i+1,j+1), REF[i][j]);
157  CPPUNIT_ASSERT_EQUAL(Bs(i+1,j+1), REF[i][j]);
158  CPPUNIT_ASSERT_EQUAL(As(j+1,i+1), REF[i][j]);
159  CPPUNIT_ASSERT_EQUAL(Bs(j+1,i+1), REF[i][j]);
160  }
161  }
162  // COLUMNVECTOR
163  Ac = Bc;
164  for (unsigned int i=0; i<r; i++){
165  CPPUNIT_ASSERT_EQUAL(Ac(i+1), REF[0][i]);
166  CPPUNIT_ASSERT_EQUAL(Bc(i+1), REF[0][i]);
167  }
168  // ROWVECTOR
169  Ar = Br;
170  for (unsigned int i=0; i<c; i++){
171  CPPUNIT_ASSERT_EQUAL(Ar(i+1), REF[0][i]);
172  CPPUNIT_ASSERT_EQUAL(Br(i+1), REF[0][i]);
173  }
174 
175 
176 
177  // test resize
178  // MATRIX
179  Matrix Km(r+2,c+2); Km = v;
180  Matrix Km_resize(r,c); Km_resize = v;
181  CPPUNIT_ASSERT_EQUAL(Km.rows(), r+2);
182  CPPUNIT_ASSERT_EQUAL(Km.columns(), c+2);
183  Km.resize(r,c);
184  CPPUNIT_ASSERT_EQUAL(Km.rows(), r);
185  CPPUNIT_ASSERT_EQUAL(Km.columns(), c);
186  CPPUNIT_ASSERT_EQUAL(Km, Km_resize);
187  // SYMMETRICMATRIX
188  SymmetricMatrix Ks(r+2); Ks = v;
189  SymmetricMatrix Ks_resize(r); Ks_resize = v;
190  CPPUNIT_ASSERT_EQUAL(Ks.rows(), r+2);
191  CPPUNIT_ASSERT_EQUAL(Ks.columns(), r+2);
192 // Ks.resize(r);
193 // CPPUNIT_ASSERT_EQUAL(Ks.rows(), r);
194 // CPPUNIT_ASSERT_EQUAL(Ks.columns(), r);
195 // CPPUNIT_ASSERT_EQUAL(Ks, Ks_resize);
196  // COLUMNVECTOR
197  ColumnVector Kc(r+2); Kc = v;
198  ColumnVector Kc_resize(r); Kc_resize = v;
199  CPPUNIT_ASSERT_EQUAL(Kc.rows(), r+2);
200  Kc.resize(r);
201  CPPUNIT_ASSERT_EQUAL(Kc.rows(), r);
202  CPPUNIT_ASSERT_EQUAL(Kc, Kc_resize);
203  // ROWVECTOR
204  RowVector Kr(c+2); Kr = v;
205  RowVector Kr_resize(c); Kr_resize = v;
206  CPPUNIT_ASSERT_EQUAL(Kr.columns(), c+2);
207  Kr.resize(c);
208  CPPUNIT_ASSERT_EQUAL(Kr.columns(), c);
209  CPPUNIT_ASSERT_EQUAL(Kr, Kr_resize);
210 
211  // test operator ==
212  // MATRIX
213  Matrix Bm_eq; Bm_eq = Bm;
214  CPPUNIT_ASSERT_EQUAL(Bm_eq == Bm, true);
215  Bm(1,1) = Bm(1,1) + v;
216  CPPUNIT_ASSERT_EQUAL(Bm == Bm_eq, false);
217  Matrix Bm_res(Bm.rows(), Bm.columns()+1);
218  CPPUNIT_ASSERT_EQUAL(Bm == Bm_res, false);
219  Bm_res = Bm;
220  CPPUNIT_ASSERT_EQUAL(Bm == Bm_res, true);
221  // SYMMETRICMATRIX
222  SymmetricMatrix Bs_eq; Bs_eq = Bs;
223  CPPUNIT_ASSERT_EQUAL(Bs_eq == Bs, true);
224  Bs(1,1) = Bs(1,1) + v;
225  CPPUNIT_ASSERT_EQUAL(Bs == Bs_eq, false);
226  SymmetricMatrix Bs_res(Bs.columns()+1);
227  CPPUNIT_ASSERT_EQUAL(Bs == Bs_res, false);
228  Bs_res = Bs;
229  CPPUNIT_ASSERT_EQUAL(Bs == Bs_res, true);
230  // COLUMNVECTOR
231  ColumnVector Bc_eq; Bc_eq = Bc;
232  CPPUNIT_ASSERT_EQUAL(Bc_eq == Bc, true);
233  Bc(1) = Bc(1) + v;
234  CPPUNIT_ASSERT_EQUAL(Bc == Bc_eq, false);
235  ColumnVector Bc_res(Bc.rows()+1);
236  CPPUNIT_ASSERT_EQUAL(Bc == Bc_res, false);
237  Bc_res = Bc;
238  CPPUNIT_ASSERT_EQUAL(Bc == Bc_res, true);
239  // ROWVECTOR
240  RowVector Br_eq; Br_eq = Br;
241  CPPUNIT_ASSERT_EQUAL(Br_eq == Br, true);
242  Br(1) = Br(1) + v;
243  CPPUNIT_ASSERT_EQUAL(Br == Br_eq, false);
244  RowVector Br_res(Br.columns()+1);
245  CPPUNIT_ASSERT_EQUAL(Br == Br_res, false);
246  Br_res = Br;
247  CPPUNIT_ASSERT_EQUAL(Br == Br_res, true);
248 
249 
250  // test transpose
251  // MATRIX
252  Matrix Am_trans = Am.transpose();
253  CPPUNIT_ASSERT_EQUAL(Am.rows(), Am_trans.columns());
254  CPPUNIT_ASSERT_EQUAL(Am.columns(), Am_trans.rows());
255  for (unsigned int i=0; i<r; i++){
256  for (unsigned int j=0; j<c; j++){
257  CPPUNIT_ASSERT_EQUAL(Am(i+1,j+1), REF[i][j]);
258  CPPUNIT_ASSERT_EQUAL(Am_trans(j+1,i+1), REF[i][j]);
259  }
260  }
261  // SYMMETRICMATRIX
262  SymmetricMatrix As_trans = As.transpose();
263  CPPUNIT_ASSERT_EQUAL(As.rows(), As_trans.columns());
264  CPPUNIT_ASSERT_EQUAL(As.columns(), As_trans.rows());
265  for (unsigned int i=0; i<r; i++){
266  for (unsigned int j=0; j<=i; j++){
267  CPPUNIT_ASSERT_EQUAL(As(i+1,j+1), REF[i][j]);
268  CPPUNIT_ASSERT_EQUAL(As_trans(j+1,i+1), REF[i][j]);
269  CPPUNIT_ASSERT_EQUAL(As(j+1,i+1), REF[i][j]);
270  CPPUNIT_ASSERT_EQUAL(As_trans(i+1,j+1), REF[i][j]);
271  }
272  }
273  // COLUMNVECTOR
274  RowVector Ac_trans = Ac.transpose();
275  CPPUNIT_ASSERT_EQUAL(Ac.rows(), Ac_trans.columns());
276  CPPUNIT_ASSERT_EQUAL(Ac.columns(), Ac_trans.rows());
277  for (unsigned int i=0; i<r; i++){
278  CPPUNIT_ASSERT_EQUAL(Ac(i+1), REF[0][i]);
279  CPPUNIT_ASSERT_EQUAL(Ac_trans(i+1), REF[0][i]);
280  }
281  // ROWVECTOR
282  ColumnVector Ar_trans = Ar.transpose();
283  CPPUNIT_ASSERT_EQUAL(Ar.rows(), Ar_trans.columns());
284  CPPUNIT_ASSERT_EQUAL(Ar.columns(), Ar_trans.rows());
285  for (unsigned int i=0; i<c; i++){
286  CPPUNIT_ASSERT_EQUAL(Ar(i+1), REF[0][i]);
287  CPPUNIT_ASSERT_EQUAL(Ar_trans(i+1), REF[0][i]);
288  }
289 
290  // test sub matrix
291  // MATRIX
292  Matrix Am_sub = Am.sub(2,r,2,c);
293  CPPUNIT_ASSERT_EQUAL(Am_sub.rows(), r-1);
294  CPPUNIT_ASSERT_EQUAL(Am_sub.columns(), c-1);
295  for (unsigned int i=0; i<c-1; i++){
296  for (unsigned int j=0; j<c-1; j++){
297  CPPUNIT_ASSERT_EQUAL(Am_sub(i+1,j+1), Am(i+2,j+2));
298  CPPUNIT_ASSERT_EQUAL(Am_sub(i+1,j+1), REF[i+1][j+1]);
299  }
300  }
301  // SYMMETRICMATRIX
302  Matrix As_sub = As.sub(2,c,2,c);
303  CPPUNIT_ASSERT_EQUAL(As_sub.rows(), c-1);
304  CPPUNIT_ASSERT_EQUAL(As_sub.columns(), c-1);
305  for (unsigned int i=0; i<c-1; i++){
306  for (unsigned int j=0; j<=i; j++){
307  CPPUNIT_ASSERT_EQUAL(As_sub(i+1,j+1), As(i+2,j+2));
308  CPPUNIT_ASSERT_EQUAL(As_sub(i+1,j+1), REF[i+1][j+1]);
309  CPPUNIT_ASSERT_EQUAL(As_sub(j+1,i+1), As(i+2,j+2));
310  CPPUNIT_ASSERT_EQUAL(As_sub(j+1,i+1), REF[i+1][j+1]);
311  }
312  }
313  // COLUMNVECTOR
314  ColumnVector Ac_sub = Ac.sub(2,c);
315  CPPUNIT_ASSERT_EQUAL(Ac_sub.rows(), c-1);
316  CPPUNIT_ASSERT_EQUAL(Ac_sub.columns(), one);
317  for (unsigned int i=0; i<c-1; i++){
318  CPPUNIT_ASSERT_EQUAL(Ac_sub(i+1), Ac(i+2));
319  CPPUNIT_ASSERT_EQUAL(Ac_sub(i+1), REF[0][i+1]);
320  }
321  // ROWVECTOR
322  RowVector Ar_sub = Ar.sub(2,r-1);
323  CPPUNIT_ASSERT_EQUAL(Ar_sub.rows(), one);
324  CPPUNIT_ASSERT_EQUAL(Ar_sub.columns(), r-2);
325  for (unsigned int i=0; i<r-2; i++){
326  CPPUNIT_ASSERT_EQUAL(Ar_sub(i+1), Ar(i+2));
327  CPPUNIT_ASSERT_EQUAL(Ar_sub(i+1), REF[0][i+1]);
328  }
329 
330  // test operator *
331  // MATRIX * MATRIX
332  Matrix Cm = Am * Am_trans;
333  Matrix Cm_check(r,c);
334  for (unsigned int i=0; i<r; i++){
335  for (unsigned int j=0; j<c; j++){
336  // test if original elements were maintained
337  CPPUNIT_ASSERT_EQUAL(Am(i+1,j+1), REF[i][j]);
338  CPPUNIT_ASSERT_EQUAL(Am_trans(j+1,i+1), REF[i][j]);
339  // test correct multiplication
340  double sum = 0.0;
341  for (unsigned int t=0; t<c; t++){
342  sum += Am(i+1,t+1) * Am_trans(t+1,j+1);
343  }
344  Cm_check(i+1,j+1) = sum;
345  CPPUNIT_ASSERT_EQUAL(Cm(i+1,j+1), Cm_check(i+1,j+1));
346  }
347  }
348 
349  // MATRIX * DOUBLE
350  Cm = Am * v;
351  for (unsigned int i=0; i<r; i++){
352  for (unsigned int j=0; j<c; j++){
353  // test if original elements were maintained
354  CPPUNIT_ASSERT_EQUAL(Am(i+1,j+1), REF[i][j]);
355  CPPUNIT_ASSERT_EQUAL(Am_trans(j+1,i+1), REF[i][j]);
356  // test correct multiplication
357  CPPUNIT_ASSERT_EQUAL(Cm(i+1,j+1), REF[i][j] * v);
358  }
359  }
360 
361  // SYMMETRICMATRIX * SYMMETRICMATRIX
362  Matrix Cs_check(r,r);
363  Matrix Cs = As * As_trans;
364  for (unsigned int i=0; i<r; i++){
365  for (unsigned int j=0; j<=i; j++){
366  // test if original elements were maintained
367  CPPUNIT_ASSERT_EQUAL(As(i+1,j+1), REF[i][j]);
368  CPPUNIT_ASSERT_EQUAL(As_trans(j+1,i+1), REF[i][j]);
369  CPPUNIT_ASSERT_EQUAL(As(j+1,i+1), REF[i][j]);
370  CPPUNIT_ASSERT_EQUAL(As_trans(i+1,j+1), REF[i][j]);
371  // test correct multiplication
372  double sum = 0.0;
373  for (unsigned int t=0; t<r; t++){
374  sum += As(i+1,t+1) * As_trans(t+1,j+1);
375  }
376  Cs_check(i+1,j+1) = sum;
377  CPPUNIT_ASSERT_EQUAL(Cs(i+1,j+1), Cs_check(i+1,j+1));
378  }
379  }
380 
381  // SYMMETRICMATRIX * DOUBLE
382  Cs = As * v;
383  for (unsigned int i=0; i<r; i++){
384  for (unsigned int j=0; j<=i; j++){
385  // test if original elements were maintained
386  CPPUNIT_ASSERT_EQUAL(As(i+1,j+1), REF[i][j]);
387  CPPUNIT_ASSERT_EQUAL(As(j+1,i+1), REF[i][j]);
388  // test correct multiplication
389  CPPUNIT_ASSERT_EQUAL(Cs(i+1,j+1), REF[i][j] * v);
390  CPPUNIT_ASSERT_EQUAL(Cs(j+1,i+1), REF[i][j] * v);
391  }
392  }
393 
394  // MATRIX * SYMMETRICMATRIX
395  // TODO: not implemented?
396 
397  // SYMMETRICMATRIX * MATRIX
398  Matrix Csm_check(r,c);
399  Matrix Csm = As * Am;
400  for (unsigned int i=0; i<r; i++){
401  for (unsigned int j=0; j<c; j++){
402  // test if original elements were maintained
403  if (j<=i){
404  CPPUNIT_ASSERT_EQUAL(As(i+1,j+1), REF[i][j]);
405  CPPUNIT_ASSERT_EQUAL(As(j+1,i+1), REF[i][j]);
406  }
407  CPPUNIT_ASSERT_EQUAL(Am(i+1,j+1), REF[i][j]);
408  // test correct multiplication
409  double sum = 0.0;
410  for (unsigned int t=0; t<r; t++){
411  sum += As(i+1,t+1) * Am(t+1,j+1);
412  }
413  Csm_check(i+1,j+1) = sum;
414  CPPUNIT_ASSERT_EQUAL(Csm(i+1,j+1), Csm_check(i+1,j+1));
415  }
416  }
417 
418  // COLUMNVECTOR * DOUBLE
419  ColumnVector Cc = Ac * v;
420  for (unsigned int i=0; i<r; i++){
421  // test if original elements were maintained
422  CPPUNIT_ASSERT_EQUAL(Ac(i+1), REF[0][i]);
423  // test correct multiplication
424  CPPUNIT_ASSERT_EQUAL(Cc(i+1), REF[0][i] * v);
425  }
426 
427  // ROWVECTOR * DOUBLE
428  RowVector Cr = Ar * v;
429  for (unsigned int i=0; i<c; i++){
430  // test if original elements were maintained
431  CPPUNIT_ASSERT_EQUAL(Ar(i+1), REF[0][i]);
432  // test correct multiplication
433  CPPUNIT_ASSERT_EQUAL(Cr(i+1), REF[0][i] * v);
434  }
435 
436  // COLUMNVECTOR * ROWVECTOR
437  Matrix Ccr = Ac * Ar;
438  Matrix Ccr_check(r,c);
439  for (unsigned int j=0; j<c; j++) CPPUNIT_ASSERT_EQUAL(Ar(j+1), REF[0][j]);
440  for (unsigned int i=0; i<r; i++){
441  // test if original elements were maintained
442  CPPUNIT_ASSERT_EQUAL(Ac(i+1), REF[0][i]);
443  for (unsigned int j=0; j<c; j++){
444  // test correct multiplication
445  Ccr_check(i+1,j+1) = Ac(i+1) * Ar(j+1);
446  CPPUNIT_ASSERT_EQUAL(Ccr(i+1,j+1), Ccr_check(i+1,j+1));
447  }
448  }
449 
450  // ROWVECTOR * COLUMNVECTOR
451  double rc = Ac_trans * Ac;
452  double rc_check;
453  // test if original elements were maintained
454  for (unsigned int j=0; j<c; j++) CPPUNIT_ASSERT_EQUAL(Ac_trans(j+1), REF[0][j]);
455  for (unsigned int i=0; i<r; i++) CPPUNIT_ASSERT_EQUAL(Ac(i+1), REF[0][i]);
456  // test correct multiplication
457  double sum = 0.0;
458  for (unsigned int t=0; t<r; t++){
459  sum += Ac_trans(t+1) * Ac(t+1);
460  }
461  rc_check = sum;
462  CPPUNIT_ASSERT_EQUAL(rc, rc_check);
463 
464  // ROWVECTOR * MATRIX
465  // TODO: only implemented for lti
466  //RowVector Cr2= Ar * Am;
467  //Matrix Cr2_check(r);
468  //for (unsigned int j=0; j<c; j++){
469  // // test if original elements were maintained
470  // CPPUNIT_ASSERT_EQUAL(Ar(j+1), REF[0][j]);
471  // for (unsigned int i=0; i<r; i++){
472  // // test if original elements were maintained
473  // CPPUNIT_ASSERT_EQUAL(Am(i+1,j+1), REF[i][j]);
474  // }
475  //}
476  //for (unsigned int i=0; i<r; i++){
477  // // test correct multiplication
478  // double sum = 0.0;
479  // for (unsigned int t=0; t<c; t++){
480  // sum += Ar(t+1) * Am(t+1,j+1);
481  // }
482  // Cr2_check(i+1) = sum;
483  // CPPUNIT_ASSERT_EQUAL(Cr2(i+1), Cr2_check(i+1));
484  //}
485 
486  // MATRIX * COLUMNVECTOR
487  ColumnVector Cc2= Am_trans * Ac;
488  ColumnVector Cc2_check(c);
489  for (unsigned int j=0; j<r; j++){
490  // test if original elements were maintained
491  CPPUNIT_ASSERT_EQUAL(Ac(j+1), REF[0][j]);
492  for (unsigned int i=0; i<c; i++){
493  // test if original elements were maintained
494  CPPUNIT_ASSERT_EQUAL(Am_trans(i+1,j+1), REF[j][i]);
495  }
496  }
497  for (unsigned int i=0; i<c; i++){
498  // test correct multiplication
499  double sum = 0.0;
500  for (unsigned int t=0; t<r; t++){
501  sum += Am_trans(i+1,t+1) * Ac(t+1);
502  }
503  Cc2_check(i+1) = sum;
504  CPPUNIT_ASSERT_EQUAL(Cc2(i+1), Cc2_check(i+1));
505  }
506 
507 
508  // test operator /
509  Cm = Am / v;
510  for (unsigned int i=0; i<r; i++){
511  for (unsigned int j=0; j<c; j++){
512  // CPPUNIT_ASSERT_EQUAL(Cm(i+1,j+1), REF[i][j] / v);
513  CPPUNIT_ASSERT_EQUAL(approxEqual(Cm(i+1,j+1), REF[i][j] / v,epsilon),true);
514  }
515  }
516  Cs = As / v;
517  for (unsigned int i=0; i<r; i++){
518  for (unsigned int j=0; j<=i; j++){
519  CPPUNIT_ASSERT_EQUAL(Cs(i+1,j+1), REF[i][j] / v);
520  CPPUNIT_ASSERT_EQUAL(Cs(j+1,i+1), REF[i][j] / v);
521  }
522  }
523  Cc = Ac / v;
524  for (unsigned int i=0; i<r; i++)
525  CPPUNIT_ASSERT_EQUAL(Cc(i+1), REF[0][i] / v);
526  Cr = Ar / v;
527  for (unsigned int i=0; i<c; i++)
528  CPPUNIT_ASSERT_EQUAL(Cr(i+1), REF[0][i] / v);
529 
530  // test inverse
531  Matrix Rm(c,c);
532  Rm(1,1) = 3; Rm(1,2) = 3; Rm(1,3) = 3;
533  Rm(2,1) = 5; Rm(2,2) = 2; Rm(2,3) = 9;
534  Rm(3,1) = 9; Rm(3,2) = 7; Rm(3,3) = 0;
535  Matrix Rm_inv = Rm.inverse();
536  Matrix Im(c,c); Im = 0;
537  for (unsigned int i=0; i<c; i++)
538  Im(i+1,i+1) = 1;
539  Matrix Im_test = Rm * Rm_inv;
540  CPPUNIT_ASSERT_EQUAL(approxEqual(Im_test, Im, epsilon),true);
541 
542  Matrix Rs(c,c);
543  Rs(1,1) = 3; Rs(1,2) = 5; Rs(1,3) = 3;
544  Rs(2,1) = 5; Rs(2,2) = 2; Rs(2,3) = 7;
545  Rs(3,1) = 3; Rs(3,2) = 7; Rs(3,3) = 0;
546  Matrix Rs_inv = Rs.inverse();
547  Matrix Is(c,c); Is = 0;
548  for (unsigned int i=0; i<c; i++)
549  Is(i+1,i+1) = 1;
550  Matrix Is_test = Rs * Rs_inv;
551  CPPUNIT_ASSERT_EQUAL(approxEqual(Is_test, Is, epsilon),true);
552 
553  // test determinant
554  CPPUNIT_ASSERT_EQUAL(approxEqual(Rm.determinant(), 105, epsilon),true);
555  CPPUNIT_ASSERT_EQUAL(approxEqual(Rs.determinant(), 45, epsilon),true);
556 
557  // test cholesky
558  SymmetricMatrix Ps(c);
559  Ps(1,1) = 3; Ps(1,2) = 2; Ps(1,3) = 1;
560  Ps(2,1) = 2; Ps(2,2) = 2; Ps(2,3) = 1;
561  Ps(3,1) = 1; Ps(3,2) = 1; Ps(3,3) = 1;
562  Matrix CHs;
563  Matrix CHs_check(c,c);
564  CHs_check(1,1) = 1.73205; CHs_check(1,2) = 0.00000; CHs_check(1,3) = 0.00000;
565  CHs_check(2,1) = 1.15470; CHs_check(2,2) = 0.81650; CHs_check(2,3) = 0.00000;
566  CHs_check(3,1) = 0.57735; CHs_check(3,2) = 0.40825; CHs_check(3,3) = 0.70711;
567  Ps.cholesky_semidefinite(CHs);
568  CPPUNIT_ASSERT_EQUAL(approxEqual(CHs, CHs_check, epsilon),true);
569 
570  // test operator - -=
571  Matrix Rm_bak; Rm_bak = Rm;
572  Matrix Dm(c,c); Dm = v;
573  Matrix Em = Rm - Dm;
574  Matrix Fm = Rm - v;
575  Matrix Gm(Rm);
576  Gm -= Dm;
577  Matrix Hm; Hm = Rm;
578  Hm -= v;
579  CPPUNIT_ASSERT_EQUAL(Rm, Rm_bak);
580  CPPUNIT_ASSERT_EQUAL( Em, Fm);
581  CPPUNIT_ASSERT_EQUAL( Fm, Gm);
582  CPPUNIT_ASSERT_EQUAL( Gm, Hm);
583  for (unsigned int i=0; i<c; i++)
584  for (unsigned int j=0; j<c; j++)
585  CPPUNIT_ASSERT_EQUAL( Dm(i+1,j+1), v);
586 
587  // test pseudoinverse
588  Matrix Bm_bak; Bm_bak = Bm;
589  Matrix Bm_pinv = Bm.pseudoinverse();
590  CPPUNIT_ASSERT_EQUAL(Bm, Bm_bak);
591  Im_test = Bm_pinv * Bm;
592  CPPUNIT_ASSERT_EQUAL(approxEqual(Im_test, Im, epsilon),true);
593 
594  // test svd
595  int rows = 4;
596  int cols = 3;
597  Matrix A_svd(rows,cols);
598  Matrix W_svd(cols,cols);
599  Matrix U_svd,V_svd;
600  W_svd = 0.0;
601 
602  ColumnVector w_svd;
603  A_svd(1,1)=1; A_svd(2,2)=2; A_svd(3,3)=3;
604  A_svd(1,2)=-0.5; A_svd(1,3)=-0.8;
605  A_svd(2,1)=-1.5; A_svd(2,3)=-2.8;
606  A_svd(3,1)=2.5; A_svd(3,2)=0.8;
607  A_svd(4,1)=0.5; A_svd(4,2)=1.8; A_svd(4,3)=1.6 ;
608 
609  A_svd.SVD(w_svd,U_svd,V_svd);
610  for (int i=1; i<=A_svd.columns() ; i++) W_svd(i,i) = w_svd(i);
611  CPPUNIT_ASSERT_EQUAL(approxEqual(A_svd, U_svd * W_svd * V_svd.transpose(), epsilon),true);
612 
613  int rows2 = 3;
614  int cols2 = 4;
615  Matrix A2_svd(rows2,cols2);
616  Matrix W2_svd(cols2,cols2);
617  Matrix U2_svd,V2_svd;
618  W2_svd = 0.0;
619 
620  ColumnVector w2_svd;
621  A2_svd(1,1)=1; A2_svd(2,2)=2; A2_svd(3,3)=3; //A(4,4)=4;
622  A2_svd(1,2)=-0.5; A2_svd(1,3)=-0.8; A2_svd(1,4)=-0.1 ;
623  A2_svd(2,1)=-1.5; A2_svd(2,3)=-2.8; A2_svd(2,4)=3.1 ;
624  A2_svd(3,1)=2.5; A2_svd(3,2)=-0.8; A2_svd(3,4)=1.1 ;
625 
626  A2_svd.SVD(w2_svd,U2_svd,V2_svd);
627  for (int i=1; i<=A2_svd.columns() ; i++) W2_svd(i,i) = w2_svd(i);
628  CPPUNIT_ASSERT_EQUAL(approxEqual(A2_svd, U2_svd * W2_svd * V2_svd.transpose(), epsilon),true);
629 
630  // TEST SPECIAL CASES
631  // Inverse for 1x1 Matrix
632  Matrix M1(1,1);
633  M1(1,1)= 1.4;
634  Matrix M1_inv = M1.inverse();
635  Matrix I1(1,1);
636  I1(1,1)= 1.0;
637  CPPUNIT_ASSERT_EQUAL(M1_inv * M1, I1);
638  // Inverse for 2x2 Matrix
639  Matrix M2(2,2);
640  M2(1,1)= 1.4;
641  M2(2,2)= 0.4;
642  M2(1,2)= 2.1;
643  M2(2,1)= -0.8;
644  Matrix M2_inv = M2.inverse();
645  Matrix I2(2,2);
646  I2=0.0;
647  I2(1,1)= 1.0;
648  I2(2,2)= 1.0;
649  CPPUNIT_ASSERT_EQUAL(approxEqual(M2_inv * M2, I2,epsilon),true);
650  // Determinant for 1x1 Matrix
651  CPPUNIT_ASSERT_EQUAL(M1.determinant(), M1(1,1));
652  // Determinant for 2x2 Matrix
653  CPPUNIT_ASSERT_EQUAL(M2.determinant(), M2(1,1)*M2(2,2)-M2(1,2)*M2(2,1));
654  // Inverse for 1x1 SymmetricMatrix
655  SymmetricMatrix SM1(1);
656  SM1(1,1)= 1.4;
657  SymmetricMatrix SM1_inv = SM1.inverse();
658  CPPUNIT_ASSERT_EQUAL(Matrix(SM1_inv * SM1), I1);
659  // Inverse for 2x2 Matrix
660  SymmetricMatrix SM2(2);
661  SM2(1,1)= 1.4;
662  SM2(2,2)= 0.4;
663  SM2(1,2)= 2.1;
664  SM2(2,1)= -0.8;
665  SymmetricMatrix SM2_inv = SM2.inverse();
666  CPPUNIT_ASSERT_EQUAL(approxEqual(Matrix(SM2_inv * SM2), I2,epsilon),true);
667  // Determinant for 1x1 Matrix
668  CPPUNIT_ASSERT_EQUAL(SM1.determinant(), SM1(1,1));
669  // Determinant for 2x2 Matrix
670  CPPUNIT_ASSERT_EQUAL(SM2.determinant(), SM2(1,1)*SM2(2,2)-SM2(1,2)*SM2(2,1));
671 
672  Matrix M3(3,3);
673  M3(1,1)=1;
674  M3(1,2)=2;
675  M3(1,3)=3;
676  M3(2,1)=4;
677  M3(2,2)=5;
678  M3(2,3)=6;
679  M3(3,1)=7;
680  M3(3,2)=8;
681  M3(3,3)=9;
682  // test rowCopy()
683  RowVector rcopy1 = M3.rowCopy(1);
684  RowVector rcopy1test(3);
685  rcopy1test(1) = M3(1,1);
686  rcopy1test(2) = M3(1,2);
687  rcopy1test(3) = M3(1,3);
688  CPPUNIT_ASSERT_EQUAL(rcopy1,rcopy1test);
689  RowVector rcopy2 = M3.rowCopy(2);
690  RowVector rcopy2test(3);
691  rcopy2test(1) = M3(2,1);
692  rcopy2test(2) = M3(2,2);
693  rcopy2test(3) = M3(2,3);
694  CPPUNIT_ASSERT_EQUAL(rcopy2,rcopy2test);
695  RowVector rcopy3 = M3.rowCopy(3);
696  RowVector rcopy3test(3);
697  rcopy3test(1) = M3(3,1);
698  rcopy3test(2) = M3(3,2);
699  rcopy3test(3) = M3(3,3);
700  CPPUNIT_ASSERT_EQUAL(rcopy3,rcopy3test);
701 
702  // test columnCopy()
703  ColumnVector copy1 = M3.columnCopy(1);
704  ColumnVector copy1test(3);
705  copy1test(1) = M3(1,1);
706  copy1test(2) = M3(2,1);
707  copy1test(3) = M3(3,1);
708  CPPUNIT_ASSERT_EQUAL(copy1,copy1test);
709  ColumnVector copy2 = M3.columnCopy(2);
710  ColumnVector copy2test(3);
711  copy2test(1) = M3(1,2);
712  copy2test(2) = M3(2,2);
713  copy2test(3) = M3(3,2);
714  CPPUNIT_ASSERT_EQUAL(copy2,copy2test);
715  ColumnVector copy3 = M3.columnCopy(3);
716  ColumnVector copy3test(3);
717  copy3test(1) = M3(1,3);
718  copy3test(2) = M3(2,3);
719  copy3test(3) = M3(3,3);
720  CPPUNIT_ASSERT_EQUAL(copy3,copy3test);
721 }
722 
723 
724 
double epsilon
CPPUNIT_TEST_SUITE_REGISTRATION(MatrixwrapperTest)
bool approxEqual(double a, double b, double epsilon)
Definition: approxEqual.cpp:20


bfl
Author(s): Klaas Gadeyne, Wim Meeussen, Tinne Delaet and many others. See web page for a full contributor list. ROS package maintained by Wim Meeussen.
autogenerated on Mon Jun 10 2019 12:47:59