Justine Devin / gammalib

Commit b4ded5bf079bf7b2c29e4bec9ef7deed4e6bc10a

Authored by Jürgen Knödlseder
1 parent a4f1132c

Improve the error messages. 

Improve some error messages so that more information is available in the report in case of failures.
Inline Side-by-side
Showing 3 changed files with 733 additions and 671 deletions
test/test_GMatrix.cpp
  View file @b4ded5b
1 1 /***************************************************************************
2   - * test_GMatrix.cpp - test GMatrix class *
  2 + * test_GMatrix.cpp - Test generic matrix class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2006-2011 by Jurgen Knodlseder *
  4 + * copyright (C) 2006-2012 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -18,778 +18,834 @@
18 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19 19 * *
20 20 ***************************************************************************/
  21 +/**
  22 + * @file test_GMatrix.cpp
  23 + * @brief Implementation of unit tests for generic matrices
  24 + * @author Juergen Knoedlseder
  25 + */
21 26  
22 27 /* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
23 31 #include <cmath>
24   -#include <iostream> // std::cout, cerr
25   -#include <stdexcept> // std::exception
26   -#include "GammaLib.hpp"
27 32 #include "test_GMatrix.hpp"
28   -/* __ Namespaces _________________________________________________________ */
29 33  
30 34 /* __ Globals ____________________________________________________________ */
  35 +double g_matrix[] = {1.0, 2.0, 3.0, 4.0,
  36 + 5.0, 6.0, 7.0, 8.0,
  37 + 9.0, 10., 11., 12.};
  38 +double g_vector[] = {1.0, 2.0, 3.0, 4.0};
  39 +int g_rows = 3;
  40 +int g_cols = 4;
31 41  
32   -void TestGMatrix::set(void){
33   - // Test name
34   - name("GMatrix");
35   -
36   - //Initialize m_test
37   - init_matrix();
38   -
39   - //Initialize m_vector
40   - init_vector();
41   -
42   - //Set parameters
43   - m_g_rows = 3;
44   - m_g_cols = 4;
45   -
46   - //add tests
47   - add_test(static_cast<pfunction>(&TestGMatrix::set_matrix),"Set matrix");
48   - add_test(static_cast<pfunction>(&TestGMatrix::set_vector),"Set Vector");
49   - add_test(static_cast<pfunction>(&TestGMatrix::set_bigger),"Set matrix bigger");
50   - add_test(static_cast<pfunction>(&TestGMatrix::test_output),"Output test matrix");
51   - add_test(static_cast<pfunction>(&TestGMatrix::test_conversion),"Matrix conversions");
52   - add_test(static_cast<pfunction>(&TestGMatrix::test_extract),"Vector extraction, insertion and addition");
53   -
54   - add_test(static_cast<pfunction>(&TestGMatrix::test1),"Test 1: Allocate zero matrix");
55   - //add_test(static_cast<pfunction>(&TestGMatrix::test2),"Test 2: Allocate too large matrix");
56   - add_test(static_cast<pfunction>(&TestGMatrix::test3),"Test3: Assign values");
57   - add_test(static_cast<pfunction>(&TestGMatrix::test4),"Test 4: Matrix copy constructor");
58   - add_test(static_cast<pfunction>(&TestGMatrix::test5),"Test 5: Matrix assignment");
59   - add_test(static_cast<pfunction>(&TestGMatrix::test6),"Test 6: Transposition");
60   - add_test(static_cast<pfunction>(&TestGMatrix::test7),"Test 7: Matrix*Vector multiplication");
61   - add_test(static_cast<pfunction>(&TestGMatrix::test8),"Test 8: Matrix*Matrix multiplication");
62   - add_test(static_cast<pfunction>(&TestGMatrix::test9),"Test 9: Assignment and arithmetics");
63   - add_test(static_cast<pfunction>(&TestGMatrix::test10),"Test 10: Matrix functions");
64   - add_test(static_cast<pfunction>(&TestGMatrix::test11),"Test 11: Matrix comparison");
65   -
66   - return;
67   -}
68 42  
69 43 /***************************************************************************
70   -* Initializer m_test *
71   -***************************************************************************/
72   -void TestGMatrix::init_matrix(void)
  44 + * @brief Set test matrix
  45 + ***************************************************************************/
  46 +GMatrix TestGMatrix::set_matrix(void) const
73 47 {
74   - double g_matrix[] = {1.0, 2.0, 3.0, 4.0,
75   - 5.0, 6.0, 7.0, 8.0,
76   - 9.0, 10., 11., 12.};
  48 + // Allocate matrix
  49 + GMatrix matrix(g_rows,g_cols);
77 50  
78   - for(int i=0;i<12;i++){
79   - m_g_matrix[i]=g_matrix[i];
  51 + // Set matrix values
  52 + for (int row = 0; row < g_rows; ++row) {
  53 + for (int col = 0; col < g_cols; ++col) {
  54 + matrix(row,col) = g_matrix[col+row*g_cols];
  55 + }
80 56 }
81 57  
82   - return;
  58 + // Return matrix
  59 + return matrix;
83 60 }
84 61  
85   -/***************************************************************************
86   -* Initializer m_vector *
87   -***************************************************************************/
88   -void TestGMatrix::init_vector(void)
89   -{
90   - double g_vector[] = {1.0, 2.0, 3.0, 4.0};
91   -
92   - for(int i=0;i<4;i++){
93   - m_g_vector[i]=g_vector[i];
94   - }
95   -
96   - return;
97   -}
98 62  
99 63 /***************************************************************************
100   - * Set test matrix *
  64 + * @brief Set test matrix with zero lines
101 65 ***************************************************************************/
102   -void TestGMatrix::set_matrix(void)
  66 +GMatrix TestGMatrix::set_matrix_zero(void) const
103 67 {
104   - GMatrix matrix(m_g_rows,m_g_cols);
105   - for (int row = 0; row < m_g_rows; ++row) {
106   - for (int col = 0; col < m_g_cols; ++col)
107   - matrix(row,col) = m_g_matrix[col+row*m_g_cols];
  68 + // Allocate matrix
  69 + GMatrix matrix(g_rows+1,g_cols+1);
  70 +
  71 + // Set matrix values
  72 + int i = 0;
  73 + for (int row = 0; row < g_rows+1; ++row) {
  74 + if (row == 2) continue;
  75 + int j = 0;
  76 + for (int col = 0; col < g_cols+1; ++col) {
  77 + if (col == 2) continue;
  78 + matrix(row,col) = g_matrix[j+i*g_cols];
  79 + j++;
  80 + }
  81 + i++;
108 82 }
109 83  
110   - m_test = matrix;
111   -
112   - return;
  84 + // Return matrix
  85 + return matrix;
113 86 }
114 87  
115 88  
116 89 /***************************************************************************
117   - * Set test vector *
  90 + * @brief Set test vector
118 91 ***************************************************************************/
119   -void TestGMatrix::set_vector(void)
  92 +GVector TestGMatrix::set_vector(void) const
120 93 {
121   - GVector vector(m_g_cols);
122   - for (int col = 0; col < m_g_cols; ++col){
123   - vector[col] = m_g_vector[col];
124   - }
  94 + // Allocate vector
  95 + GVector vector(g_cols);
125 96  
126   - v_test = vector;
  97 + // Set vector values
  98 + for (int col = 0; col < g_cols; ++col) {
  99 + vector[col] = g_vector[col];
  100 + }
127 101  
128   - return;
  102 + // Return vector
  103 + return vector;
129 104 }
130 105  
131   -void TestGMatrix::set_bigger(void)
132   -{
133   - bigger= GMatrix(m_g_rows+1,m_g_cols+1);
134   -}
  106 +
135 107 /***************************************************************************
136   - * Check test matrix *
  108 + * @brief Check if matrix corresponds to test matrix
  109 + *
  110 + * @param[in] matrix Matrix.
  111 + * @param[in] scale scale factor (default: 1.0).
  112 + * @param[in] offset offset value (default: 0.0).
  113 + *
  114 + * Checks if a matrix corresponds to the test matrix. Optionally, the test
  115 + * matrix can be scaled and an offset can be added.
137 116 ***************************************************************************/
138   -int TestGMatrix::check_matrix(const GMatrix& m, const double scale, const double add)
  117 +bool TestGMatrix::check_matrix(const GMatrix& matrix,
  118 + const double& scale,
  119 + const double& offset) const
139 120 {
140   - int result = 1;
141   -
142   - for (int row = 0; row < m_g_rows; ++row) {
143   - for (int col = 0; col < m_g_cols; ++col) {
144   - double value = m_g_matrix[col+row*m_g_cols] * scale + add;
145   - if (fabs(m(row,col)-value) > 1.0e-15) {
146   - result = 0;
  121 + // Initialise check with true
  122 + bool result = true;
  123 +
  124 + // Compare all matrix elements
  125 + for (int row = 0; row < g_rows; ++row) {
  126 + for (int col = 0; col < g_cols; ++col) {
  127 + double value = g_matrix[col+row*g_cols] * scale + offset;
  128 + if (std::abs(matrix(row,col)-value) > 1.0e-15) {
  129 + result = false;
147 130 break;
148   - }
149   - }
  131 + }
150 132 }
  133 + }
151 134  
  135 + // Return result
152 136 return result;
153 137 }
154 138  
155 139  
156 140 /***************************************************************************
157   - * Check transposed test matrix *
  141 + * @brief Check if transposed matrix corresponds to test matrix
  142 + *
  143 + * @param[in] matrix Trasposed matrix.
  144 + * @param[in] scale scale factor (default: 1.0).
  145 + * @param[in] offset offset value (default: 0.0).
  146 + *
  147 + * Checks if a matrix corresponds to the test matrix. Optionally, the test
  148 + * matrix can be scaled and an offset can be added.
158 149 ***************************************************************************/
159   -int TestGMatrix::check_transpose_matrix(const GMatrix& m, const double scale, const double add)
  150 +bool TestGMatrix::check_matrix_trans(const GMatrix& matrix,
  151 + const double& scale,
  152 + const double& offset) const
160 153 {
161   - int result = 1;
162   -
163   - for (int row = 0; row < m_g_rows; ++row) {
164   - for (int col = 0; col < m_g_cols; ++col) {
165   - double value = m_g_matrix[col+row*m_g_cols] * scale + add;
166   - if (m(col,row) != value) {
167   - result = 0;
  154 + // Initialise check with true
  155 + bool result = true;
  156 +
  157 + // Compare all matrix elements
  158 + for (int row = 0; row < g_rows; ++row) {
  159 + for (int col = 0; col < g_cols; ++col) {
  160 + double value = g_matrix[col+row*g_cols] * scale + offset;
  161 + if (std::abs(matrix(col,row)-value) > 1.0e-15) {
  162 + result = false;
168 163 break;
169   - }
170   - }
  164 + }
171 165 }
  166 + }
172 167  
  168 + // Return result
173 169 return result;
174 170 }
175 171  
176 172  
177 173 /***************************************************************************
178   - * Check full test matrix (lower triangle) *
  174 + * @brief Check lower triangle general matrix
  175 + *
  176 + * @param[in] matrix Matrix.
  177 + * @param[in] ref Reference matrix.
179 178 ***************************************************************************/
180   -int TestGMatrix::check_matrix_lt(const GMatrix& m, const double scale, const double add)
  179 +bool TestGMatrix::check_matrix_lt(const GMatrix& matrix,
  180 + const GMatrix& ref) const
181 181 {
182   - int result = 1;
183   - test_try("Check full test matrix (lower triangle)");
184   - try {
185   - for (int row = 0; row < m_g_rows; ++row) {
186   - for (int col = 0; col < m_g_cols; ++col) {
187   - double value = (col <= row) ? m_g_matrix[col+row*m_g_cols] * scale + add : 0.0;
188   - if (m(row,col) != value) {
189   - result = 0;
190   - break;
191   - }
192   - }
  182 + // Initialise check with true
  183 + bool result = true;
  184 +
  185 + // Compare all matrix elements
  186 + for (int row = 0; row < ref.rows(); ++row) {
  187 + for (int col = 0; col < ref.cols(); ++col) {
  188 + double value = (col <= row) ? ref(row,col) : 0.0;
  189 + if (matrix(row,col) != value) {
  190 + result = false;
  191 + break;
193 192 }
194   - test_try_success();
195   - }
196   - catch (std::exception &e) {
197   - test_try_failure(e);
198   - }
199   - return result;
200   -}
201   -
  193 + }
  194 + }
202 195  
203   -/***************************************************************************
204   - * Check full test matrix (upper triangle) *
205   - ***************************************************************************/
206   -int TestGMatrix::check_matrix_ut(const GMatrix& m, const double scale, const double add)
207   -{
208   - int result = 1;
209   - test_try("Check full test matrix (upper triangle)");
210   - try {
211   - for (int row = 0; row < m_g_rows; ++row) {
212   - for (int col = 0; col < m_g_cols; ++col) {
213   - double value = (col >= row) ? m_g_matrix[col+row*m_g_cols] * scale + add : 0.0;
214   - if (m(row,col) != value) {
215   - result = 0;
216   - break;
217   - }
218   - }
219   - }
220   - test_try_success();
221   - }
222   - catch (std::exception &e) {
223   - test_try_failure(e);
224   - }
  196 + // Return result
225 197 return result;
226 198 }
227 199  
228 200  
229 201 /***************************************************************************
230   - * Check matrix*vector *
  202 + * @brief Check upper triangle general matrix
  203 + *
  204 + * @param[in] matrix Matrix.
  205 + * @param[in] ref Reference matrix.
231 206 ***************************************************************************/
232   -int TestGMatrix::check_matrix_vector(const GVector& v)
  207 +bool TestGMatrix::check_matrix_ut(const GMatrix& matrix,
  208 + const GMatrix& ref) const
233 209 {
234   - int result = 1;
235   - test_try(" Check matrix*vector");
236   - try {
237   - for (int row = 0; row < v.size(); ++row) {
238   - double value = 0.0;
239   - for (int col = 0; col < m_g_cols; ++col)
240   - value += m_g_matrix[col+row*m_g_cols] * m_g_vector[col];
241   - if (v[row] != value) {
242   - result = 0;
  210 + // Initialise check with true
  211 + bool result = true;
  212 +
  213 + // Compare all matrix elements
  214 + for (int row = 0; row < ref.rows(); ++row) {
  215 + for (int col = 0; col < ref.cols(); ++col) {
  216 + double value = (col >= row) ? ref(row,col) : 0.0;
  217 + if (matrix(row,col) != value) {
  218 + result = false;
243 219 break;
244 220 }
245 221 }
246   - test_try_success();
247   - }
248   - catch (std::exception &e) {
249   - test_try_failure(e);
250   - }
251   - return result;
252   -}
253   -
254   -
255   -/***************************************************************************
256   - * Check matrix*matrix *
257   - ***************************************************************************/
258   -int TestGMatrix::check_matrix_matrix(const GMatrix& m)
259   -{
260   - if (m.rows() != m_g_rows || m.cols() != m_g_rows)
261   - return 0;
262   - int result = 1;
263   - test_try("Check matrix*matrix");
264   - try {
265   - for (int row = 0; row < m.rows(); ++row) {
266   - for (int col = 0; col < m.cols(); ++col) {
267   - double value = 0.0;
268   - for (int i = 0; i < m_g_cols; ++i)
269   - value += m_g_matrix[i+row*m_g_cols] * m_g_matrix[i+col*m_g_cols];
270   - if (m(row,col) != value) {
271   - result = 0;
272   - break;
273   - }
274   - }
275 222 }
276   - test_try_success();
277   - }
278   - catch (std::exception &e) {
279   - test_try_failure(e);
280   - }
281   - return result;
282   -}
283   -
284   -
285   -/***************************************************************************
286   - * Check matrix min *
287   - ***************************************************************************/
288   -int TestGMatrix::check_matrix_min(const double min)
289   -{
290   - double value = m_g_matrix[0];
291   - for (int row = 0; row < m_g_rows; ++row) {
292   - for (int col = 0; col < m_g_cols; ++col) {
293   - if (m_g_matrix[col+row*m_g_cols] < value)
294   - value = m_g_matrix[col+row*m_g_cols];
295   - }
296   - }
297   - test_assert(min==value, "Check matrix min");
298   - return (min == value);
299   -}
300 223  
301   -
302   -/***************************************************************************
303   - * Check matrix max *
304   - ***************************************************************************/
305   -int TestGMatrix::check_matrix_max(const double max)
306   -{
307   - double value = m_g_matrix[0];
308   - for (int row = 0; row < m_g_rows; ++row) {
309   - for (int col = 0; col < m_g_cols; ++col) {
310   - if (m_g_matrix[col+row*m_g_cols] > value)
311   - value = m_g_matrix[col+row*m_g_cols];
312   - }
313   - }
314   - test_assert(max==value,"Check matrix max");
315   - return (max == value);
  224 + // Return result
  225 + return result;
316 226 }
317 227  
318 228  
319   -/***************************************************************************
320   - * Check matrix sum *
321   - ***************************************************************************/
322   -int TestGMatrix::check_matrix_sum(const double sum)
  229 +/***********************************************************************//**
  230 + * @brief Set parameters and tests
  231 + **************************************************************************/
  232 +void TestGMatrix::set(void)
323 233 {
324   - double value = 0.0;
325   - for (int row = 0; row < m_g_rows; ++row) {
326   - for (int col = 0; col < m_g_cols; ++col)
327   - value += m_g_matrix[col+row*m_g_cols];
328   - }
329   - test_assert(sum==value,"Check matrix sum");
330   - return (sum == value);
331   -}
  234 + // Set test name
  235 + name("GMatrix");
332 236  
  237 + // Add tests
  238 + append(static_cast<pfunction>(&TestGMatrix::alloc_matrix), "Test matrix allocation");
  239 + append(static_cast<pfunction>(&TestGMatrix::assign_values), "Test value assignment");
  240 + append(static_cast<pfunction>(&TestGMatrix::copy_matrix), "Test matrix copying");
  241 + append(static_cast<pfunction>(&TestGMatrix::matrix_operations), "Test matrix operations");
  242 + append(static_cast<pfunction>(&TestGMatrix::matrix_arithmetics), "Test matrix arithmetics");
  243 + append(static_cast<pfunction>(&TestGMatrix::matrix_functions), "Test matrix functions");
  244 + append(static_cast<pfunction>(&TestGMatrix::matrix_compare), "Test matrix comparisons");
  245 + //append(static_cast<pfunction>(&TestGMatrix::matrix_cholesky), "Test matrix Cholesky decomposition");
  246 + append(static_cast<pfunction>(&TestGMatrix::matrix_print), "Test matrix printing");
  247 +
  248 + // Set members
  249 + m_test = set_matrix();
  250 + v_test = set_vector();
  251 + m_bigger = GMatrix(g_rows+1, g_cols+1);
333 252  
334   -/***************************************************************************
335   - * Test: Output *
336   - ***************************************************************************/
337   -void TestGMatrix::test_output(void)
338   -{
339   - std::cout << "Output test matrix"<<std::endl<< m_test << std::endl;
  253 + // Return
  254 + return;
340 255 }
341 256  
342 257  
343   -/***************************************************************************
344   - * Test: Conversion between matrix types *
  258 +/***********************************************************************//**
  259 + * @brief Test matrix allocation
345 260 ***************************************************************************/
346   -void TestGMatrix::test_conversion(void)
  261 +void TestGMatrix::alloc_matrix(void)
347 262 {
348   - //
349   - // Setup a symmetric matrix
350   - int num = 10;
351   - GMatrix symmetric(num,num);
352   - for (int i = 0; i < num; ++i) {
353   - for (int j = 0; j < num; ++j)
354   - symmetric(i,j) = (i+j+1)/2.0;
355   - }
356   -
357   - //
358   - // Convert symmetric matrix into GSymMatrix object
359   - // GSymMatrix converted = sym_matrix(symmetric);
360   - GSymMatrix converted = GSymMatrix(symmetric);
361   -
362   - //
363   - // Convert GSymMatrix back to full matrix
364   - GMatrix back_convert = matrix(converted);
365   -
366   - //
367   - // Compare back converted matrix to original one. They should be identical
368   - test_assert(symmetric == back_convert,"Compare back converted matrix to original one");
369   - //
370   - // Determine the fill of the matrix. It should be 1.0
371   - double fill = back_convert.fill();
372   - test_assert(fabs(fill-1.0) <= 1.0e-15,"Test the fill of the matrix","Bad fill "+str(fill)+" determined (expected 1.0).");
373   - //
374   - // Extract lower triangle and check values
375   - GMatrix lower = symmetric.extract_lower_triangle();
376   - int ok = 1;
377   - for (int j = 1; j < num; ++j) {
378   - for (int i = 0; i < j; ++i) {
379   - if (lower(i,j) != 0.0)
380   - ok = 0;
381   - }
382   - }
383   - for (int j = 0; j < num; ++j) {
384   - for (int i = j; i < num; ++i) {
385   - if (lower(i,j) != symmetric(i,j))
386   - ok = 0;
387   - }
388   - }
389   - test_assert(ok,"Extract lower triangle");
390   -
391   - //
392   - // Extract upper triangle and check values
393   - GMatrix upper = symmetric.extract_upper_triangle();
394   - ok = 1;
395   - for (int j = 0; j < num; ++j) {
396   - for (int i = j+1; i < num; ++i) {
397   - if (upper(i,j) != 0.0)
398   - ok = 0;
399   - }
400   - }
401   - for (int j = 0; j < num; ++j) {
402   - for (int i = 0; i <= j; ++i) {
403   - if (upper(i,j) != symmetric(i,j))
404   - ok = 0;
405   - }
406   - }
407   -
408   - test_assert(ok,"Extract upper triangle","Corrupt extract_upper_triangle.");
409   -
410   - //
411   - // Now make the matrix unsymmetric
412   - symmetric(0,num-1) = 1000.0;
413   - //
414   - // Try converting now into GSymMatrix object (this should fail)
415   - test_try("Try converting now into GSymMatrix object");
  263 + // Allocate zero matrix. The allocation should fail.
  264 + test_try("Allocate zero matrix");
416 265 try {
417   - converted = sym_matrix(symmetric);
418   - test_try_failure();
  266 + GMatrix test(0,0);
  267 + test_try_failure("Expected GException::empty exception.");
419 268 }
420   - catch (GException::matrix_not_symmetric &e) {
  269 + catch (GException::empty &e) {
421 270 test_try_success();
422 271 }
423 272 catch (std::exception &e) {
424 273 test_try_failure(e);
425 274 }
426   - //
427   - // Now zero some elements to emulate a sparse matrix
428   - symmetric(0,num-1) = 0.0;
429   - for (int i = 5; i < 8; ++i) {
430   - for (int j = i-2; j < i+2; ++j)
431   - symmetric(i,j) = 0.0;
432   - }
433   -
434   - //
435   - // Convert symmetric matrix into GSparseMatrix object
436   - GSparseMatrix sparse = sparse_matrix(symmetric);
437   - //
438   - // Convert GSparseMatrix back to full matrix
439   - back_convert = matrix(sparse);
440   - //
441   - // Compare back converted matrix to original one. They should be identical
442   - test_assert(symmetric == back_convert,
443   - "Compare back converted matrix to original one",
444   - "Unable to convert matrixes (sparse)");
445 275  
446 276 // Return
447 277 return;
448 278 }
449 279  
450 280  
451   -/***************************************************************************
452   - * Test: extraction and insertion *
  281 +/***********************************************************************//**
  282 + * @brief Test value assignment
453 283 ***************************************************************************/
454   -void TestGMatrix::test_extract(void)
  284 +void TestGMatrix::assign_values(void)
455 285 {
456   -
457   - //
458   - // Set-up test matrix
459   - int rows = 10;
460   - int cols = 20;
461   - GMatrix test(rows, cols);
462   -
463   - //
464   - // Add and extract column vectors
465   - test_try("Add and extract column vectors");
466   - try
467   - {
468   - for (int col = 0; col < cols; ++col) {
469   - GVector column(rows);
470   - for (int row = 0; row < rows; ++row)
471   - column[row] = (col+1)*100.0 + (row+1)*1.0;
472   - test.add_col(column, col);
473   - GVector check = test.extract_col(col);
474   - if (check != column) {
475   - throw exception_failure("Unable to add and extract columns.\nAdded column ...: "+column.print()+"\nExtracted column: "+check.print());
476   - }
  286 + // Setup 3x3 matrix
  287 + GMatrix test(3,3);
  288 +
  289 + // Assignment individual values
  290 + for (int i = 0; i < 3; ++i) {
  291 + for (int k = 0; k < 3; ++k) {
  292 + double value = i*2.0 + k*2.0;
  293 + test(i,k) = value;
477 294 }
478   - test_try_success();
479   - }
480   - catch(std::exception& e)
481   - {
482   - test_try_failure(e);
483 295 }
484 296  
485   - //
486   - // Insert and extract column vectors
487   - test_try("Insert and extract column vectors");
488   - try
489   - {
490   - for (int col = 0; col < cols; ++col) {
491   - GVector column(rows);
492   - for (int row = 0; row < rows; ++row)
493   - column[row] = (col+1)*100.0 + (row+1)*1.0;
494   - test.insert_col(column, col);
495   - GVector check = test.extract_col(col);
496   - if (check != column) {
497   - throw exception_failure("Unable to insert and extract columns.\nIserted column ...: "+column.print()+"\nExtracted column: "+check.print());
498   - }
  297 + // Check assignment of individual values
  298 + for (int i = 0; i < 3; ++i) {
  299 + for (int k = 0; k < 3; ++k) {
  300 + double value = i*2.0 + k*2.0;
  301 + test_value(test(i,k), value, 1.0e-10, "Test matrix element assignment");
499 302 }
500   - test_try_success();
501   - }
502   - catch(std::exception& e)
503   - {
504   - test_try_failure(e);
505 303 }
506 304  
507   - //
508   - // Extract rows
509   - test_try("Extract rows");
510   - try
511   - {
512   - for (int row = 0; row < rows; ++row) {
513   - GVector v_row(cols);
514   - for (int col = 0; col < cols; ++col)
515   - v_row[col] = (col+1)*100.0 + (row+1)*1.0;
516   - GVector check = test.extract_row(row);
517   - if (check != v_row) {
518   - throw exception_failure("Unable to extract rows.\nIserted row ...: "+v_row.print()+"\nExtracted row: "+check.print());
519   - }
520   - }
521   -
  305 + // Verify range checking
  306 + #ifdef G_RANGE_CHECK
  307 + test_try("Verify range checking");
  308 + try {
  309 + test(3,3) = 1.0;
  310 + test_try_failure("Expected GException::out_of_range exception.");
  311 + }
  312 + catch (GException::out_of_range &e) {
522 313 test_try_success();
523 314 }
524   - catch(std::exception& e)
525   - {
  315 + catch (std::exception &e) {
526 316 test_try_failure(e);
527 317 }
  318 + #endif
528 319  
529   -}
530   -
531   -// Test 1: Allocate zero matrix
532   -void TestGMatrix::test1(void)
533   -{
534   - GMatrix test1(0,0);
535   -
  320 + // Return
536 321 return;
537 322 }
538 323  
539   -// Test 2: Allocate too large matrix
540   -void TestGMatrix::test2(void)
  324 +
  325 +/***********************************************************************//**
  326 + * @brief Test matrix copy
  327 + ***************************************************************************/
  328 +void TestGMatrix::copy_matrix(void)
541 329 {
542   - //GMatrix test2(100000,100000);
  330 + // Copy matrix
  331 + GMatrix test = m_test;
  332 +
  333 + // Test if original and compied matrices are correct
  334 + test_assert(check_matrix(m_test), "Test source matrix");
  335 + test_assert(check_matrix(test), "Test matrix copy operator",
  336 + "Found:\n"+test.print()+"\nExpected:\n"+m_test.print());
543 337  
  338 + // Return
544 339 return;
545 340 }
546 341  
547   -// Test 3: Assign values
548   -void TestGMatrix::test3(void)
  342 +
  343 +/***********************************************************************//**
  344 + * @brief Test matrix operations
  345 + *
  346 + * Tests matrix*vector and matrix*matrix multiplication operations.
  347 + ***************************************************************************/
  348 +void TestGMatrix::matrix_operations(void)
549 349 {
550   - GMatrix test3(3,3);
551   - test3(1,1) = 1.0;
552   - test_assert(test3(0,0) == 0.0 && test3(1,0) == 0.0 && test3(2,0) == 0.0 &&
553   - test3(0,1) == 0.0 && test3(1,1) == 1.0 && test3(2,1) == 0.0 &&
554   - test3(0,2) == 0.0 && test3(1,2) == 0.0 && test3(2,2) == 0.0,"Assign value");
555   -
556   - #if defined(G_RANGE_CHECK)
557   - test_try("Range check");
558   - try {
559   - GMatrix test3(3,3);
560   - test3(3,3) = 1.0;
561   - test_try_failure();
562   - }
563   - catch (GException::out_of_range &e) {
564   - test_try_success();
  350 + // Perform vector multiplication
  351 + GVector test1 = m_test * v_test;
  352 +
  353 + // Check if the result vector is as expected
  354 + GVector ref1 = test1;
  355 + bool result = true;
  356 + for (int row = 0; row < g_rows; ++row) {
  357 + double value = 0.0;
  358 + for (int col = 0; col < g_cols; ++col) {
  359 + value += g_matrix[col+row*g_cols] * g_vector[col];
565 360 }
566   - catch (std::exception &e) {
567   - test_try_failure(e);
  361 + ref1[row] = value;
  362 + if (test1[row] != value) {
  363 + result = false;
  364 + break;
568 365 }
569   - #endif
  366 + }
570 367  
571   - return;
572   -}
  368 + // Test if original matrix and result vector are correct
  369 + test_assert(check_matrix(m_test), "Test source matrix");
  370 + test_assert(result, "Test matrix*vector multiplication",
  371 + "Found:\n"+test1.print()+"\nExpected:\n"+ref1.print());
573 372  
574   -// Test 4: Matrix copy constructor
575   -void TestGMatrix::test4(void)
576   -{
577   - GMatrix test4 = m_test;
578   - if (!check_matrix(test4, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {
579   - throw exception_failure("Corrupt copy constructor");
  373 + // Test incompatible vector multiplication
  374 + test_try("Test incompatible matrix*vector multiplication");
  375 + try {
  376 + GVector test2 = m_bigger * v_test;
  377 + test_try_failure("Expected GException::matrix_vector_mismatch exception.");
  378 + }
  379 + catch (GException::matrix_vector_mismatch &e) {
  380 + test_try_success();
  381 + }
  382 + catch (std::exception &e) {
  383 + test_try_failure(e);
  384 + }
  385 +
  386 + // Test matrix multiplication
  387 + GMatrix test3 = m_test * transpose(m_test);
  388 +
  389 + // Check if the result matrix is as expected
  390 + GMatrix ref3 = test3;
  391 + result = true;
  392 + for (int row = 0; row < test3.rows(); ++row) {
  393 + for (int col = 0; col < test3.cols(); ++col) {
  394 + double value = 0.0;
  395 + for (int i = 0; i < g_cols; ++i) {
  396 + value += g_matrix[i+row*g_cols] * g_matrix[i+col*g_cols];
  397 + }
  398 + ref3(row,col) = value;
  399 + if (test3(row,col) != value) {
  400 + result = false;
  401 + break;
  402 + }
  403 + }
580 404 }
581 405  
582   - return;
583   -}
  406 + // Test if original matrix and result matrix are correct
  407 + test_assert(check_matrix(m_test), "Test source matrix");
  408 + test_assert(result, "Test matrix multiplication",
  409 + "Found:\n"+test3.print()+"\nExpected:\n"+ref3.print());
  410 + test_value(test3.rows(), g_rows, "Test number of rows of result matrix");
  411 + test_value(test3.cols(), g_rows, "Test number of columns of result matrix");
584 412  
585   -// Test 5: Matrix assignment
586   -void TestGMatrix::test5(void)
587   -{
588   - //
589   - // GMatrix = GMatrix
590   - result = m_test;
591   - if (!check_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {
592   - throw exception_failure("Corrupt assignment");
  413 + // Test incompatible matrix multiplication
  414 + test_try("Test incompatible matrix multiplication");
  415 + try {
  416 + GMatrix test4 = m_bigger * m_test;
  417 + test_try_failure("Expected GException::matrix_mismatch exception.");
  418 + }
  419 + catch (GException::matrix_mismatch &e) {
  420 + test_try_success();
  421 + }
  422 + catch (std::exception &e) {
  423 + test_try_failure(e);
593 424 }
594   - //
595   - // GMatrix = GMatrix (bigger matrix)
596   - result = bigger;
597 425  
  426 + // Test another incompatible matrix multiplication
  427 + test_try("Test incompatible matrix multiplication");
  428 + try {
  429 + GMatrix test5 = m_bigger * m_test;
  430 + test_try_failure("Expected GException::matrix_mismatch exception.");
  431 + }
  432 + catch (GException::matrix_mismatch &e) {
  433 + test_try_success();
  434 + }
  435 + catch (std::exception &e) {
  436 + test_try_failure(e);
  437 + }
  438 +
  439 + // Return
598 440 return;
599 441 }
600 442  
601   -// Test 6: Transposition
602   -void TestGMatrix::test6(void)
603   -{
604   - //
605   - // transpose(GMatrix)
606   - result = transpose(m_test);
607   - test_assert(check_transpose_matrix(result, 1.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"transpose(GMatrix)","Corrupt transpose(GMatrix) function.");
608   - //
609   - // GMatrix.transpose()
610   - result = m_test;
611   - result.transpose();
612   - test_assert(check_transpose_matrix(result, 1.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix.transpose()","Corrupt GMatrix.transpose() function.");
613   -}
614 443  
615   -// Test 7: Matrix*Vector multiplication
616   -void TestGMatrix::test7(void)
  444 +/***********************************************************************//**
  445 + * @brief Test matrix arithmetics
  446 + *
  447 + * Tests matrix arithmetics.
  448 + ***************************************************************************/
  449 +void TestGMatrix::matrix_arithmetics(void)
617 450 {
618   - GVector v_test7 = m_test*v_test;
619   - test_assert(check_matrix_vector(v_test7) && check_matrix(m_test, 1.0, 0.0),"","Corrupt Matrix*Vector multiplication");
620   -
621   - test_try("Bigger*Vector");
  451 + // -GMatrix
  452 + GMatrix test = -m_test;
  453 + test_assert(check_matrix(m_test), "Test source matrix");
  454 + test_assert(check_matrix(test, -1.0, 0.0), "Test -GMatrix",
  455 + "Unexpected result matrix:\n"+test.print());
  456 +
  457 + // GMatrix += GMatrix
  458 + test = m_test;
  459 + test += m_test;
  460 + test_assert(check_matrix(m_test), "Test source matrix");
  461 + test_assert(check_matrix(test, 2.0, 0.0), "Test GMatrix += GMatrix",
  462 + "Unexpected result matrix:\n"+test.print());
  463 +
  464 + // GMatrix -= GMatrix
  465 + test = m_test;
  466 + test -= m_test;
  467 + test_assert(check_matrix(m_test), "Test source matrix");
  468 + test_assert(check_matrix(test, 0.0, 0.0), "Test GMatrix -= GMatrix",
  469 + "Unexpected result matrix:\n"+test.print());
  470 +
  471 + // GMatrix *= 3.0
  472 + test = m_test;
  473 + test *= 3.0;
  474 + test_assert(check_matrix(m_test), "Test source matrix");
  475 + test_assert(check_matrix(test, 3.0, 0.0), "Test GMatrix *= 3.0",
  476 + "Unexpected result matrix:\n"+test.print());
  477 +
  478 + // GMatrix /= 3.0
  479 + test = m_test;
  480 + test /= 3.0;
  481 + test_assert(check_matrix(m_test), "Test source matrix");
  482 + test_assert(check_matrix(test, 1.0/3.0, 0.0), "Test GMatrix /= 3.0",
  483 + "Unexpected result matrix:\n"+test.print());
  484 +
  485 + // GMatrix + GMatrix
  486 + test = m_test + m_test;
  487 + test_assert(check_matrix(m_test), "Test source matrix");
  488 + test_assert(check_matrix(test, 2.0, 0.0), "Test GMatrix + GMatrix",
  489 + "Unexpected result matrix:\n"+test.print());
  490 +
  491 + // GMatrix - GMatrix
  492 + test = m_test - m_test;
  493 + test_assert(check_matrix(m_test), "Test source matrix");
  494 + test_assert(check_matrix(test, 0.0, 0.0), "Test GMatrix - GMatrix",
  495 + "Unexpected result matrix:\n"+test.print());
  496 +
  497 + // GMatrix * 3.0
  498 + test = m_test * 3.0;
  499 + test_assert(check_matrix(m_test), "Test source matrix");
  500 + test_assert(check_matrix(test, 3.0, 0.0), "Test GMatrix * 3.0",
  501 + "Unexpected result matrix:\n"+test.print());
  502 +
  503 + // 3.0 * GMatrix
  504 + test = 3.0 * m_test;
  505 + test_assert(check_matrix(m_test), "Test source matrix");
  506 + test_assert(check_matrix(test, 3.0, 0.0), "Test 3.0 * GMatrix",
  507 + "Unexpected result matrix:\n"+test.print());
  508 +
  509 + // GMatrix / 3.0
  510 + test = m_test / 3.0;
  511 + test_assert(check_matrix(m_test), "Test source matrix");
  512 + test_assert(check_matrix(test, 1.0/3.0, 0.0), "Test GMatrix / 3.0",
  513 + "Unexpected result matrix:\n"+test.print());
  514 +
  515 + // Test invalid matrix addition
  516 + test_try("Test invalid matrix addition");
622 517 try {
623   - GVector v_test7 = bigger*v_test;
624   - test_try_failure();
  518 + test = m_test;
  519 + test += m_bigger;
  520 + test_try_failure("Expected GException::matrix_mismatch exception.");
625 521 }
626   - catch (GException::matrix_vector_mismatch &e) {
  522 + catch (GException::matrix_mismatch &e) {
627 523 test_try_success();
628 524 }
629 525 catch (std::exception &e) {
630 526 test_try_failure(e);
631 527 }
632 528  
  529 + // Return
633 530 return;
634 531 }
635 532  
636   -// Test 8: Matrix*Matrix multiplication
637   -void TestGMatrix::test8(void)
  533 +
  534 +/***********************************************************************//**
  535 + * @brief Test matrix functions
  536 + *
  537 + * Tests matrix functions.
  538 + ***************************************************************************/
  539 +void TestGMatrix::matrix_functions(void)
638 540 {
639   - GMatrix m_test8 = m_test * transpose(m_test);
640   - test_assert(check_matrix_matrix(m_test8) && check_matrix(m_test, 1.0, 0.0),"","Corrupt Matrix*Matrix multiplication");
  541 + // Minimum
  542 + double min = m_test.min();
  543 +
  544 + // Check mimimum
  545 + double value = g_matrix[0];
  546 + for (int row = 0; row < g_rows; ++row) {
  547 + for (int col = 0; col < g_cols; ++col) {
  548 + if (g_matrix[col+row*g_cols] < value) {
  549 + value = g_matrix[col+row*g_cols];
  550 + }
  551 + }
  552 + }
  553 + test_value(min, value, 0.0, "Test minimum function");
  554 +
  555 + // Maximum
  556 + double max = m_test.max();
  557 +
  558 + // Check maximum
  559 + value = g_matrix[0];
  560 + for (int row = 0; row < g_rows; ++row) {
  561 + for (int col = 0; col < g_cols; ++col) {
  562 + if (g_matrix[col+row*g_cols] > value) {
  563 + value = g_matrix[col+row*g_cols];
  564 + }
  565 + }
  566 + }
  567 + test_value(max, value, 0.0, "Test maximum function");
641 568  
642   - test_try("Bigger*Matrix multiplication");
  569 + // Sum
  570 + double sum = m_test.sum();
  571 +
  572 + // Check sum
  573 + value = 0.0;
  574 + for (int row = 0; row < g_rows; ++row) {
  575 + for (int col = 0; col < g_cols; ++col) {
  576 + value += g_matrix[col+row*g_cols];
  577 + }
  578 + }
  579 + test_value(sum, value, 1.0e-20, "Test sum function");
  580 +
  581 + // Transpose function
  582 + GMatrix test1 = transpose(m_test);
  583 + test_assert(check_matrix(m_test), "Test source matrix");
  584 + test_assert(check_matrix_trans(test1, 1.0, 0.0),
  585 + "Test transpose(GMatrix) function",
  586 + "Unexpected transposed matrix:\n"+test1.print());
  587 +
  588 + // Transpose method
  589 + test1 = m_test;
  590 + test1.transpose();
  591 + test_assert(check_matrix(m_test), "Test source matrix");
  592 + test_assert(check_matrix_trans(test1, 1.0, 0.0),
  593 + "Test GMatrix.transpose() method",
  594 + "Unexpected transposed matrix:\n"+test1.print());
  595 +
  596 + // Convert to general matrix
  597 + GMatrix test2 = GMatrix(m_test);
  598 + test_assert(check_matrix(m_test), "Test source matrix");
  599 + test_assert(check_matrix(test2, 1.0, 0.0),
  600 + "Test GMatrix(GMatrix) constructor",
  601 + "Unexpected GMatrix:\n"+test2.print());
  602 +
  603 + // Try extracting of lower triangle
  604 + test_try("Try extracting of lower triangle");
643 605 try {
644   - GMatrix m_test8 = m_test*bigger;
  606 + test2 = m_test.extract_lower_triangle();
  607 + test_try_failure("Expected GException::matrix_not_square exception.");
  608 + }
  609 + catch (GException::matrix_not_square &e) {
645 610 test_try_success();
646 611 }
647 612 catch (std::exception &e) {
648 613 test_try_failure(e);
649 614 }
650 615  
651   - test_try("Bigger*Matrix multiplication 2");
  616 + // Try extracting of upper triangle
  617 + test_try("Try extracting of lower triangle");
652 618 try {
653   - GMatrix m_test8 = bigger*m_test;
654   - test_try_failure();
  619 + test2 = m_test.extract_upper_triangle();
  620 + test_try_failure("Expected GException::matrix_not_square exception.");
655 621 }
656   - catch (GException::matrix_mismatch &e) {
  622 + catch (GException::matrix_not_square &e) {
657 623 test_try_success();
658 624 }
659 625 catch (std::exception &e) {
660 626 test_try_failure(e);
661 627 }
  628 +
  629 + // Extract lower triangle for square matrix
  630 + GMatrix test3 = m_test * transpose(m_test);
  631 + GMatrix test4 = test3.extract_lower_triangle();
  632 + test_assert(check_matrix_lt(test4, test3),
  633 + "Test GMatrix.extract_lower_triangle() method",
  634 + "Unexpected lower triangle:\n"+test4.print());
  635 +
  636 + // Extract upper triangle
  637 + test4 = test3.extract_upper_triangle();
  638 + test_assert(check_matrix_ut(test4, test3),
  639 + "Test GMatrix.extract_upper_triangle() method",
  640 + "Unexpected lower triangle:\n"+test4.print());
662 641  
  642 + // Return
663 643 return;
664 644 }
665 645  
666   -// Test 9: Assignment and arithmetics
667   -void TestGMatrix::test9(void)
  646 +
  647 +/***********************************************************************//**
  648 + * @brief Test matrix comparisons
  649 + ***************************************************************************/
  650 +void TestGMatrix::matrix_compare(void)
668 651 {
669   - //
670   - // -GMatrix
671   - result = -m_test;
672   - test_assert(check_matrix(result, -1.0, 0.0)&&check_matrix(m_test, 1.0, 0.0),"-GMatrix","Corrupt -GMatrix operator.");
673   -
674   - //
675   - // GMatrix += GMatrix
676   - result = m_test;
677   - result += m_test;
678   - test_assert(check_matrix(result, 2.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix += GMatrix","Corrupt GMatrix += GMatrix operator.");
679   -
680   - //
681   - // GMatrix -= GMatrix
682   - result = m_test;
683   - result -= m_test;
684   - test_assert(check_matrix(result, 0.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix -= GMatrix","Corrupt GMatrix -= GMatrix operator.");
685   -
686   - //
687   - // GMatrix *= 3.0
688   - result = m_test;
689   - result *= 3.0;
690   - test_assert(check_matrix(result, 3.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix *= 3.0","Corrupt GMatrix *= double operator.");
691   -
692   - //
693   - // GMatrix /= 3.0
694   - result = m_test;
695   - result /= 3.0;
696   - test_assert(check_matrix(result, 1.0/3.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix /= 3.0","Corrupt GMatrix /= double operator.");
697   -
698   - //
699   - // GMatrix + GMatrix
700   - result = m_test + m_test;
701   - test_assert(check_matrix(result, 2.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix + GMatrix","Corrupt GMatrix + GMatrix operator.");
702   -
703   - //
704   - // GMatrix - GMatrix
705   - result = m_test - m_test;
706   - test_assert(check_matrix(result, 0.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix - GMatrix","Corrupt GMatrix - GMatrix operator.");
707   -
708   - // GMatrix * 3.0
709   - result = m_test * 3.0;
710   - test_assert(check_matrix(result, 3.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix * 3.0","Corrupt GMatrix * double operator.");
711   -
712   - //
713   - // 3.0 * GMatrix
714   - result = 3.0 * m_test;
715   - test_assert(check_matrix(result, 3.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"3.0 * GMatrix","Corrupt double * GMatrix operator.");
716   -
717   - //
718   - // GMatrix / 3.0
719   - result = m_test / 3.0;
720   - test_assert(check_matrix(result, 1.0/3.0, 0.0) && check_matrix(m_test, 1.0, 0.0),"GMatrix / 3.0","Corrupt GMatrix / double operator.");
721   -
722   -
723   - test_try("Test += with bigger");
724   - try {
725   - result = m_test;
726   - result += bigger;
727   - test_try_failure();
728   - }
729   - catch (GException::matrix_mismatch &e) {
730   - test_try_success();
731   - }
732   - catch (std::exception &e) {
733   - test_try_failure(e);
734   - }
  652 + // Allocate an empty matrix
  653 + GMatrix empty(g_rows,g_cols);
  654 +
  655 + // Test operators
  656 + test_assert((m_test == m_test), "Test == operator");
  657 + test_assert(!(m_test == empty), "Test == operator");
  658 + test_assert(!(m_test == m_bigger), "Test == operator");
  659 + test_assert(!(m_test != m_test), "Test != operator");
  660 + test_assert((m_test != empty), "Test != operator");
  661 + test_assert((m_test != m_bigger), "Test != operator");
735 662  
  663 + // Return
736 664 return;
737 665 }
738 666  
739   -// Test 10: Matrix functions
740   -void TestGMatrix::test10(void){
741   - //
742   - double min = m_test.min();
743   - test_assert(check_matrix_min(min),"Test min","Corrupt GMatrix.min() function.");
744 667  
745   - //
746   - double max = m_test.max();
747   - test_assert(check_matrix_max(max),"Test max","Corrupt GMatrix.max() function.");
748   -
749   - //
750   - double sum = m_test.sum();
751   - test_assert(check_matrix_sum(sum),"Test sum","Corrupt GMatrix.sum() function.");
  668 +/***********************************************************************//**
  669 + * @brief Test Cholesky decomposition
  670 + ***************************************************************************/
  671 +/*
  672 +void TestGMatrix::matrix_cholesky(void)
  673 +{
  674 + // Test Cholesky decomposition
  675 + GMatrix cd = cholesky_decompose(m_test);
  676 + GMatrix cd_lower = cd.extract_lower_triangle();
  677 + GMatrix cd_upper = transpose(cd_lower);
  678 + GMatrix cd_product = cd_lower * cd_upper;
  679 + GMatrix cd_residuals = matrix(m_test) - cd_product;
  680 + double res = (abs(cd_residuals)).max();
  681 + test_value(res, 0.0, 1.0e-15, "Test cholesky_decompose() method");
  682 +
  683 + // Test compressed Cholesky decomposition
  684 + GMatrix test_zero = set_matrix_zero();
  685 + GMatrix cd_zero = cholesky_decompose(test_zero);
  686 + GMatrix cd_zero_lower = cd_zero.extract_lower_triangle();
  687 + GMatrix cd_zero_upper = transpose(cd_zero_lower);
  688 + GMatrix cd_zero_product = cd_zero_lower * cd_zero_upper;
  689 + GMatrix cd_zero_residuals = matrix(test_zero) - cd_zero_product;
  690 + res = (abs(cd_zero_residuals)).max();
  691 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_decompose() method");
  692 +
  693 + // Test Cholesky inplace decomposition
  694 + GMatrix test = m_test;
  695 + test.cholesky_decompose();
  696 + GMatrix cd_lower2 = test.extract_lower_triangle();
  697 + test_assert((cd_lower2 == cd_lower), "Test inplace cholesky_decompose() method");
  698 +
  699 + // Test Cholesky solver (first test)
  700 + GVector e0(g_rows);
  701 + GVector a0(g_rows);
  702 + e0[0] = 1.0;
  703 + e0[1] = 0.0;
  704 + e0[2] = 0.0;
  705 + a0[0] = g_matrix[0];
  706 + a0[1] = g_matrix[3];
  707 + a0[2] = g_matrix[6];
  708 + GVector s0 = cd.cholesky_solver(a0) - e0;
  709 + res = max(abs(s0));
  710 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method");
  711 +
  712 + // Test Cholesky solver (second test)
  713 + e0[0] = 0.0;
  714 + e0[1] = 1.0;
  715 + e0[2] = 0.0;
  716 + a0[0] = g_matrix[1];
  717 + a0[1] = g_matrix[4];
  718 + a0[2] = g_matrix[7];
  719 + s0 = cd.cholesky_solver(a0) - e0;
  720 + res = max(abs(s0));
  721 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method");
  722 +
  723 + // Test Cholesky solver (third test)
  724 + e0[0] = 0.0;
  725 + e0[1] = 0.0;
  726 + e0[2] = 1.0;
  727 + a0[0] = g_matrix[2];
  728 + a0[1] = g_matrix[5];
  729 + a0[2] = g_matrix[8];
  730 + s0 = cd.cholesky_solver(a0) - e0;
  731 + res = max(abs(s0));
  732 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method");
  733 +
  734 + // Test compressed Cholesky solver (first test)
  735 + e0 = GVector(g_rows+1);
  736 + a0 = GVector(g_rows+1);
  737 + e0[0] = 1.0;
  738 + e0[1] = 0.0;
  739 + e0[2] = 0.0;
  740 + e0[3] = 0.0;
  741 + a0[0] = g_matrix[0];
  742 + a0[1] = g_matrix[3];
  743 + a0[2] = 0.0;
  744 + a0[3] = g_matrix[6];
  745 + s0 = cd_zero.cholesky_solver(a0) - e0;
  746 + res = max(abs(s0));
  747 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method");
  748 +
  749 + // Test compressed Cholesky solver (second test)
  750 + e0[0] = 0.0;
  751 + e0[1] = 1.0;
  752 + e0[2] = 0.0;
  753 + e0[3] = 0.0;
  754 + a0[0] = g_matrix[1];
  755 + a0[1] = g_matrix[4];
  756 + a0[2] = 0.0;
  757 + a0[3] = g_matrix[7];
  758 + s0 = cd_zero.cholesky_solver(a0) - e0;
  759 + res = max(abs(s0));
  760 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method");
  761 +
  762 + // Test compressed Cholesky solver (third test)
  763 + e0[0] = 0.0;
  764 + e0[1] = 0.0;
  765 + e0[2] = 0.0;
  766 + e0[3] = 1.0;
  767 + a0[0] = g_matrix[2];
  768 + a0[1] = g_matrix[5];
  769 + a0[2] = 0.0;
  770 + a0[3] = g_matrix[8];
  771 + s0 = cd_zero.cholesky_solver(a0) - e0;
  772 + res = max(abs(s0));
  773 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method");
  774 +
  775 + // Test Cholesky inverter
  776 + GMatrix unit(g_rows,g_cols);
  777 + unit(0,0) = unit(1,1) = unit(2,2) = 1.0;
  778 + GMatrix test_inv = m_test;
  779 + test_inv.cholesky_invert();
  780 + GMatrix ci_product = m_test * test_inv;
  781 + GMatrix ci_residuals = ci_product - unit;
  782 + res = (abs(ci_residuals)).max();
  783 + test_value(res, 0.0, 1.0e-15, "Test cholesky_invert method");
  784 +
  785 + // Test Cholesky inverter for compressed matrix
  786 + unit = GMatrix(4,4);
  787 + unit(0,0) = unit(1,1) = unit(3,3) = 1.0;
  788 + GMatrix test_zero_inv = test_zero;
  789 + test_zero_inv.cholesky_invert();
  790 + GMatrix ciz_product = test_zero * test_zero_inv;
  791 + GMatrix ciz_residuals = ciz_product - unit;
  792 + res = (abs(ciz_residuals)).max();
  793 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_invert method");
752 794  
  795 + // Return
753 796 return;
754 797 }
  798 +*/
755 799  
756   -// Test 11: Matrix comparison
757   -void TestGMatrix::test11(void){
758   -
759   - test_assert(((m_test == m_test) == 1),"GMatrix == GMatrix","Corrupt GMatrix == GMatrix operator.");
760   - //
761   - GMatrix m_test10(m_g_rows,m_g_cols);
762   - test_assert(((m_test == m_test10) == 0),"GMatrix == GMatrix 2","Corrupt GMatrix == GMatrix operator.");
763   - //
764   - test_assert(((m_test == bigger) == 0),"GMatrix == GMatrix 3","Corrupt GMatrix == GMatrix operator.");
765   - //
766   - test_assert(((m_test != m_test) == 0),"GMatrix != GMatrix","Corrupt GMatrix != GMatrix operator.");
767   - //
768   - test_assert(((m_test != m_test10) == 1),"GMatrix != GMatrix 2","Corrupt GMatrix != GMatrix operator.");
769   - //
770   - test_assert(((m_test != bigger) == 1),"GMatrix != GMatrix 3","Corrupt GMatrix != GMatrix operator.");
  800 +/***************************************************************************
  801 + * @brief Test matrix printing
  802 + ***************************************************************************/
  803 +void TestGMatrix::matrix_print(void)
  804 +{
  805 + // Set reference
  806 + std::string reference;
  807 + reference.append("=== GMatrix ===\n");
  808 + reference.append(" Number of rows ............: 3\n");
  809 + reference.append(" Number of columns .........: 4\n");
  810 + reference.append(" Number of elements ........: 12\n");
  811 + reference.append(" Number of allocated cells .: 12\n");
  812 + reference.append(" 1, 2, 3, 4\n");
  813 + reference.append(" 5, 6, 7, 8\n");
  814 + reference.append(" 9, 10, 11, 12");
  815 +
  816 + // Test the print method
  817 + std::string output = m_test.print();
  818 + test_assert((output == reference), "Test matrix printing",
  819 + "Unexpected print() output:\n"+output);
771 820  
  821 + // Return
772 822 return;
773 823 }
774 824  
  825 +
  826 +/***************************************************************************
  827 + * @brief Main entry point for test executable
  828 + ***************************************************************************/
775 829 int main(void)
776 830 {
777   - GTestSuites testsuite("GMatrix");
  831 + // Allocate test suit container
  832 + GTestSuites testsuites("GMatrix class testing");
778 833  
779   - bool was_successful=true;
  834 + // Initially assume that we pass all tests
  835 + bool success = true;
780 836  
781   - //Create a test suite
  837 + // Create a test suite
782 838 TestGMatrix test;
783 839  
784   - //Append to the container
785   - testsuite.append(test);
  840 + // Append test suite to the container
  841 + testsuites.append(test);
786 842  
787   - //Run
788   - was_successful=testsuite.run();
  843 + // Run the testsuites
  844 + success = testsuites.run();
789 845  
790   - //save xml report
791   - testsuite.save("reports/GMatrix.xml");
  846 + // Save test report
  847 + testsuites.save("reports/GMatrix.xml");
792 848  
793   - // Return
794   - return was_successful ? 0:1;
  849 + // Return success status
  850 + return (success ? 0 : 1);
795 851 }
... ...
test/test_GMatrix.hpp
  View file @b4ded5b
1 1 /***************************************************************************
2   - * test_GMatrix.hpp - test matrix class *
  2 + * test_GMatrix.hpp - Test generic matrix class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012 by Jean-Baptiste Cayrou *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -18,70 +18,65 @@
18 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19 19 * *
20 20 ***************************************************************************/
  21 +/**
  22 + * @file test_GMatrix.hpp
  23 + * @brief Definition of unit tests for generic matrices
  24 + * @author Juergen Knoedlseder
  25 + */
21 26  
22 27 #ifndef TEST_GMATRIX_HPP
23 28 #define TEST_GMATRIX_HPP
24 29  
25 30 /* __ Includes ___________________________________________________________ */
26 31 #include "GammaLib.hpp"
27   -#include <iostream> // cout, cerr
28 32  
29   -class TestGMatrix : public GTestSuite
30   -{
31   - public:
32   - // Constructors and destructors
33   - TestGMatrix(void) : GTestSuite(){ return; }
34   - virtual ~TestGMatrix(void){ return; }
35 33  
36   - // Methods
37   - virtual void set(void);
38   -
39   - void set_matrix(void);
40   - void set_vector(void);
41   - void set_bigger(void);
42   - void test_output(void);
43   - void test_conversion(void);
44   - void test_extract(void);
45   - void test1(void);
46   - void test2(void);
47   - void test3(void);
48   - void test4(void);
49   - void test5(void);
50   - void test6(void);
51   - void test7(void);
52   - void test8(void);
53   - void test9(void);
54   - void test10(void);
55   - void test11(void);
56   -
57   - // Private methods
58   - private:
59   - void init_matrix(void);
60   - void init_vector(void);
61   -
62   - int check_matrix(const GMatrix& m, const double scale, const double add);
63   - int check_transpose_matrix(const GMatrix& m, const double scale, const double add);
64   - int check_matrix_lt(const GMatrix& m, const double scale, const double add);
65   - int check_matrix_ut(const GMatrix& m, const double scale, const double add);
66   - int check_matrix_vector(const GVector& v);
67   - int check_matrix_matrix(const GMatrix& m);
68   - int check_matrix_min(const double min);
69   - int check_matrix_max(const double max);
70   - int check_matrix_sum(const double sum);
  34 +/***********************************************************************//**
  35 + * @class TestGMatrix
  36 + *
  37 + * @brief Test suite for generic matrix class testing
  38 + *
  39 + * This test suite performs a unit test for the generic matrix class. Tests
  40 + * that match logically are gathered in a number of class methods, and the
  41 + * class methods are executed when running the test.
  42 + ***************************************************************************/
  43 +class TestGMatrix : public GTestSuite {
71 44  
  45 +public:
  46 + // Constructors and destructors
  47 + TestGMatrix(void) : GTestSuite() {}
  48 + virtual ~TestGMatrix(void) {}
72 49  
  50 + // Methods
  51 + virtual void set(void);
  52 + void alloc_matrix(void);
  53 + void assign_values(void);
  54 + void copy_matrix(void);
  55 + void matrix_operations(void);
  56 + void matrix_arithmetics(void);
  57 + void matrix_functions(void);
  58 + void matrix_compare(void);
  59 + //void matrix_cholesky(void);
  60 + void matrix_print(void);
73 61  
  62 +private:
  63 + // Private methods
  64 + GMatrix set_matrix(void) const;
  65 + GMatrix set_matrix_zero(void) const;
  66 + GVector set_vector(void) const;
  67 + bool check_matrix(const GMatrix& matrix,
  68 + const double& scale = 1.0,
  69 + const double& offset = 0.0) const;
  70 + bool check_matrix_trans(const GMatrix& matrix,
  71 + const double& scale,
  72 + const double& offset) const;
  73 + bool check_matrix_lt(const GMatrix& matrix, const GMatrix& ref) const;
  74 + bool check_matrix_ut(const GMatrix& matrix, const GMatrix& ref) const;
74 75  
75   - // Private attributes
76   - GMatrix m_test;
77   - GMatrix bigger;
78   - GMatrix result;
79   - GVector v_test;
80   - double m_g_matrix[12];//
81   - double m_g_vector[4];
82   - int m_g_rows;
83   - int m_g_cols;
84   -
  76 + // Private members;
  77 + GMatrix m_test; //!< Test matrix
  78 + GMatrix m_bigger; //!< Bigger test matrix (for collisions)
  79 + GVector v_test; //!< Test vector
85 80 };
86 81  
87 82 #endif /* TEST_GMATRIX_HPP */
... ...
test/test_GSymMatrix.cpp
  View file @b4ded5b
... ... @@ -268,7 +268,7 @@ void TestGSymMatrix::alloc_matrix(void)
268 268 test_try("Allocate zero matrix");
269 269 try {
270 270 GSymMatrix test(0,0);
271   - test_try_failure();
  271 + test_try_failure("Expected GException::empty exception.");
272 272 }
273 273 catch (GException::empty &e) {
274 274 test_try_success();
... ... @@ -311,7 +311,7 @@ void TestGSymMatrix::assign_values(void)
311 311 test_try("Verify range checking");
312 312 try {
313 313 test(3,3) = 1.0;
314   - test_try_failure();
  314 + test_try_failure("Expected GException::out_of_range exception.");
315 315 }
316 316 catch (GException::out_of_range &e) {
317 317 test_try_success();
... ... @@ -336,7 +336,8 @@ void TestGSymMatrix::copy_matrix(void)
336 336  
337 337 // Test if original and compied matrices are correct
338 338 test_assert(check_matrix(m_test), "Test source matrix");
339   - test_assert(check_matrix(test), "Test matrix copy operator", test.print());
  339 + test_assert(check_matrix(test), "Test matrix copy operator",
  340 + "Found:\n"+test.print()+"\nExpected:\n"+m_test.print());
340 341  
341 342 // Return
342 343 return;
... ... @@ -354,12 +355,14 @@ void TestGSymMatrix::matrix_operations(void)
354 355 GVector test1 = m_test * v_test;
355 356  
356 357 // Check if the result vector is as expected
  358 + GVector ref1 = test1;
357 359 bool result = true;
358 360 for (int row = 0; row < g_rows; ++row) {
359 361 double value = 0.0;
360 362 for (int col = 0; col < g_cols; ++col) {
361 363 value += g_matrix[col+row*g_cols] * g_vector[col];
362 364 }
  365 + ref1[row] = value;
363 366 if (test1[row] != value) {
364 367 result = false;
365 368 break;
... ... @@ -368,7 +371,8 @@ void TestGSymMatrix::matrix_operations(void)
368 371  
369 372 // Test if original matrix and result vector are correct
370 373 test_assert(check_matrix(m_test), "Test source matrix");
371   - test_assert(result, "Test matrix*vector multiplication", test1.print());
  374 + test_assert(result, "Test matrix*vector multiplication",
  375 + "Found:\n"+test1.print()+"\nExpected:\n"+ref1.print());
372 376  
373 377 // Test incompatible vector multiplication
374 378 test_try("Test incompatible matrix*vector multiplication");
... ... @@ -387,6 +391,7 @@ void TestGSymMatrix::matrix_operations(void)
387 391 GSymMatrix test3 = m_test * m_test;
388 392  
389 393 // Check if the result matrix is as expected
  394 + GSymMatrix ref3 = test3;
390 395 result = true;
391 396 for (int row = 0; row < test3.rows(); ++row) {
392 397 for (int col = 0; col < test3.cols(); ++col) {
... ... @@ -394,6 +399,7 @@ void TestGSymMatrix::matrix_operations(void)
394 399 for (int i = 0; i < g_cols; ++i) {
395 400 value += g_matrix[i+row*g_cols] * g_matrix[i+col*g_cols];
396 401 }
  402 + ref3(row,col) = value;
397 403 if (test3(row,col) != value) {
398 404 result = false;
399 405 break;
... ... @@ -403,7 +409,8 @@ void TestGSymMatrix::matrix_operations(void)
403 409  
404 410 // Test if original matrix and result matrix are correct
405 411 test_assert(check_matrix(m_test), "Test source matrix");
406   - test_assert(result, "Test matrix multiplication", test3.print());
  412 + test_assert(result, "Test matrix multiplication",
  413 + "Found:\n"+test3.print()+"\nExpected:\n"+ref3.print());
407 414 test_assert(test3.rows() == g_rows, "Test number of rows of result matrix");
408 415 test_assert(test3.cols() == g_cols, "Test number of columns of result matrix");
409 416  
... ... @@ -411,7 +418,7 @@ void TestGSymMatrix::matrix_operations(void)
411 418 test_try("Test incompatible matrix multiplication");
412 419 try {
413 420 GSymMatrix test4 = m_test * m_bigger;
414   - test_try_failure();
  421 + test_try_failure("Expected GException::matrix_mismatch exception.");
415 422 }
416 423 catch (GException::matrix_mismatch &e) {
417 424 test_try_success();
... ... @@ -424,7 +431,7 @@ void TestGSymMatrix::matrix_operations(void)
424 431 test_try("Test incompatible matrix multiplication");
425 432 try {
426 433 GSymMatrix test5 = m_bigger * m_test;
427   - test_try_failure();
  434 + test_try_failure("Expected GException::matrix_mismatch exception.");
428 435 }
429 436 catch (GException::matrix_mismatch &e) {
430 437 test_try_success();
... ... @@ -514,7 +521,7 @@ void TestGSymMatrix::matrix_arithmetics(void)
514 521 try {
515 522 test = m_test;
516 523 test += m_bigger;
517   - test_try_failure();
  524 + test_try_failure("Expected GException::matrix_mismatch exception.");
518 525 }
519 526 catch (GException::matrix_mismatch &e) {
520 527 test_try_success();
... ... @@ -578,35 +585,38 @@ void TestGSymMatrix::matrix_functions(void)
578 585 // Transpose function
579 586 GSymMatrix test1 = transpose(m_test);
580 587 test_assert(check_matrix(m_test), "Test source matrix");
581   - test_assert(check_matrix(test1, 1.0, 0.0), "Test transpose(GSymMatrix) function",
582   - test1.print());
  588 + test_assert(check_matrix(test1, 1.0, 0.0),
  589 + "Test transpose(GSymMatrix) function",
  590 + "Unexpected transposed matrix:\n"+test1.print());
583 591  
584 592 // Transpose method
585 593 test1 = m_test;
586 594 test1.transpose();
587 595 test_assert(check_matrix(m_test), "Test source matrix");
588   - test_assert(check_matrix(test1, 1.0, 0.0), "Test GSymMatrix.transpose() method",
589   - test1.print());
  596 + test_assert(check_matrix(test1, 1.0, 0.0),
  597 + "Test GSymMatrix.transpose() method",
  598 + "Unexpected transposed matrix:\n"+test1.print());
590 599  
591 600 // Convert to general matrix
592 601 GMatrix test2 = GMatrix(m_test);
593 602 test_assert(check_matrix(m_test), "Test source matrix");
594   - test_assert(check_matrix(test2, 1.0, 0.0), "Test GMatrix(GSymMatrix) constructor",
595   - test2.print());
  603 + test_assert(check_matrix(test2, 1.0, 0.0),
  604 + "Test GMatrix(GSymMatrix) constructor",
  605 + "Unexpected GMatrix:\n"+test2.print());
596 606  
597 607 // Extract lower triangle
598 608 test2 = m_test.extract_lower_triangle();
599 609 test_assert(check_matrix(m_test), "Test source matrix");
600 610 test_assert(check_matrix_lt(test2, 1.0, 0.0),
601 611 "Test GSymMatrix.extract_lower_triangle() method",
602   - test2.print());
  612 + "Unexpected GMatrix:\n"+test2.print());
603 613  
604 614 // Extract upper triangle
605 615 test2 = m_test.extract_upper_triangle();
606 616 test_assert(check_matrix(m_test), "Test source matrix");
607 617 test_assert(check_matrix_ut(test2, 1.0, 0.0),
608 618 "Test GSymMatrix.extract_upper_triangle() method",
609   - test2.print());
  619 + "Unexpected GMatrix:\n"+test2.print());
610 620  
611 621 // Return
612 622 return;
... ... @@ -644,7 +654,7 @@ void TestGSymMatrix::matrix_cholesky(void)
644 654 GMatrix cd_lower = cd.extract_lower_triangle();
645 655 GMatrix cd_upper = transpose(cd_lower);
646 656 GMatrix cd_product = cd_lower * cd_upper;
647   - GMatrix cd_residuals = matrix(m_test) - cd_product;
  657 + GMatrix cd_residuals = GMatrix(m_test) - cd_product;
648 658 double res = (abs(cd_residuals)).max();
649 659 test_value(res, 0.0, 1.0e-15, "Test cholesky_decompose() method");
650 660  
... ... @@ -654,7 +664,7 @@ void TestGSymMatrix::matrix_cholesky(void)
654 664 GMatrix cd_zero_lower = cd_zero.extract_lower_triangle();
655 665 GMatrix cd_zero_upper = transpose(cd_zero_lower);
656 666 GMatrix cd_zero_product = cd_zero_lower * cd_zero_upper;
657   - GMatrix cd_zero_residuals = matrix(test_zero) - cd_zero_product;
  667 + GMatrix cd_zero_residuals = GMatrix(test_zero) - cd_zero_product;
658 668 res = (abs(cd_zero_residuals)).max();
659 669 test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_decompose() method");
660 670  
... ... @@ -745,8 +755,8 @@ void TestGSymMatrix::matrix_cholesky(void)
745 755 unit(0,0) = unit(1,1) = unit(2,2) = 1.0;
746 756 GSymMatrix test_inv = m_test;
747 757 test_inv.cholesky_invert();
748   - GSymMatrix ci_product = m_test * test_inv;
749   - GSymMatrix ci_residuals = ci_product - unit;
  758 + GMatrix ci_product = m_test * test_inv;
  759 + GMatrix ci_residuals = ci_product - unit;
750 760 res = (abs(ci_residuals)).max();
751 761 test_value(res, 0.0, 1.0e-15, "Test cholesky_invert method");
752 762  
... ... @@ -755,8 +765,8 @@ void TestGSymMatrix::matrix_cholesky(void)
755 765 unit(0,0) = unit(1,1) = unit(3,3) = 1.0;
756 766 GSymMatrix test_zero_inv = test_zero;
757 767 test_zero_inv.cholesky_invert();
758   - GSymMatrix ciz_product = test_zero * test_zero_inv;
759   - GSymMatrix ciz_residuals = ciz_product - unit;
  768 + GMatrix ciz_product = test_zero * test_zero_inv;
  769 + GMatrix ciz_residuals = ciz_product - unit;
760 770 res = (abs(ciz_residuals)).max();
761 771 test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_invert method");
762 772  
... ... @@ -783,7 +793,8 @@ void TestGSymMatrix::matrix_print(void)
783 793  
784 794 // Test the print method
785 795 std::string output = m_test.print();
786   - test_assert((output == reference), "Test matrix printing", output);
  796 + test_assert((output == reference), "Test matrix printing",
  797 + "Unexpected print() output:\n"+output);
787 798  
788 799 // Return
789 800 return;
... ...