Jürgen Knödlseder / gammalib

Commit 973bbf0b466cf93b8671dafedffb59fed43416f7

Authored by Jürgen Knödlseder
1 parent ada01abc

Use unit testsuite for testing. 

Adapt the unit test for the usage of the new GammaLib test suite. Test reports are now reported in a JUnit compliant XML file.
Inline Side-by-side
Showing 2 changed files with 1154 additions and 1563 deletions
test/test_GSparseMatrix.cpp
  View file @973bbf0
1 /*************************************************************************** 1 /***************************************************************************
2 - * test_GSparseMatrix.cpp - test GSparseMatrix class * 2 + * test_GSparseMatrix.cpp - Test sparse 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 * This program is free software: you can redistribute it and/or modify * 7 * This program is free software: you can redistribute it and/or modify *
@@ -20,23 +20,16 @@ @@ -20,23 +20,16 @@
20 ***************************************************************************/ 20 ***************************************************************************/
21 /** 21 /**
22 * @file test_GSparseMatrix.cpp 22 * @file test_GSparseMatrix.cpp
23 - * @brief Testing of GSparseMatrix class  
24 - * @author J. Knodlseder 23 + * @brief Implementation of unit tests for sparse matrices
  24 + * @author Juergen Knoedlseder
25 */ 25 */
26 26
27 /* __ Includes ___________________________________________________________ */ 27 /* __ Includes ___________________________________________________________ */
28 -#include <cstdlib> 28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
29 #include <cmath> 31 #include <cmath>
30 -#include <iostream> // cout, cerr  
31 -#include <stdexcept> // std::exception  
32 -#include "GammaLib.hpp"  
33 -  
34 -/* __ Namespaces _________________________________________________________ */  
35 -using namespace std;  
36 -  
37 -/* __ Definitions ________________________________________________________ */  
38 -//#define DUMP_RESIDUALS // Dump residuals  
39 -#define DUMP_TIMING // Dump timing 32 +#include "test_GSparseMatrix.hpp"
40 33
41 /* __ Globals ____________________________________________________________ */ 34 /* __ Globals ____________________________________________________________ */
42 double g_matrix[] = {1.0, 7.0, 3.0, 2.0, 4.0, 8.0, 5.0, 6.0, 9.0}; 35 double g_matrix[] = {1.0, 7.0, 3.0, 2.0, 4.0, 8.0, 5.0, 6.0, 9.0};
@@ -49,1207 +42,1033 @@ int g_cols = 5; @@ -49,1207 +42,1033 @@ int g_cols = 5;
49 42
50 43
51 /*************************************************************************** 44 /***************************************************************************
52 - * Set test matrix * 45 + * @brief Set test matrix
53 ***************************************************************************/ 46 ***************************************************************************/
54 -GSparseMatrix set_matrix(void) 47 +GSparseMatrix TestGSparseMatrix::set_matrix(void) const
55 { 48 {
56 - try {  
57 - GSparseMatrix matrix(g_rows,g_cols);  
58 - for (int i = 0; i < g_elements; ++i)  
59 - matrix(g_row[i],g_col[i]) = g_matrix[i];  
60 - return matrix;  
61 - }  
62 - catch (exception &e) {  
63 - cout << "TEST ERROR: Unable to set test matrix." << endl;  
64 - throw;  
65 - }  
66 - return GSparseMatrix(g_rows,g_cols);  
67 -} 49 + // Allocate matrix
  50 + GSparseMatrix matrix(g_rows, g_cols);
68 51
  52 + // Set matrix values
  53 + for (int i = 0; i < g_elements; ++i) {
  54 + matrix(g_row[i],g_col[i]) = g_matrix[i];
  55 + }
69 56
70 -/***************************************************************************  
71 - * Set test vector *  
72 - ***************************************************************************/  
73 -GVector set_vector(void)  
74 -{  
75 - try {  
76 - GVector vector(g_cols);  
77 - for (int col = 0; col < g_cols; ++col)  
78 - vector[col] = g_vector[col];  
79 - return vector;  
80 - }  
81 - catch (exception &e) {  
82 - cout << "TEST ERROR: Unable to set test vector." << endl;  
83 - throw;  
84 - }  
85 - return GVector(g_cols); 57 + // Return matrix
  58 + return matrix;
86 } 59 }
87 60
88 61
89 /*************************************************************************** 62 /***************************************************************************
90 - * Check test matrix * 63 + * @brief Set test matrix with zero lines
91 ***************************************************************************/ 64 ***************************************************************************/
92 -int check_matrix(const GSparseMatrix& m, const double scale, const double add) 65 +/*
  66 +GSparseMatrix TestGSparseMatrix::set_matrix_zero(void) const
93 { 67 {
94 - int result = 1;  
95 - try {  
96 - for (int row = 0; row < g_rows; ++row) {  
97 - for (int col = 0; col < g_cols; ++col) {  
98 - int i;  
99 - double ref_value = add;  
100 - for (i = 0; i < g_elements; ++i) {  
101 - if (g_row[i] == row && g_col[i] == col) {  
102 - ref_value = g_matrix[i] * scale + add;  
103 - break;  
104 - }  
105 - }  
106 - if (abs(m(row,col)-ref_value) > 1.0e-15) {  
107 - result = 0;  
108 - break;  
109 - }  
110 - }  
111 - }  
112 - }  
113 - catch (exception &e) {  
114 - cout << "TEST ERROR: Unable to check test matrix." << endl;  
115 - throw;  
116 - }  
117 - return result;  
118 -} 68 + // Allocate matrix
  69 + GSparseMatrix 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++;
  82 + }
119 83
  84 + // Return matrix
  85 + return matrix;
  86 +}
  87 +*/
120 88
121 /*************************************************************************** 89 /***************************************************************************
122 - * Check matrix*vector * 90 + * @brief Set test vector
123 ***************************************************************************/ 91 ***************************************************************************/
124 -int check_matrix_vector(const GVector& v) 92 +GVector TestGSparseMatrix::set_vector(void) const
125 { 93 {
126 - int result = 1;  
127 - try {  
128 - for (int row = 0; row < v.size(); ++row) {  
129 - double value = 0.0;  
130 - for (int col = 0; col < g_cols; ++col) {  
131 - int i;  
132 - double ref_value = 0.0;  
133 - for (i = 0; i < g_elements; ++i) {  
134 - if (g_row[i] == row && g_col[i] == col) {  
135 - ref_value = g_matrix[i];  
136 - break;  
137 - }  
138 - }  
139 - value += ref_value * g_vector[col];  
140 - }  
141 - if (v[row] != value) {  
142 - result = 0;  
143 - break;  
144 - } 94 + // Allocate vector
  95 + GVector vector(g_cols);
  96 +
  97 + // Set vector values
  98 + for (int col = 0; col < g_cols; ++col) {
  99 + vector[col] = g_vector[col];
145 } 100 }
146 - }  
147 - catch (exception &e) {  
148 - cout << "TEST ERROR: Unable to check matrix*vector product." << endl;  
149 - throw;  
150 - }  
151 - return result; 101 +
  102 + // Return vector
  103 + return vector;
152 } 104 }
153 105
154 106
155 /*************************************************************************** 107 /***************************************************************************
156 - * Check matrix*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.
157 ***************************************************************************/ 116 ***************************************************************************/
158 -int check_matrix_matrix(const GSparseMatrix& m) 117 +bool TestGSparseMatrix::check_matrix(const GSparseMatrix& matrix,
  118 + const double& scale,
  119 + const double& offset) const
159 { 120 {
160 - if (m.rows() != g_rows || m.cols() != g_rows)  
161 - return 0;  
162 - int result = 1;  
163 - try {  
164 - for (int row = 0; row < m.rows(); ++row) {  
165 - for (int col = 0; col < m.cols(); ++col) {  
166 - double value = 0.0;  
167 - for (int i = 0; i < g_cols; ++i) {  
168 - int k;  
169 - double ref_value_1 = 0.0;  
170 - double ref_value_2 = 0.0;  
171 - for (k = 0; k < g_elements; ++k) {  
172 - if (g_row[k] == row && g_col[k] == i) {  
173 - ref_value_1 = g_matrix[k];  
174 - break;  
175 - }  
176 - }  
177 - for (k = 0; k < g_elements; ++k) {  
178 - if (g_row[k] == col && g_col[k] == i) {  
179 - ref_value_2 = g_matrix[k];  
180 - break;  
181 - }  
182 - }  
183 - value += ref_value_1 * ref_value_2;  
184 - }  
185 -//cout << row << "," << col << ": " << abs(m(row,col)-value) << endl;  
186 - if (m(row,col) != value) {  
187 - result = 0;  
188 - break;  
189 - }  
190 - }  
191 - }  
192 - }  
193 - catch (exception &e) {  
194 - cout << "TEST ERROR: Unable to check matrix*matrix product." << endl;  
195 - throw;  
196 - }  
197 - return result; 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 +
  128 + // Set reference value
  129 + double ref_value = offset;
  130 + for (int i = 0; i < g_elements; ++i) {
  131 + if (g_row[i] == row && g_col[i] == col) {
  132 + ref_value = g_matrix[i] * scale + offset;
  133 + break;
  134 + }
  135 + }
  136 +
  137 + // Check matrix
  138 + if (std::abs(matrix(row,col)-ref_value) > 1.0e-15) {
  139 + result = false;
  140 + break;
  141 + }
  142 +
  143 + } // endfor: looped over columns
  144 + } // endfor: looped over rows
  145 +
  146 + // Return result
  147 + return result;
198 } 148 }
199 149
200 150
201 /*************************************************************************** 151 /***************************************************************************
202 - * Check transposed test matrix * 152 + * @brief Check if transposed matrix corresponds to test matrix
  153 + *
  154 + * @param[in] matrix Trasposed matrix.
  155 + * @param[in] scale scale factor (default: 1.0).
  156 + * @param[in] offset offset value (default: 0.0).
  157 + *
  158 + * Checks if a matrix corresponds to the test matrix. Optionally, the test
  159 + * matrix can be scaled and an offset can be added.
203 ***************************************************************************/ 160 ***************************************************************************/
204 -int check_transpose_matrix(const GSparseMatrix& m, const double scale, const double add) 161 +bool TestGSparseMatrix::check_matrix_trans(const GSparseMatrix& matrix,
  162 + const double& scale,
  163 + const double& offset) const
205 { 164 {
206 - int result = 1;  
207 - try { 165 + // Initialise check with true
  166 + bool result = true;
  167 +
  168 + // Compare all matrix elements
208 for (int row = 0; row < g_rows; ++row) { 169 for (int row = 0; row < g_rows; ++row) {
209 - for (int col = 0; col < g_cols; ++col) {  
210 - int i;  
211 - double ref_value = 0.0;  
212 - for (i = 0; i < g_elements; ++i) {  
213 - if (g_row[i] == row && g_col[i] == col) {  
214 - ref_value = g_matrix[i] * scale + add;  
215 - break;  
216 - }  
217 - }  
218 - if (abs(m(col,row)-ref_value) > 1.0e-15) {  
219 - result = 0;  
220 - break;  
221 - }  
222 - }  
223 - }  
224 - }  
225 - catch (exception &e) {  
226 - cout << "TEST ERROR: Unable to check test matrix." << endl;  
227 - throw;  
228 - }  
229 - return result;  
230 -} 170 + for (int col = 0; col < g_cols; ++col) {
  171 +
  172 + // Set reference value
  173 + double ref_value = offset;
  174 + for (int i = 0; i < g_elements; ++i) {
  175 + if (g_row[i] == row && g_col[i] == col) {
  176 + ref_value = g_matrix[i] * scale + offset;
  177 + break;
  178 + }
  179 + }
231 180
  181 + // Check matrix
  182 + if (std::abs(matrix(col,row)-ref_value) > 1.0e-15) {
  183 + result = false;
  184 + break;
  185 + }
  186 +
  187 + } // endfor: looped over columns
  188 + } // endfor: looped over rows
232 189
233 -/***************************************************************************  
234 - * Check matrix min *  
235 - ***************************************************************************/  
236 -int check_matrix_min(const double min)  
237 -{  
238 - double value = g_matrix[0];  
239 - for (int i = 0; i < g_elements; ++i) {  
240 - if (g_matrix[i] < value)  
241 - value = g_matrix[i];  
242 - }  
243 -// return (min == value);  
244 - return (min == 0.0); 190 + // Return result
  191 + return result;
245 } 192 }
246 193
247 194
248 /*************************************************************************** 195 /***************************************************************************
249 - * Check matrix max * 196 + * @brief Check lower triangle general matrix
  197 + *
  198 + * @param[in] matrix Matrix.
  199 + * @param[in] ref Reference matrix.
250 ***************************************************************************/ 200 ***************************************************************************/
251 -int check_matrix_max(const double max) 201 +bool TestGSparseMatrix::check_matrix_lt(const GSparseMatrix& matrix,
  202 + const GSparseMatrix& ref) const
252 { 203 {
253 - double value = g_matrix[0];  
254 - for (int i = 0; i < g_elements; ++i) {  
255 - if (g_matrix[i] > value)  
256 - value = g_matrix[i];  
257 - }  
258 - return (max == value); 204 + // Initialise check with true
  205 + bool result = true;
  206 +
  207 + // Compare all matrix elements
  208 + for (int row = 0; row < g_rows; ++row) {
  209 + for (int col = 0; col < g_cols; ++col) {
  210 +
  211 + // Set reference value
  212 + double ref_value = 0.0;
  213 + if (col <= row) {
  214 + for (int i = 0; i < g_elements; ++i) {
  215 + if (g_row[i] == row && g_col[i] == col) {
  216 + ref_value = g_matrix[i];
  217 + break;
  218 + }
  219 + }
  220 + }
  221 +
  222 + // Check matrix
  223 + if (std::abs(matrix(row,col)-ref_value) > 1.0e-15) {
  224 + result = false;
  225 + break;
  226 + }
  227 +
  228 + } // endfor: looped over columns
  229 + } // endfor: looped over rows
  230 +
  231 + // Return result
  232 + return result;
259 } 233 }
260 234
261 235
262 /*************************************************************************** 236 /***************************************************************************
263 - * Check matrix sum * 237 + * @brief Check upper triangle general matrix
  238 + *
  239 + * @param[in] matrix Matrix.
  240 + * @param[in] ref Reference matrix.
264 ***************************************************************************/ 241 ***************************************************************************/
265 -int check_matrix_sum(const double sum) 242 +bool TestGSparseMatrix::check_matrix_ut(const GSparseMatrix& matrix,
  243 + const GSparseMatrix& ref) const
266 { 244 {
267 - double value = 0.0;  
268 - for (int i = 0; i < g_elements; ++i)  
269 - value += g_matrix[i];  
270 - return (sum == value); 245 + // Initialise check with true
  246 + bool result = true;
  247 +
  248 + // Compare all matrix elements
  249 + for (int row = 0; row < g_rows; ++row) {
  250 + for (int col = 0; col < g_cols; ++col) {
  251 +
  252 + // Set reference value
  253 + double ref_value = 0.0;
  254 + if (col >= row) {
  255 + for (int i = 0; i < g_elements; ++i) {
  256 + if (g_row[i] == row && g_col[i] == col) {
  257 + ref_value = g_matrix[i];
  258 + break;
  259 + }
  260 + }
  261 + }
  262 +
  263 + // Check matrix
  264 + if (std::abs(matrix(row,col)-ref_value) > 1.0e-15) {
  265 + result = false;
  266 + break;
  267 + }
  268 +
  269 + } // endfor: looped over columns
  270 + } // endfor: looped over rows
  271 +
  272 + // Return result
  273 + return result;
271 } 274 }
272 275
273 276
274 -/***************************************************************************  
275 - * Test: Output *  
276 - ***************************************************************************/  
277 -void test_output(const GSparseMatrix& m_test) 277 +/***********************************************************************//**
  278 + * @brief Set parameters and tests
  279 + **************************************************************************/
  280 +void TestGSparseMatrix::set(void)
278 { 281 {
279 - cout << "Test GSparseMatrix: Output test matrix: " << endl;  
280 - try {  
281 - cout << m_test << endl;  
282 - }  
283 - catch (exception &e) {  
284 - cout << endl << "TEST ERROR: Unable to output matrix." << endl;  
285 - cout << e.what() << endl;  
286 - throw;  
287 - } 282 + // Set test name
  283 + name("GSparseMatrix");
  284 +
  285 + // Add tests
  286 + append(static_cast<pfunction>(&TestGSparseMatrix::alloc_matrix), "Test matrix allocation");
  287 + append(static_cast<pfunction>(&TestGSparseMatrix::assign_values), "Test value assignment");
  288 + append(static_cast<pfunction>(&TestGSparseMatrix::copy_matrix), "Test matrix copying");
  289 + append(static_cast<pfunction>(&TestGSparseMatrix::matrix_operations), "Test matrix operations");
  290 + append(static_cast<pfunction>(&TestGSparseMatrix::matrix_arithmetics), "Test matrix arithmetics");
  291 + append(static_cast<pfunction>(&TestGSparseMatrix::matrix_functions), "Test matrix functions");
  292 + append(static_cast<pfunction>(&TestGSparseMatrix::matrix_compare), "Test matrix comparisons");
  293 + append(static_cast<pfunction>(&TestGSparseMatrix::matrix_cholesky), "Test matrix Cholesky decomposition");
  294 + append(static_cast<pfunction>(&TestGSparseMatrix::matrix_print), "Test matrix printing");
  295 +
  296 + // Set members
  297 + m_test = set_matrix();
  298 + v_test = set_vector();
  299 + m_bigger = GSparseMatrix(g_rows+1, g_cols+1);
  300 +
  301 + // Return
  302 + return;
288 } 303 }
289 304
290 305
291 -/***************************************************************************  
292 - * Test: Allocate zero matrix * 306 +/***********************************************************************//**
  307 + * @brief Test matrix allocation
293 ***************************************************************************/ 308 ***************************************************************************/
294 -void test_allocate(void) 309 +void TestGSparseMatrix::alloc_matrix(void)
295 { 310 {
296 - cout << "Test GSparseMatrix: Allocate zero matrix: ";  
297 - try {  
298 - GSparseMatrix test(0,0);  
299 - }  
300 - catch (GException::empty &e) {  
301 - cout << ". ok." << endl;  
302 - }  
303 - catch (exception &e) {  
304 - cout << endl << "TEST ERROR: Did not signal zero matrix allocation." << endl;  
305 - cout << e.what() << endl;  
306 - throw;  
307 - }  
308 - cout << ". ok." << endl; 311 + // Allocate zero matrix. The allocation should fail.
  312 + test_try("Allocate zero matrix");
  313 + try {
  314 + GSparseMatrix test(0,0);
  315 + test_try_failure("Expected GException::empty exception.");
  316 + }
  317 + catch (GException::empty &e) {
  318 + test_try_success();
  319 + }
  320 + catch (std::exception &e) {
  321 + test_try_failure(e);
  322 + }
  323 +
  324 + // Setup a symmetric sparse matrix
  325 + int size = 30;
  326 + GSparseMatrix symmetric(size,size);
  327 + for (int i = 0; i < size; i+=2) {
  328 + for (int j = 0; j < size; j+=2) {
  329 + symmetric(i,j) = 1.0+i+j;
  330 + }
  331 + }
  332 +
  333 + // Convert to GMatrix
  334 + test_try("Test symmetric GMatrix conversion");
  335 + try {
  336 + GMatrix cnv_matrix = GMatrix(symmetric);
  337 + GSparseMatrix back_matrix = GSparseMatrix(cnv_matrix);
  338 + test_assert((symmetric == back_matrix),
  339 + "Test symmetric GSparseMatrix <-> GMatrix conversion",
  340 + "Found:\n"+back_matrix.print()+"\nExpected:\n"+symmetric.print());
  341 + test_try_success();
  342 + }
  343 + catch (GException::empty &e) {
  344 + test_try_success();
  345 + }
  346 + catch (std::exception &e) {
  347 + test_try_failure(e);
  348 + }
  349 +
  350 + // Test GMatrix conversion
  351 + test_try("Test GMatrix conversion");
  352 + try {
  353 + GMatrix cnv_matrix = GMatrix(m_test);
  354 + GSparseMatrix back_matrix = GSparseMatrix(cnv_matrix);
  355 + test_assert((m_test == back_matrix),
  356 + "Test GSparseMatrix <-> GMatrix conversion",
  357 + "Found:\n"+back_matrix.print()+"\nExpected:\n"+m_test.print());
  358 + test_try_success();
  359 + }
  360 + catch (std::exception &e) {
  361 + test_try_failure(e);
  362 + }
  363 +
  364 + // Test GSparseMatrix <-> GSymMatrix conversion
  365 + test_try("Test GSymMatrix conversion");
  366 + try {
  367 + GSymMatrix cnv_sym = GSymMatrix(symmetric);
  368 + GSparseMatrix back_sym = GSparseMatrix(cnv_sym);
  369 + test_assert((symmetric == back_sym),
  370 + "Test GSparseMatrix <-> GSymMatrix conversion",
  371 + "Found:\n"+back_sym.print()+"\nExpected:\n"+symmetric.print());
  372 + test_try_success();
  373 + }
  374 + catch (std::exception &e) {
  375 + test_try_failure(e);
  376 + }
  377 +
  378 + // Test invalid GSparseMatrix <-> GSymMatrix conversion
  379 + test_try("Test invalid GSymMatrix conversion");
  380 + try {
  381 + GSymMatrix bad_sym = GSymMatrix(m_test);
  382 + test_try_failure("Expected GException::matrix_not_symmetric exception.");
  383 + }
  384 + catch (GException::matrix_not_symmetric &e) {
  385 + test_try_success();
  386 + }
  387 + catch (std::exception &e) {
  388 + test_try_failure(e);
  389 + }
  390 +
  391 + // Return
  392 + return;
309 } 393 }
310 394
311 395
312 -/***************************************************************************  
313 - * Test: Assign values * 396 +/***********************************************************************//**
  397 + * @brief Test value assignment
314 ***************************************************************************/ 398 ***************************************************************************/
315 -void test_assign_values(const GSparseMatrix& m_test) 399 +void TestGSparseMatrix::assign_values(void)
316 { 400 {
317 - //  
318 - cout << "Test GSparseMatrix: Assign matrix values: ";  
319 - try {  
320 - // 401 + // Setup 3x3 matrix
321 GSparseMatrix test(3,3); 402 GSparseMatrix test(3,3);
322 - test(1,1) = 1.0;  
323 - if (test(0,0) != 0.0 || test(1,0) != 0.0 || test(2,0) != 0.0 ||  
324 - test(0,1) != 0.0 || test(1,1) != 1.0 || test(2,1) != 0.0 ||  
325 - test(0,2) != 0.0 || test(1,2) != 0.0 || test(2,2) != 0.0) {  
326 - cout << endl << "TEST ERROR: Unable to set matrix values." << endl;  
327 - cout << test << endl;  
328 - throw;  
329 - }  
330 - cout << ".";  
331 - //  
332 - GSparseMatrix result = m_test;  
333 - result(0,0) = 0.0;  
334 - result(0,0) = 1.0;  
335 - if (!check_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
336 - cout << endl << "TEST ERROR: Bad M(row,col) assignment." << endl;  
337 - cout << "result:" << endl << result << endl;  
338 - cout << "m_test:" << endl << m_test << endl;  
339 - throw;  
340 - }  
341 - cout << ".";  
342 - //  
343 - result = m_test;  
344 - result(0,0) = 0.0;  
345 - result(0,4) = 0.0;  
346 - result(0,0) = 1.0;  
347 - result(0,4) = 0.0;  
348 - result(0,4) = 7.0;  
349 - if (!check_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
350 - cout << endl << "TEST ERROR: Bad M(row,col) assignment." << endl;  
351 - cout << "result:" << endl << result << endl;  
352 - cout << "m_test:" << endl << m_test << endl;  
353 - throw;  
354 - }  
355 - cout << ".";  
356 - //  
357 - result = m_test;  
358 - result(0,1) = 2.0;  
359 - result(0,2) = 2.0;  
360 - result(0,3) = 2.0;  
361 - result(0,1) = 0.0;  
362 - result(0,2) = 0.0;  
363 - result(0,3) = 0.0;  
364 - if (!check_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
365 - cout << endl << "TEST ERROR: Bad M(row,col) assignment." << endl;  
366 - cout << "result:" << endl << result << endl;  
367 - cout << "m_test:" << endl << m_test << endl;  
368 - throw;  
369 - }  
370 - cout << ".";  
371 - //  
372 - result.clear();  
373 - if (!check_matrix(result, 0.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
374 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.clear()." << endl;  
375 - cout << "result:" << endl << result << endl;  
376 - cout << "m_test:" << endl << m_test << endl;  
377 - throw;  
378 - }  
379 - cout << ".";  
380 - }  
381 - catch (exception &e) {  
382 - cout << e.what() << endl;  
383 - throw;  
384 - }  
385 - #if defined(G_RANGE_CHECK) 403 +
  404 + // Assignment individual values
  405 + for (int i = 0; i < 3; ++i) {
  406 + for (int k = 0; k < 3; ++k) {
  407 + double value = i*2.0 + k*2.0;
  408 + test(i,k) = value;
  409 + }
  410 + }
  411 +
  412 + // Check assignment of individual values
  413 + for (int i = 0; i < 3; ++i) {
  414 + for (int k = 0; k < 3; ++k) {
  415 + double value = i*2.0 + k*2.0;
  416 + test_value(test(i,k), value, 1.0e-10, "Test matrix element assignment");
  417 + }
  418 + }
  419 +
  420 + // Verify range checking
  421 + #ifdef G_RANGE_CHECK
  422 + test_try("Verify range checking");
386 try { 423 try {
387 - GSparseMatrix test(3,3);  
388 test(3,3) = 1.0; 424 test(3,3) = 1.0;
  425 + test_try_failure("Expected GException::out_of_range exception.");
389 } 426 }
390 catch (GException::out_of_range &e) { 427 catch (GException::out_of_range &e) {
  428 + test_try_success();
391 } 429 }
392 - catch (exception &e) {  
393 - cout << e.what() << endl;  
394 - throw; 430 + catch (std::exception &e) {
  431 + test_try_failure(e);
395 } 432 }
396 - cout << ".";  
397 #endif 433 #endif
398 - //  
399 - // Insert column into 10 x 10 matrix using matrix stack and the  
400 - // compressed vector format  
401 - try {  
402 - //  
403 - // Set-up sparse matrix for adding 434 +
  435 + // Setup 10x10 matrix and keep for reference
404 GSparseMatrix sparse(10,10); 436 GSparseMatrix sparse(10,10);
405 - for (int i = 3; i < 5; ++i)  
406 - sparse(i,i) = 5.0; 437 + for (int i = 3; i < 5; ++i) {
  438 + sparse(i,i) = 5.0;
  439 + }
  440 + GSparseMatrix initial = sparse;
407 GSparseMatrix reference = sparse; 441 GSparseMatrix reference = sparse;
408 - //  
409 - // Initialise vector for column adding  
410 - GVector column(10);  
411 - //  
412 - // Initialise matrix stack 442 +
  443 + // Insert column into 10 x 10 matrix using large matrix stack and the
  444 + // add_col(GVector) method
413 sparse.stack_init(100,50); 445 sparse.stack_init(100,50);
414 - //  
415 - // Fill matrix using add_col 446 + GVector column(10);
416 for (int j = 0; j < 10; ++j) { 447 for (int j = 0; j < 10; ++j) {
417 - //  
418 - // Set column index  
419 - int col = int(0.8 * j + 0.5); // This allows that some columns are twice  
420 - if (col > 9) col -= 10; // This avoids overflow  
421 - //  
422 - // Set-up vector  
423 - int i_min = (j < 2) ? 0 : j-2;  
424 - int i_max = (j > 8) ? 10 : j+2;  
425 - column = 0.0;  
426 - for (int i = i_min; i < i_max; ++i) {  
427 - column[i] = (i+1)*1;  
428 - reference(i,col) += column[i];  
429 - }  
430 - //  
431 - // Add vector  
432 - sparse.add_col(column, col); 448 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  449 + if (col > 9) col -= 10; // This avoids overflow
  450 + int i_min = (j < 2) ? 0 : j-2;
  451 + int i_max = (j > 8) ? 10 : j+2;
  452 + column = 0.0;
  453 + for (int i = i_min; i < i_max; ++i) {
  454 + column[i] = (i+1)*1;
  455 + reference(i,col) += column[i];
  456 + }
  457 + sparse.add_col(column, col);
433 } 458 }
434 - //  
435 - // Destroy stack  
436 sparse.stack_destroy(); 459 sparse.stack_destroy();
437 - if (sparse != reference) {  
438 - cout << endl << "TEST ERROR: Corrupt add_col function (vector version)." << endl;  
439 - cout << "Sparse " << sparse << endl;  
440 - cout << "Reference " << reference << endl;  
441 - throw; 460 + test_assert((sparse == reference),
  461 + "Test stack fill with large stack using add_col(GVector) method",
  462 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  463 +
  464 + // Insert column into 10 x 10 matrix using small matrix stack and the
  465 + // add_col(GVector) method
  466 + sparse = initial;
  467 + sparse.stack_init(100,3);
  468 + for (int j = 0; j < 10; ++j) {
  469 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  470 + if (col > 9) col -= 10; // This avoids overflow
  471 + int i_min = (j < 2) ? 0 : j-2;
  472 + int i_max = (j > 8) ? 10 : j+2;
  473 + column = 0.0;
  474 + for (int i = i_min; i < i_max; ++i) {
  475 + column[i] = (i+1)*1;
  476 + }
  477 + sparse.add_col(column, col);
442 } 478 }
443 - }  
444 - catch (exception &e) {  
445 - cout << e.what() << endl;  
446 - throw;  
447 - }  
448 - //  
449 - // Insert column into 10 x 10 matrix using matrix stack and the  
450 - // compressed vector format  
451 - try {  
452 - //  
453 - // Set-up sparse matrix for adding  
454 - GSparseMatrix sparse(10,10);  
455 - sparse.clear();  
456 - for (int i = 3; i < 5; ++i)  
457 - sparse(i,i) = 5.0;  
458 - GSparseMatrix reference = sparse;  
459 - //  
460 - // Set-up workspace for column adding  
461 - double* wrk_data = new double[10];  
462 - int* wrk_row = new int[10];  
463 - //  
464 - // Initialise matrix stack  
465 - sparse.stack_init(100,50);  
466 - //  
467 - // Fill matrix using add_col 479 + sparse.stack_destroy();
  480 + test_assert((sparse == reference),
  481 + "Test stack fill with small stack using add_col(GVector) method",
  482 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  483 +
  484 + // Insert column into 10 x 10 matrix using tiny matrix stack and the
  485 + // add_col(GVector) method
  486 + sparse = initial;
  487 + sparse.stack_init(8,3);
468 for (int j = 0; j < 10; ++j) { 488 for (int j = 0; j < 10; ++j) {
469 - //  
470 - // Set column index  
471 - int col = int(0.8 * j + 0.5); // This allows that some columns are twice  
472 - if (col > 9) col -= 10; // This avoids overflow  
473 - //  
474 - // Set-up vector  
475 - int inx = 0;  
476 - int i_min = (j < 2) ? 0 : j-2;  
477 - int i_max = (j > 8) ? 10 : j+2;  
478 - for (int i = i_min; i < i_max; ++i) {  
479 - wrk_data[inx] = (i+1)*1;  
480 - wrk_row[inx] = i;  
481 - reference(i,col) += wrk_data[inx];  
482 - inx++;  
483 - }  
484 - //  
485 - // Add vector  
486 - sparse.add_col(wrk_data, wrk_row, inx, col); 489 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  490 + if (col > 9) col -= 10; // This avoids overflow
  491 + int i_min = (j < 2) ? 0 : j-2;
  492 + int i_max = (j > 8) ? 10 : j+2;
  493 + column = 0.0;
  494 + for (int i = i_min; i < i_max; ++i) {
  495 + column[i] = (i+1)*1;
  496 + }
  497 + sparse.add_col(column, col);
487 } 498 }
488 - //  
489 - // Destroy stack  
490 sparse.stack_destroy(); 499 sparse.stack_destroy();
491 - if (sparse != reference) {  
492 - cout << endl << "TEST ERROR: Corrupt add_col function (compressed array version)." << endl;  
493 - cout << "Sparse " << sparse << endl;  
494 - cout << "Reference " << reference << endl;  
495 - throw; 500 + test_assert((sparse == reference),
  501 + "Test stack fill with tiny stack using add_col(GVector) method",
  502 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  503 +
  504 + // Insert column into 10 x 10 matrix using no matrix stack and the
  505 + // add_col(GVector) method
  506 + sparse = initial;
  507 + for (int j = 0; j < 10; ++j) {
  508 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  509 + if (col > 9) col -= 10; // This avoids overflow
  510 + int i_min = (j < 2) ? 0 : j-2;
  511 + int i_max = (j > 8) ? 10 : j+2;
  512 + column = 0.0;
  513 + for (int i = i_min; i < i_max; ++i) {
  514 + column[i] = (i+1)*1;
  515 + }
  516 + sparse.add_col(column, col);
496 } 517 }
497 - // 518 + sparse.stack_destroy();
  519 + test_assert((sparse == reference),
  520 + "Test fill using add_col(GVector) method",
  521 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  522 +
  523 + // Set-up workspace for compressed column adding
  524 + double* wrk_data = new double[10];
  525 + int* wrk_row = new int[10];
  526 +
  527 + // Compressed tests
  528 + test_try("Verify compressed column add_col() method");
  529 + try {
  530 +
  531 + // Insert column into 10 x 10 matrix using large matrix stack and the
  532 + // compressed column add_col() method
  533 + sparse = initial;
  534 + reference = initial;
  535 + sparse.stack_init(100,50);
  536 + for (int j = 0; j < 10; ++j) {
  537 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  538 + if (col > 9) col -= 10; // This avoids overflow
  539 + int inx = 0;
  540 + int i_min = (j < 3) ? 0 : j-3;
  541 + int i_max = (j > 8) ? 10 : j+2;
  542 + for (int i = i_min; i < i_max; ++i) {
  543 + wrk_data[inx] = (i+1)*3.7;
  544 + wrk_row[inx] = i;
  545 + reference(i,col) += wrk_data[inx];
  546 + inx++;
  547 + }
  548 + sparse.add_col(wrk_data, wrk_row, inx, col);
  549 + }
  550 + sparse.stack_destroy();
  551 + test_assert((sparse == reference),
  552 + "Test stack fill with large stack using compressed add_col() method",
  553 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  554 +
  555 + // Insert column into 10 x 10 matrix using small matrix stack and the
  556 + // compressed column add_col() method
  557 + sparse = initial;
  558 + sparse.stack_init(100,2);
  559 + for (int j = 0; j < 10; ++j) {
  560 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  561 + if (col > 9) col -= 10; // This avoids overflow
  562 + int inx = 0;
  563 + int i_min = (j < 3) ? 0 : j-3;
  564 + int i_max = (j > 8) ? 10 : j+2;
  565 + for (int i = i_min; i < i_max; ++i) {
  566 + wrk_data[inx] = (i+1)*3.7;
  567 + wrk_row[inx] = i;
  568 + inx++;
  569 + }
  570 + sparse.add_col(wrk_data, wrk_row, inx, col);
  571 + }
  572 + sparse.stack_destroy();
  573 + test_assert((sparse == reference),
  574 + "Test stack fill with small stack using compressed add_col() method",
  575 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  576 +
  577 + // Insert column into 10 x 10 matrix using tiny matrix stack and the
  578 + // compressed column add_col() method
  579 + sparse = initial;
  580 + sparse.stack_init(3,2);
  581 + for (int j = 0; j < 10; ++j) {
  582 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  583 + if (col > 9) col -= 10; // This avoids overflow
  584 + int inx = 0;
  585 + int i_min = (j < 3) ? 0 : j-3;
  586 + int i_max = (j > 8) ? 10 : j+2;
  587 + for (int i = i_min; i < i_max; ++i) {
  588 + wrk_data[inx] = (i+1)*3.7;
  589 + wrk_row[inx] = i;
  590 + inx++;
  591 + }
  592 + sparse.add_col(wrk_data, wrk_row, inx, col);
  593 + }
  594 + sparse.stack_destroy();
  595 + test_assert((sparse == reference),
  596 + "Test stack fill with tiny stack using compressed add_col() method",
  597 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  598 +
  599 + // Insert column into 10 x 10 matrix using no matrix stack and the
  600 + // compressed column add_col() method
  601 + sparse = initial;
  602 + for (int j = 0; j < 10; ++j) {
  603 + int col = int(0.8 * j + 0.5); // This allows that some columns are twice
  604 + if (col > 9) col -= 10; // This avoids overflow
  605 + int inx = 0;
  606 + int i_min = (j < 3) ? 0 : j-3;
  607 + int i_max = (j > 8) ? 10 : j+2;
  608 + for (int i = i_min; i < i_max; ++i) {
  609 + wrk_data[inx] = (i+1)*3.7;
  610 + wrk_row[inx] = i;
  611 + inx++;
  612 + }
  613 + sparse.add_col(wrk_data, wrk_row, inx, col);
  614 + }
  615 + sparse.stack_destroy();
  616 + test_assert((sparse == reference),
  617 + "Test fill using compressed add_col() method",
  618 + "Found:\n"+sparse.print()+"\nExpected:\n"+reference.print());
  619 +
  620 + // Signal success
  621 + test_try_success();
  622 + }
  623 + catch (std::exception &e) {
  624 + test_try_failure(e);
  625 + }
  626 +
498 // Free workspace 627 // Free workspace
499 delete [] wrk_data; 628 delete [] wrk_data;
500 delete [] wrk_row; 629 delete [] wrk_row;
501 - }  
502 - catch (exception &e) {  
503 - cout << e.what() << endl;  
504 - throw;  
505 - }  
506 - cout << " ok." << endl; 630 +
  631 + // Return
  632 + return;
507 } 633 }
508 634
509 635
510 -/***************************************************************************  
511 - * Test: Matrix copy constructor * 636 +/***********************************************************************//**
  637 + * @brief Test matrix copy
512 ***************************************************************************/ 638 ***************************************************************************/
513 -void test_copy_constructor(const GSparseMatrix& m_test) 639 +void TestGSparseMatrix::copy_matrix(void)
514 { 640 {
515 - cout << "Test GSparseMatrix: Define matrix using copy constructor: ";  
516 - try {  
517 - // 641 + // Copy matrix
518 GSparseMatrix test = m_test; 642 GSparseMatrix test = m_test;
519 - if (!check_matrix(test, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
520 - cout << endl << "TEST ERROR: Corrupt copy constructor." << endl;  
521 - cout << test << endl;  
522 - throw;  
523 - }  
524 - }  
525 - catch (exception &e) {  
526 - cout << e.what() << endl;  
527 - throw;  
528 - }  
529 - cout << ". ok." << endl; 643 +
  644 + // Test if original and compied matrices are correct
  645 + test_assert(check_matrix(m_test), "Test source matrix");
  646 + test_assert(check_matrix(test), "Test matrix copy operator",
  647 + "Found:\n"+test.print()+"\nExpected:\n"+m_test.print());
  648 +
  649 + // Return
  650 + return;
530 } 651 }
531 652
532 653
533 -/***************************************************************************  
534 - * Test: Matrix assignment * 654 +/***********************************************************************//**
  655 + * @brief Test matrix operations
  656 + *
  657 + * Tests matrix*vector and matrix*matrix multiplication operations.
535 ***************************************************************************/ 658 ***************************************************************************/
536 -void test_assign_matrix(const GSparseMatrix& m_test) 659 +void TestGSparseMatrix::matrix_operations(void)
537 { 660 {
538 - cout << "Test GSparseMatrix: Matrix assignment: ";  
539 - try {  
540 - //  
541 - // GSparseMatrix = GSparseMatrix  
542 - GSparseMatrix result(1,1);  
543 - result = m_test;  
544 - if (!check_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
545 - cout << endl << "TEST ERROR: Corrupt assignment." << endl;  
546 - cout << result << endl;  
547 - throw;  
548 - }  
549 - cout << ".";  
550 - //  
551 - // GSparseMatrix = GSparseMatrix (bigger matrix)  
552 - GSparseMatrix bigger(1000,1000);  
553 - for (int i = 0; i < 1000; ++i)  
554 - bigger(i,i) = (i+1)*1.1;  
555 - result = bigger;  
556 - cout << ".";  
557 - }  
558 - catch (exception &e) {  
559 - cout << e.what() << endl;  
560 - throw;  
561 - }  
562 - cout << " ok." << endl;  
563 -} 661 + // Perform vector multiplication
  662 + GVector test1 = m_test * v_test;
564 663
  664 + // Check result
  665 + GVector ref1(g_rows);
  666 + for (int i = 0; i < g_elements; ++i) {
  667 + ref1[g_row[i]] += g_matrix[i] * v_test[g_col[i]];
  668 + }
  669 + bool result = true;
  670 + for (int i = 0; i < g_rows; ++i) {
  671 + if (ref1[i] != test1[i]) {
  672 + result = false;
  673 + break;
  674 + }
  675 + }
565 676
566 -/***************************************************************************  
567 - * Test: Matrix transpose *  
568 - ***************************************************************************/  
569 -void test_transpose(const GSparseMatrix& m_test)  
570 -{  
571 - cout << "Test GSparseMatrix: Matrix transpose: ";  
572 - try {  
573 - //  
574 - // transpose(GSparseMatrix)  
575 - GSparseMatrix result(1,1);  
576 - result = transpose(m_test);  
577 - if (!check_transpose_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
578 - cout << endl << "TEST ERROR: Corrupt transpose(GSparseMatrix) function." << endl;  
579 - cout << result << endl;  
580 - throw;  
581 - }  
582 - cout << ".";  
583 - //  
584 - // GSparseMatrix.transpose()  
585 - result = m_test;  
586 - result.transpose();  
587 - if (!check_transpose_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
588 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.transpose() function." << endl;  
589 - cout << result << endl;  
590 - throw;  
591 - }  
592 - cout << ".";  
593 - }  
594 - catch (exception &e) {  
595 - cout << e.what() << endl;  
596 - throw;  
597 - }  
598 - cout << " ok." << endl;  
599 -} 677 + // Test if original matrix and result vector are correct
  678 + test_assert(check_matrix(m_test), "Test source matrix");
  679 + test_assert(result, "Test matrix*vector multiplication",
  680 + "Found:\n"+test1.print()+"\nExpected:\n"+ref1.print());
600 681
  682 + // Test incompatible vector multiplication
  683 + test_try("Test incompatible matrix*vector multiplication");
  684 + try {
  685 + GVector test2 = m_bigger * v_test;
  686 + test_try_failure("Expected GException::matrix_vector_mismatch exception.");
  687 + }
  688 + catch (GException::matrix_vector_mismatch &e) {
  689 + test_try_success();
  690 + }
  691 + catch (std::exception &e) {
  692 + test_try_failure(e);
  693 + }
  694 +
  695 + // Test matrix multiplication
  696 + GSparseMatrix test3 = m_test * transpose(m_test);
601 697
602 -/***************************************************************************  
603 - * Test: Matrix*Vector multiplication *  
604 - ***************************************************************************/  
605 -void test_matrix_vector(const GSparseMatrix& m_test, const GVector& v_test)  
606 -{  
607 - cout << "Test GSparseMatrix: Matrix*Vector multiplication: ";  
608 - try {  
609 - //  
610 - GVector test = m_test*v_test;  
611 - if (!check_matrix_vector(test) || !check_matrix(m_test, 1.0, 0.0)) {  
612 - cout << endl << "TEST ERROR: Corrupt Matrix*Vector multiplication." << endl;  
613 - cout << test << endl;  
614 - throw;  
615 - }  
616 - }  
617 - catch (exception &e) {  
618 - cout << e.what() << endl;  
619 - throw;  
620 - }  
621 - cout << ".";  
622 - try {  
623 - GSparseMatrix bigger(1000,1000);  
624 - for (int i = 0; i < 1000; ++i)  
625 - bigger(i,i) = (i+1)*1.1;  
626 - GVector test = bigger*v_test;  
627 - }  
628 - catch (GException::matrix_vector_mismatch &e) {  
629 - }  
630 - catch (exception &e) {  
631 - cout << e.what() << endl;  
632 - throw;  
633 - }  
634 - cout << ". ok." << endl;  
635 -} 698 + // Check if the result matrix is as expected
  699 + GSparseMatrix ref3(g_rows, g_rows);
  700 + for (int row = 0; row < g_rows; ++row) {
  701 + for (int col = 0; col < g_rows; ++col) {
  702 + double value = 0.0;
  703 + for (int i = 0; i < g_cols; ++i) {
  704 + double ref_value_1 = 0.0;
  705 + double ref_value_2 = 0.0;
  706 + for (int k = 0; k < g_elements; ++k) {
  707 + if (g_row[k] == row && g_col[k] == i) {
  708 + ref_value_1 = g_matrix[k];
  709 + break;
  710 + }
  711 + }
  712 + for (int k = 0; k < g_elements; ++k) {
  713 + if (g_row[k] == col && g_col[k] == i) {
  714 + ref_value_2 = g_matrix[k];
  715 + break;
  716 + }
  717 + }
  718 + value += ref_value_1 * ref_value_2;
  719 + }
  720 + ref3(row,col) = value;
  721 + if (test3(row,col) != value) {
  722 + result = false;
  723 + break;
  724 + }
  725 + }
  726 + }
636 727
  728 + // Test if original matrix and result matrix are correct
  729 + test_assert(check_matrix(m_test), "Test source matrix");
  730 + test_assert(result, "Test matrix multiplication",
  731 + "Found:\n"+test3.print()+"\nExpected:\n"+ref3.print());
  732 + test_value(test3.rows(), g_rows, "Test number of rows of result matrix");
  733 + test_value(test3.cols(), g_rows, "Test number of columns of result matrix");
637 734
638 -/***************************************************************************  
639 - * Test: Matrix*Vector multiplication *  
640 - ***************************************************************************/  
641 -void test_matrix_matrix(const GSparseMatrix& m_test)  
642 -{  
643 - cout << "Test GSparseMatrix: Matrix*Matrix multiplication: ";  
644 - try {  
645 - //  
646 - GSparseMatrix test = m_test * transpose(m_test);  
647 - if (!check_matrix_matrix(test) || !check_matrix(m_test, 1.0, 0.0)) {  
648 - cout << endl << "TEST ERROR: Corrupt Matrix*Matrix multiplication." << endl;  
649 - cout << "test:" << endl << test << endl;  
650 - cout << "m_test:" << endl << m_test << endl;  
651 - throw;  
652 - }  
653 - cout << ".";  
654 - //  
655 - GSparseMatrix bigger(1000,1000);  
656 - for (int i = 0; i < 1000; ++i)  
657 - bigger(i,i) = (i+1)*1.1; 735 + // Test incompatible matrix multiplication
  736 + test_try("Test incompatible matrix multiplication");
658 try { 737 try {
659 - GSparseMatrix test = m_test * bigger; 738 + GSparseMatrix test4 = m_bigger * m_test;
  739 + test_try_failure("Expected GException::matrix_mismatch exception.");
660 } 740 }
661 catch (GException::matrix_mismatch &e) { 741 catch (GException::matrix_mismatch &e) {
662 - cout << "."; 742 + test_try_success();
663 } 743 }
664 - catch (exception &e) {  
665 - cout << e.what() << endl;  
666 - throw; 744 + catch (std::exception &e) {
  745 + test_try_failure(e);
667 } 746 }
668 - // 747 +
  748 + // Test another incompatible matrix multiplication
  749 + test_try("Test incompatible matrix multiplication");
669 try { 750 try {
670 - GSparseMatrix test7 = bigger * m_test; 751 + GSparseMatrix test5 = m_bigger * m_test;
  752 + test_try_failure("Expected GException::matrix_mismatch exception.");
671 } 753 }
672 catch (GException::matrix_mismatch &e) { 754 catch (GException::matrix_mismatch &e) {
673 - cout << "."; 755 + test_try_success();
674 } 756 }
675 - catch (exception &e) {  
676 - cout << e.what() << endl;  
677 - throw; 757 + catch (std::exception &e) {
  758 + test_try_failure(e);
678 } 759 }
679 - }  
680 - catch (exception &e) {  
681 - cout << e.what() << endl;  
682 - throw;  
683 - }  
684 - cout << " ok." << endl; 760 +
  761 + // Return
  762 + return;
685 } 763 }
686 764
687 765
688 -/***************************************************************************  
689 - * Test: Arithmetics * 766 +/***********************************************************************//**
  767 + * @brief Test matrix arithmetics
  768 + *
  769 + * Tests matrix arithmetics.
690 ***************************************************************************/ 770 ***************************************************************************/
691 -void test_arithmetics(const GSparseMatrix& m_test) 771 +void TestGSparseMatrix::matrix_arithmetics(void)
692 { 772 {
693 - cout << "Test GSparseMatrix: Arithmetics: ";  
694 - try {  
695 - //  
696 // -GSparseMatrix 773 // -GSparseMatrix
697 - GSparseMatrix result(1,1);  
698 - result = -m_test;  
699 - if (!check_matrix(result, -1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
700 - cout << endl << "TEST ERROR: Corrupt -GSparseMatrix operator." << endl;  
701 - cout << "result " << result << endl;  
702 - cout << "m_test " << m_test << endl;  
703 - throw;  
704 - }  
705 - cout << ".";  
706 - //  
707 - // abs(GSparseMatrix)  
708 - result = abs(m_test);  
709 - if (!check_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
710 - cout << endl << "TEST ERROR: Corrupt abs(GSparseMatrix) function." << endl;  
711 - cout << "result " << result << endl;  
712 - cout << "m_test " << m_test << endl;  
713 - throw;  
714 - }  
715 - cout << ".";  
716 - // 774 + GSparseMatrix test = -m_test;
  775 + test_assert(check_matrix(m_test), "Test source matrix");
  776 + test_assert(check_matrix(test, -1.0, 0.0), "Test -GSparseMatrix",
  777 + "Unexpected result matrix:\n"+test.print());
  778 +
717 // GSparseMatrix += GSparseMatrix 779 // GSparseMatrix += GSparseMatrix
718 - result = m_test;  
719 - result += m_test;  
720 - if (!check_matrix(result, 2.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
721 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix += GSparseMatrix operator." << endl;  
722 - cout << "result " << result << endl;  
723 - cout << "m_test " << m_test << endl;  
724 - throw;  
725 - }  
726 - cout << ".";  
727 - // 780 + test = m_test;
  781 + test += m_test;
  782 + test_assert(check_matrix(m_test), "Test source matrix");
  783 + test_assert(check_matrix(test, 2.0, 0.0), "Test GSparseMatrix += GSparseMatrix",
  784 + "Unexpected result matrix:\n"+test.print());
  785 +
728 // GSparseMatrix -= GSparseMatrix 786 // GSparseMatrix -= GSparseMatrix
729 - result = m_test;  
730 - result -= m_test;  
731 - if (!check_matrix(result, 0.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
732 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix -= GSparseMatrix operator." << endl;  
733 - cout << "result " << result << endl;  
734 - cout << "m_test " << m_test << endl;  
735 - throw;  
736 - }  
737 - cout << ".";  
738 - // 787 + test = m_test;
  788 + test -= m_test;
  789 + test_assert(check_matrix(m_test), "Test source matrix");
  790 + test_assert(check_matrix(test, 0.0, 0.0), "Test GSparseMatrix -= GSparseMatrix",
  791 + "Unexpected result matrix:\n"+test.print());
  792 +
739 // GSparseMatrix *= 3.0 793 // GSparseMatrix *= 3.0
740 - result = m_test;  
741 - result *= 3.0;  
742 - if (!check_matrix(result, 3.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
743 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix *= double operator." << endl;  
744 - cout << "result " << result << endl;  
745 - cout << "m_test " << m_test << endl;  
746 - throw;  
747 - }  
748 - cout << ".";  
749 - // 794 + test = m_test;
  795 + test *= 3.0;
  796 + test_assert(check_matrix(m_test), "Test source matrix");
  797 + test_assert(check_matrix(test, 3.0, 0.0), "Test GSparseMatrix *= 3.0",
  798 + "Unexpected result matrix:\n"+test.print());
  799 +
750 // GSparseMatrix /= 3.0 800 // GSparseMatrix /= 3.0
751 - result = m_test;  
752 - result /= 3.0;  
753 - if (!check_matrix(result, 1.0/3.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
754 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix /= double operator." << endl;  
755 - cout << "result " << result << endl;  
756 - cout << "m_test " << m_test << endl;  
757 - throw;  
758 - }  
759 - cout << ".";  
760 - // 801 + test = m_test;
  802 + test /= 3.0;
  803 + test_assert(check_matrix(m_test), "Test source matrix");
  804 + test_assert(check_matrix(test, 1.0/3.0, 0.0), "Test GSparseMatrix /= 3.0",
  805 + "Unexpected result matrix:\n"+test.print());
  806 +
761 // GSparseMatrix + GSparseMatrix 807 // GSparseMatrix + GSparseMatrix
762 - result = m_test + m_test;  
763 - if (!check_matrix(result, 2.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
764 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix + GSparseMatrix operator." << endl;  
765 - cout << "result " << result << endl;  
766 - cout << "m_test " << m_test << endl;  
767 - throw;  
768 - }  
769 - cout << ".";  
770 - //  
771 - // GSparseMatrix - GSparseMatrix  
772 - result = m_test - m_test;  
773 - if (!check_matrix(result, 0.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
774 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix - GSparseMatrix operator." << endl;  
775 - cout << "result " << result << endl;  
776 - cout << "m_test " << m_test << endl;  
777 - throw;  
778 - }  
779 - cout << ".";  
780 - // 808 + test = m_test + m_test;
  809 + test_assert(check_matrix(m_test), "Test source matrix");
  810 + test_assert(check_matrix(test, 2.0, 0.0), "Test GSparseMatrix + GSparseMatrix",
  811 + "Unexpected result matrix:\n"+test.print());
  812 +
781 // GSparseMatrix - GSparseMatrix 813 // GSparseMatrix - GSparseMatrix
782 - result = m_test - m_test;  
783 - if (!check_matrix(result, 0.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
784 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix - GSparseMatrix operator." << endl;  
785 - cout << "result " << result << endl;  
786 - cout << "m_test " << m_test << endl;  
787 - throw;  
788 - }  
789 - cout << ".";  
790 - // 814 + test = m_test - m_test;
  815 + test_assert(check_matrix(m_test), "Test source matrix");
  816 + test_assert(check_matrix(test, 0.0, 0.0), "Test GSparseMatrix - GSparseMatrix",
  817 + "Unexpected result matrix:\n"+test.print());
  818 +
791 // GSparseMatrix * 3.0 819 // GSparseMatrix * 3.0
792 - result = m_test * 3.0;  
793 - if (!check_matrix(result, 3.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
794 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix * double operator." << endl;  
795 - cout << "result " << result << endl;  
796 - cout << "m_test " << m_test << endl;  
797 - throw;  
798 - }  
799 - cout << ".";  
800 - //  
801 - // GSparseMatrix * -7.12  
802 - result = m_test * -7.12;  
803 - if (!check_matrix(result, -7.12, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
804 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix * double operator." << endl;  
805 - cout << "result " << result << endl;  
806 - cout << "m_test " << m_test << endl;  
807 - throw;  
808 - }  
809 - cout << ".";  
810 - // 820 + test = m_test * 3.0;
  821 + test_assert(check_matrix(m_test), "Test source matrix");
  822 + test_assert(check_matrix(test, 3.0, 0.0), "Test GSparseMatrix * 3.0",
  823 + "Unexpected result matrix:\n"+test.print());
  824 +
811 // 3.0 * GSparseMatrix 825 // 3.0 * GSparseMatrix
812 - result = 3.0 * m_test;  
813 - if (!check_matrix(result, 3.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
814 - cout << endl << "TEST ERROR: Corrupt double * GSparseMatrix operator." << endl;  
815 - cout << "result " << result << endl;  
816 - cout << "m_test " << m_test << endl;  
817 - throw;  
818 - }  
819 - cout << ".";  
820 - // 826 + test = 3.0 * m_test;
  827 + test_assert(check_matrix(m_test), "Test source matrix");
  828 + test_assert(check_matrix(test, 3.0, 0.0), "Test 3.0 * GSparseMatrix",
  829 + "Unexpected result matrix:\n"+test.print());
  830 +
821 // GSparseMatrix / 3.0 831 // GSparseMatrix / 3.0
822 - result = m_test / 3.0;  
823 - if (!check_matrix(result, 1.0/3.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
824 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix / double operator." << endl;  
825 - cout << "result " << result << endl;  
826 - cout << "m_test " << m_test << endl;  
827 - throw;  
828 - }  
829 - cout << ".";  
830 - //  
831 - // GSparseMatrix / 0.0  
832 - double zero = 0.0;  
833 - result = m_test / zero;  
834 - if (!check_matrix(result, 1.0/zero, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
835 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix / 0.0." << endl;  
836 - cout << "result " << result << endl;  
837 - cout << "m_test " << m_test << endl;  
838 - throw;  
839 - cout << ".";  
840 - }  
841 - //  
842 - // GSparseMatrix.add_col()  
843 - result = m_test;  
844 - GVector v_add(g_rows);  
845 - for (int i = 0; i < g_rows; ++i)  
846 - v_add[i] = 7.0;  
847 - for (int col = 0; col < g_cols; ++col)  
848 - result.add_col(v_add, col);  
849 - if (!check_matrix(result, 1.0, 7.0) || !check_matrix(m_test, 1.0, 0.0)) {  
850 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.add_col(+v)." << endl;  
851 - cout << "result " << result << endl;  
852 - cout << "m_test " << m_test << endl;  
853 - throw;  
854 - }  
855 - cout << ".";  
856 - //  
857 - // GSparseMatrix.add_col()  
858 - for (int col = 0; col < g_cols; ++col)  
859 - result.add_col(-v_add, col);  
860 - if (!check_matrix(result, 1.0, 0.0) || !check_matrix(m_test, 1.0, 0.0)) {  
861 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.add_col(-v)." << endl;  
862 - cout << "result " << result << endl;  
863 - cout << "m_test " << m_test << endl;  
864 - throw;  
865 - }  
866 - cout << ".";  
867 - //  
868 - // GSparseMatrix.add_col() using a compressed array  
869 - GSparseMatrix compare = m_test;  
870 - GVector v_sparse(g_rows);  
871 - for (int i = 0; i < g_rows; i += 2) // Build vector for comparison  
872 - v_sparse[i] = 7.0;  
873 - for (int col = 0; col < g_cols; ++col)  
874 - compare.add_col(v_sparse, col);  
875 - result = m_test;  
876 - double* values = new double[g_rows];  
877 - int* rows = new int[g_rows];  
878 - int number = 0;  
879 - for (int i = 0; i < g_rows; i += 2, ++number) {  
880 - values[number] = 7.0;  
881 - rows[number] = i;  
882 - }  
883 - for (int col = 0; col < g_cols; ++col)  
884 - result.add_col(values, rows, number, col);  
885 - delete [] values;  
886 - delete [] rows;  
887 - if (result != compare) {  
888 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.add_col(compressed array)." << endl;  
889 - cout << "m_test " << m_test << endl;  
890 - for (int i = 0; i < g_rows; ++i) {  
891 - for (int j = 0; j < g_cols; ++j) {  
892 - std::cout << " " << m_test(i,j);  
893 - }  
894 - std::cout << std::endl;  
895 - }  
896 - cout << "result " << result << endl;  
897 - for (int i = 0; i < g_rows; ++i) {  
898 - for (int j = 0; j < g_cols; ++j) {  
899 - std::cout << " " << result(i,j);  
900 - }  
901 - std::cout << std::endl;  
902 - }  
903 - cout << "compare " << compare << endl;  
904 - for (int i = 0; i < g_rows; ++i) {  
905 - for (int j = 0; j < g_cols; ++j) {  
906 - std::cout << " " << compare(i,j);  
907 - }  
908 - std::cout << std::endl;  
909 - }  
910 - throw;  
911 - }  
912 - cout << ".";  
913 - }  
914 - catch (exception &e) {  
915 - cout << e.what() << endl;  
916 - throw;  
917 - }  
918 - try {  
919 - GSparseMatrix bigger(1000,1000);  
920 - GSparseMatrix result(1,1);  
921 - for (int i = 0; i < 1000; ++i)  
922 - bigger(i,i) = (i+1)*1.1;  
923 - result = m_test;  
924 - result += bigger;  
925 - }  
926 - catch (GException::matrix_mismatch &e) {  
927 - cout << ".";  
928 - }  
929 - catch (exception &e) {  
930 - cout << e.what() << endl;  
931 - throw;  
932 - }  
933 - cout << " ok." << endl;  
934 -} 832 + test = m_test / 3.0;
  833 + test_assert(check_matrix(m_test), "Test source matrix");
  834 + test_assert(check_matrix(test, 1.0/3.0, 0.0), "Test GSparseMatrix / 3.0",
  835 + "Unexpected result matrix:\n"+test.print());
935 836
  837 + // Test invalid matrix addition
  838 + test_try("Test invalid matrix addition");
  839 + try {
  840 + test = m_test;
  841 + test += m_bigger;
  842 + test_try_failure("Expected GException::matrix_mismatch exception.");
  843 + }
  844 + catch (GException::matrix_mismatch &e) {
  845 + test_try_success();
  846 + }
  847 + catch (std::exception &e) {
  848 + test_try_failure(e);
  849 + }
936 850
937 -/***************************************************************************  
938 - * Test: Matrix functions *  
939 - ***************************************************************************/  
940 -void test_functions(const GSparseMatrix& m_test)  
941 -{  
942 - cout << "Test GSparseMatrix: Matrix functions: ";  
943 - try {  
944 - //  
945 - double min = m_test.min();  
946 - if (!check_matrix_min(min)) {  
947 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.min() function." << endl;  
948 - cout << "min:" << min << endl;  
949 - throw;  
950 - }  
951 - cout << ".";  
952 - //  
953 - double max = m_test.max();  
954 - if (!check_matrix_max(max)) {  
955 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.max() function." << endl;  
956 - cout << "max:" << max << endl;  
957 - throw;  
958 - }  
959 - cout << ".";  
960 - //  
961 - double sum = m_test.sum();  
962 - if (!check_matrix_sum(sum)) {  
963 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.sum() function." << endl;  
964 - cout << "sum:" << sum << endl;  
965 - throw;  
966 - }  
967 - cout << ".";  
968 - }  
969 - catch (exception &e) {  
970 - cout << e.what() << endl;  
971 - }  
972 - cout << " ok." << endl; 851 + // Return
  852 + return;
973 } 853 }
974 854
975 855
976 -/***************************************************************************  
977 - * Test: Matrix comparison * 856 +/***********************************************************************//**
  857 + * @brief Test matrix functions
  858 + *
  859 + * Tests matrix functions.
978 ***************************************************************************/ 860 ***************************************************************************/
979 -void test_comparison(const GSparseMatrix& m_test) 861 +void TestGSparseMatrix::matrix_functions(void)
980 { 862 {
981 - cout << "Test GSparseMatrix: Comparison: ";  
982 - try {  
983 - GSparseMatrix bigger(1000,1000);  
984 - for (int i = 0; i < 1000; ++i)  
985 - bigger(i,i) = (i+1)*1.1;  
986 - //  
987 - if ((m_test == m_test) != 1) {  
988 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix == GSparseMatrix operator." << endl;  
989 - throw;  
990 - }  
991 - cout << ".";  
992 - //  
993 - GSparseMatrix test(g_rows,g_cols);  
994 - if ((m_test == test) != 0) {  
995 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix == GSparseMatrix operator." << endl;  
996 - throw;  
997 - }  
998 - cout << ".";  
999 - //  
1000 - if ((m_test == bigger) != 0) {  
1001 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix == GSparseMatrix operator." << endl;  
1002 - throw;  
1003 - }  
1004 - cout << ".";  
1005 - //  
1006 - if ((m_test != m_test) != 0) {  
1007 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix != GSparseMatrix operator." << endl;  
1008 - throw;  
1009 - }  
1010 - cout << ".";  
1011 - //  
1012 - if ((m_test != test) != 1) {  
1013 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix != GSparseMatrix operator." << endl;  
1014 - throw;  
1015 - }  
1016 - cout << ".";  
1017 - //  
1018 - if ((m_test != bigger) != 1) {  
1019 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix != GSparseMatrix operator." << endl;  
1020 - throw;  
1021 - }  
1022 - cout << ".";  
1023 - }  
1024 - catch (exception &e) {  
1025 - cout << e.what() << endl;  
1026 - throw;  
1027 - }  
1028 - cout << " ok." << endl;  
1029 -} 863 + // Minimum
  864 + double min = m_test.min();
1030 865
  866 + // Check mimimum
  867 + double value = g_matrix[0];
  868 + for (int row = 0; row < g_rows; ++row) {
  869 + for (int col = 0; col < g_cols; ++col) {
  870 + if (g_matrix[col+row*g_cols] < value) {
  871 + value = g_matrix[col+row*g_cols];
  872 + }
  873 + }
  874 + }
  875 + test_value(min, value, 0.0, "Test minimum function");
1031 876
1032 -/***************************************************************************  
1033 - * Test: Conversion between matrix types *  
1034 - ***************************************************************************/  
1035 -void test_conversion(void)  
1036 -{  
1037 - // Dump header  
1038 - std::cout << "Test GSparseMatrix: Matrix conversions: "; 877 + // Maximum
  878 + double max = m_test.max();
1039 879
1040 - // Test 1  
1041 - try {  
1042 - //  
1043 - // Setup a symmetric sparse matrix  
1044 - int num = 10;  
1045 - GSparseMatrix sparse(num,num);  
1046 - for (int i = 0; i < num; ++i)  
1047 - sparse(i,i) = (i+1);  
1048 - for (int j = 1; j < num; ++j) {  
1049 - sparse(0,j) = (j+1);  
1050 - sparse(j,0) = (j+1);  
1051 - }  
1052 - std::cout << ".";  
1053 - //  
1054 - // Convert GSparseMatrix matrix into GSymMatrix object  
1055 -// GSymMatrix converted = sym_matrix(sparse);  
1056 - GSymMatrix converted = GSymMatrix(sparse);  
1057 - std::cout << ".";  
1058 - //  
1059 - // Convert GSymMatrix back to GSparseMatrix matrix  
1060 - GSparseMatrix back_convert = GSparseMatrix(converted);  
1061 - std::cout << ".";  
1062 - //  
1063 - // Compare back converted matrix to original one. They should be identical  
1064 - if (sparse != back_convert) {  
1065 - std::cout << std::endl  
1066 - << "TEST ERROR: Unable to convert matrixes (symmetric)."  
1067 - << std::endl;  
1068 - std::cout << "Original matrix " << sparse << std::endl;  
1069 - std::cout << "GSymMatrix matrix " << converted << std::endl;  
1070 - std::cout << "Back converted matrix " << back_convert << std::endl;  
1071 - throw;  
1072 - }  
1073 - //  
1074 - // Determine the fill of the matrix. It should be 0.28  
1075 - double fill = back_convert.fill();  
1076 - if (abs(fill-0.28) > 1.0e-15) {  
1077 - cout << endl << "TEST ERROR: Bad fill " << fill << " determined (expected 1.0)." <<  
1078 - endl;  
1079 - throw;  
1080 - }  
1081 - cout << ".";  
1082 - //  
1083 - // Extract lower triangle and check values  
1084 - #if 0  
1085 - GMatrix lower = sparse.extract_lower_triangle();  
1086 - int ok = 1;  
1087 - for (int j = 1; j < num; ++j) {  
1088 - for (int i = 0; i < j; ++i) {  
1089 - if (lower(i,j) != 0.0)  
1090 - ok = 0;  
1091 - }  
1092 - }  
1093 - for (int j = 0; j < num; ++j) {  
1094 - for (int i = j; i < num; ++i) {  
1095 - if (lower(i,j) != sparse(i,j))  
1096 - ok = 0;  
1097 - } 880 + // Check maximum
  881 + value = g_matrix[0];
  882 + for (int row = 0; row < g_rows; ++row) {
  883 + for (int col = 0; col < g_cols; ++col) {
  884 + if (g_matrix[col+row*g_cols] > value) {
  885 + value = g_matrix[col+row*g_cols];
  886 + }
  887 + }
1098 } 888 }
1099 - if (!ok) {  
1100 - cout << endl << "TEST ERROR: Corrupt extract_lower_triangle." << endl;  
1101 - cout << "Original matrix " << sparse << endl;  
1102 - cout << "Lower triangle matrix " << lower << endl;  
1103 - throw; 889 + test_value(max, value, 0.0, "Test maximum function");
  890 +
  891 + // Sum
  892 + double sum = m_test.sum();
  893 +
  894 + // Check sum
  895 + value = 0.0;
  896 + for (int row = 0; row < g_rows; ++row) {
  897 + for (int col = 0; col < g_cols; ++col) {
  898 + value += g_matrix[col+row*g_cols];
  899 + }
1104 } 900 }
1105 - cout << ".";  
1106 - //  
1107 - // Extract upper triangle and check values  
1108 - GMatrix upper = sparse.extract_upper_triangle();  
1109 - ok = 1;  
1110 - for (int j = 0; j < num; ++j) {  
1111 - for (int i = j+1; i < num; ++i) {  
1112 - if (upper(i,j) != 0.0)  
1113 - ok = 0;  
1114 - } 901 + test_value(sum, value, 1.0e-20, "Test sum function");
  902 +
  903 + // Transpose function
  904 + GSparseMatrix test1 = transpose(m_test);
  905 + test_assert(check_matrix(m_test), "Test source matrix");
  906 + test_assert(check_matrix_trans(test1, 1.0, 0.0),
  907 + "Test transpose(GSparseMatrix) function",
  908 + "Unexpected transposed matrix:\n"+test1.print());
  909 +
  910 + // Transpose method
  911 + test1 = m_test;
  912 + test1.transpose();
  913 + test_assert(check_matrix(m_test), "Test source matrix");
  914 + test_assert(check_matrix_trans(test1, 1.0, 0.0),
  915 + "Test GSparseMatrix.transpose() method",
  916 + "Unexpected transposed matrix:\n"+test1.print());
  917 +
  918 + // Convert to general matrix
  919 + GSparseMatrix test2 = GSparseMatrix(m_test);
  920 + test_assert(check_matrix(m_test), "Test source matrix");
  921 + test_assert(check_matrix(test2, 1.0, 0.0),
  922 + "Test GSparseMatrix(GSparseMatrix) constructor",
  923 + "Unexpected GSparseMatrix:\n"+test2.print());
  924 +
  925 + // Try extracting of lower triangle
  926 + /*
  927 + test_try("Try extracting of lower triangle");
  928 + try {
  929 + test2 = m_test.extract_lower_triangle();
  930 + test_try_failure("Expected GException::matrix_not_square exception.");
1115 } 931 }
1116 - for (int j = 0; j < num; ++j) {  
1117 - for (int i = 0; i <= j; ++i) {  
1118 - if (upper(i,j) != sparse(i,j))  
1119 - ok = 0;  
1120 - } 932 + catch (GException::matrix_not_square &e) {
  933 + test_try_success();
1121 } 934 }
1122 - if (!ok) {  
1123 - cout << endl << "TEST ERROR: Corrupt extract_upper_triangle." << endl;  
1124 - cout << "Original matrix " << sparse << endl;  
1125 - cout << "Upper triangle matrix " << upper << endl;  
1126 - throw; 935 + catch (std::exception &e) {
  936 + test_try_failure(e);
1127 } 937 }
1128 - cout << ".";  
1129 - #endif  
1130 - //  
1131 - // Now make the matrix unsymmetric  
1132 - sparse(0,num-1) = 1000.0;  
1133 - //  
1134 - // Try converting now into GSymMatrix object (this should fail) 938 + */
  939 +
  940 + // Try extracting of upper triangle
  941 + /*
  942 + test_try("Try extracting of lower triangle");
1135 try { 943 try {
1136 - converted = GSymMatrix(sparse);  
1137 - cout << endl << "TEST ERROR: Sparse->Sym conversion should have failed." << endl;  
1138 - throw;  
1139 - }  
1140 - catch (GException::matrix_not_symmetric &e) {  
1141 - cout << "."; 944 + test2 = m_test.extract_upper_triangle();
  945 + test_try_failure("Expected GException::matrix_not_square exception.");
1142 } 946 }
1143 - catch (exception &e) {  
1144 - cout << e.what() << endl;  
1145 - throw; 947 + catch (GException::matrix_not_square &e) {
  948 + test_try_success();
1146 } 949 }
1147 - //  
1148 - // Convert matrix into GMatrix object  
1149 - GMatrix full = GMatrix(sparse);  
1150 - cout << ".";  
1151 - //  
1152 - // Convert full matrix back to GSparseMatrix  
1153 - back_convert = GSparseMatrix(full);  
1154 - cout << ".";  
1155 - //  
1156 - // Compare back converted matrix to original one. They should be identical  
1157 - if (sparse != back_convert) {  
1158 - cout << endl << "TEST ERROR: Unable to convert matrixes (sparse)." << endl;  
1159 - cout << "Original matrix " << sparse << endl;  
1160 - cout << "GMatrix matrix " << full << endl;  
1161 - cout << "Back converted matrix " << back_convert << endl;  
1162 - throw; 950 + catch (std::exception &e) {
  951 + test_try_failure(e);
1163 } 952 }
1164 - cout << ".";  
1165 - }  
1166 - catch (exception &e) {  
1167 - cout << e.what() << endl;  
1168 - throw;  
1169 - }  
1170 - cout << " ok." << endl;  
1171 -  
1172 - // Return  
1173 - return; 953 + */
  954 +
  955 + // Extract lower triangle for square matrix
  956 + /*
  957 + GSparseMatrix test3 = m_test * transpose(m_test);
  958 + GSparseMatrix test4 = test3.extract_lower_triangle();
  959 + test_assert(check_matrix_lt(test4, test3),
  960 + "Test GSparseMatrix.extract_lower_triangle() method",
  961 + "Unexpected lower triangle:\n"+test4.print());
  962 +
  963 + // Extract upper triangle
  964 + test4 = test3.extract_upper_triangle();
  965 + test_assert(check_matrix_ut(test4, test3),
  966 + "Test GSparseMatrix.extract_upper_triangle() method",
  967 + "Unexpected lower triangle:\n"+test4.print());
  968 + */
  969 +
  970 + // Return
  971 + return;
1174 } 972 }
1175 973
1176 974
1177 -/***************************************************************************  
1178 - * Test: Cholesky factorisation * 975 +/***********************************************************************//**
  976 + * @brief Test matrix comparisons
1179 ***************************************************************************/ 977 ***************************************************************************/
1180 -void test_cholesky(void) 978 +void TestGSparseMatrix::matrix_compare(void)
1181 { 979 {
1182 - // Dump header  
1183 - std::cout << "Test GSparseMatrix: Cholesky decomposition, solver and inverter: "; 980 + // Allocate an empty matrix
  981 + GSparseMatrix empty(g_rows, g_cols);
  982 +
  983 + // Allocate a full matrix
  984 + GSparseMatrix full(g_rows, g_cols);
  985 + for (int i = 0; i < g_rows; ++i) {
  986 + for (int j = 0; j < g_rows; ++j) {
  987 + full(i,j) = 1.0 + i*1000.0 + j*10.0;
  988 + }
  989 + }
  990 +
  991 + // Test == operator
  992 + test_assert((m_test == m_test), "Test == operator",
  993 + "Identity test faileded on:\n"+m_test.print());
  994 + test_assert(!(m_test == empty), "Test == operator",
  995 + "Matrix\n"+m_test.print()+"\nshould not be equal to\n"+empty.print());
  996 + test_assert(!(m_test == full), "Test == operator",
  997 + "Matrix\n"+m_test.print()+"\nshould not be equal to\n"+full.print());
  998 + test_assert(!(m_test == m_bigger), "Test == operator",
  999 + "Matrix\n"+m_test.print()+"\nshould not be equal to\n"+m_bigger.print());
  1000 +
  1001 + // Test != operator
  1002 + test_assert(!(m_test != m_test), "Test != operator",
  1003 + "Unequality test faileded on:\n"+m_test.print());
  1004 + test_assert((m_test != empty), "Test != operator",
  1005 + "Matrix\n"+m_test.print()+"\nshould differ from\n"+empty.print());
  1006 + test_assert((m_test != full), "Test != operator",
  1007 + "Matrix\n"+m_test.print()+"\nshould differ from\n"+full.print());
  1008 + test_assert((m_test != m_bigger), "Test != operator",
  1009 + "Matrix\n"+m_test.print()+"\nshould differ from\n"+m_bigger.print());
  1010 +
  1011 + // Return
  1012 + return;
  1013 +}
1184 1014
1185 - // Test 1 1015 +
  1016 +/***********************************************************************//**
  1017 + * @brief Test Cholesky decomposition
  1018 + ***************************************************************************/
  1019 +void TestGSparseMatrix::matrix_cholesky(void)
  1020 +{
  1021 + // Setup matrix for Cholesky decomposition
  1022 + GSparseMatrix chol_test(5,5);
  1023 + chol_test(0,0) = 1.0;
  1024 + chol_test(0,1) = 0.2;
  1025 + chol_test(0,2) = 0.2;
  1026 + chol_test(0,3) = 0.2;
  1027 + chol_test(0,4) = 0.2;
  1028 + chol_test(1,0) = 0.2;
  1029 + chol_test(2,0) = 0.2;
  1030 + chol_test(3,0) = 0.2;
  1031 + chol_test(4,0) = 0.2;
  1032 + chol_test(1,1) = 1.0;
  1033 + chol_test(2,2) = 1.0;
  1034 + chol_test(3,3) = 1.0;
  1035 + chol_test(4,4) = 1.0;
  1036 +
  1037 + // Try to solve now (should not work)
  1038 + test_try("Try Cholesky solver without factorisation");
1186 try { 1039 try {
1187 - //  
1188 - // Setup matrix for Cholesky decomposition  
1189 - GSparseMatrix m_chol_test(5,5);  
1190 - m_chol_test(0,0) = 1.0;  
1191 - m_chol_test(0,1) = 0.2;  
1192 - m_chol_test(0,2) = 0.2;  
1193 - m_chol_test(0,3) = 0.2;  
1194 - m_chol_test(0,4) = 0.2;  
1195 - m_chol_test(1,0) = 0.2;  
1196 - m_chol_test(2,0) = 0.2;  
1197 - m_chol_test(3,0) = 0.2;  
1198 - m_chol_test(4,0) = 0.2;  
1199 - m_chol_test(1,1) = 1.0;  
1200 - m_chol_test(2,2) = 1.0;  
1201 - m_chol_test(3,3) = 1.0;  
1202 - m_chol_test(4,4) = 1.0;  
1203 - //  
1204 - // Try to solve now (should not work)  
1205 - try {  
1206 - GVector v_test12(5);  
1207 - v_test12 = m_chol_test.cholesky_solver(v_test12);  
1208 - }  
1209 - catch (GException::matrix_not_factorised) {  
1210 - std::cout << ".";  
1211 - }  
1212 - catch (exception &e) {  
1213 - std::cout << e.what() << std::endl;  
1214 - throw;  
1215 - }  
1216 - //  
1217 - // Test Cholesky decomposition  
1218 - GSparseMatrix cd = cholesky_decompose(m_chol_test);  
1219 - std::cout << ".";  
1220 - //  
1221 - // Test inplace Cholesky decomposition  
1222 - cd = m_chol_test;  
1223 - cd.cholesky_decompose();  
1224 - std::cout << ".";  
1225 - //  
1226 - // Test Cholesky solver  
1227 - GVector e0(5);  
1228 - GVector a0(5);  
1229 - e0[0] = 1.0;  
1230 - a0[0] = 1.0;  
1231 - a0[1] = 0.2;  
1232 - a0[2] = 0.2;  
1233 - a0[3] = 0.2;  
1234 - a0[4] = 0.2;  
1235 - GVector s0 = cd.cholesky_solver(a0);  
1236 - double res = max(abs(s0-e0));  
1237 - if (res < 1.0e-15) {  
1238 - #if defined(DUMP_RESIDUALS)  
1239 - std::cout << " Res(S0)=" << res << " ";  
1240 - #else  
1241 - std::cout << ".";  
1242 - #endif  
1243 - }  
1244 - else {  
1245 - std::cout << std::endl  
1246 - << "TEST ERROR: Corrupt cholesky_solver(GVector) function."  
1247 - << std::endl;  
1248 - std::cout << "Residual vector (all elements should be zero):" << std::endl;  
1249 - std::cout << s0 << std::endl;  
1250 - throw;  
1251 - }  
1252 - // 1040 + GVector vector(5);
  1041 + vector = chol_test.cholesky_solver(vector);
  1042 + test_try_failure("Expected GException::matrix_not_factorised exception.");
  1043 + }
  1044 + catch (GException::matrix_not_factorised) {
  1045 + test_try_success();
  1046 + }
  1047 + catch (std::exception &e) {
  1048 + test_try_failure(e);
  1049 + }
  1050 +
  1051 + // Perform Cholesky decomposition
  1052 + GSparseMatrix cd = cholesky_decompose(chol_test);
  1053 +
  1054 + // Perform inplace Cholesky decomposition
  1055 + cd = chol_test;
  1056 + cd.cholesky_decompose();
  1057 +
  1058 + // Test Cholesky solver (first test)
  1059 + GVector e0(5);
  1060 + GVector a0(5);
  1061 + e0[0] = 1.0;
  1062 + a0[0] = 1.0;
  1063 + a0[1] = 0.2;
  1064 + a0[2] = 0.2;
  1065 + a0[3] = 0.2;
  1066 + a0[4] = 0.2;
  1067 + GVector s0 = cd.cholesky_solver(a0);
  1068 + double res = max(abs(s0-e0));
  1069 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method - 1");
  1070 +
  1071 + // Test Cholesky solver (second test)
1253 e0[0] = 0.0; 1072 e0[0] = 0.0;
1254 e0[1] = 1.0; 1073 e0[1] = 1.0;
1255 a0[0] = 0.2; 1074 a0[0] = 0.2;
@@ -1259,97 +1078,56 @@ void test_cholesky(void) @@ -1259,97 +1078,56 @@ void test_cholesky(void)
1259 a0[4] = 0.0; 1078 a0[4] = 0.0;
1260 s0 = cd.cholesky_solver(a0); 1079 s0 = cd.cholesky_solver(a0);
1261 res = max(abs(s0-e0)); 1080 res = max(abs(s0-e0));
1262 - if (res < 1.0e-15)  
1263 - #if defined(DUMP_RESIDUALS)  
1264 - cout << "Res(S1)=" << res << " ";  
1265 - #else  
1266 - cout << ".";  
1267 - #endif  
1268 - else {  
1269 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1270 - cout << "Residual vector (all elements should be zero):" << endl;  
1271 - cout << s0 << endl;  
1272 - throw;  
1273 - }  
1274 - // 1081 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method - 2");
  1082 +
  1083 + // Test Cholesky solver (third test)
1275 e0[1] = 0.0; 1084 e0[1] = 0.0;
1276 e0[2] = 1.0; 1085 e0[2] = 1.0;
1277 a0[1] = 0.0; 1086 a0[1] = 0.0;
1278 a0[2] = 1.0; 1087 a0[2] = 1.0;
1279 s0 = cd.cholesky_solver(a0); 1088 s0 = cd.cholesky_solver(a0);
1280 res = max(abs(s0-e0)); 1089 res = max(abs(s0-e0));
1281 - if (res < 1.0e-15)  
1282 - #if defined(DUMP_RESIDUALS)  
1283 - cout << "Res(S2)=" << res << " ";  
1284 - #else  
1285 - cout << ".";  
1286 - #endif  
1287 - else {  
1288 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1289 - cout << "Residual vector (all elements should be zero):" << endl;  
1290 - cout << s0 << endl;  
1291 - throw;  
1292 - }  
1293 - // 1090 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method - 3");
  1091 +
  1092 + // Test Cholesky solver (forth test)
1294 e0[2] = 0.0; 1093 e0[2] = 0.0;
1295 e0[3] = 1.0; 1094 e0[3] = 1.0;
1296 a0[2] = 0.0; 1095 a0[2] = 0.0;
1297 a0[3] = 1.0; 1096 a0[3] = 1.0;
1298 s0 = cd.cholesky_solver(a0); 1097 s0 = cd.cholesky_solver(a0);
1299 res = max(abs(s0-e0)); 1098 res = max(abs(s0-e0));
1300 - if (res < 1.0e-15)  
1301 - #if defined(DUMP_RESIDUALS)  
1302 - cout << "Res(S3)=" << res << " ";  
1303 - #else  
1304 - cout << ".";  
1305 - #endif  
1306 - else {  
1307 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1308 - cout << "Residual vector (all elements should be zero):" << endl;  
1309 - cout << s0 << endl;  
1310 - throw;  
1311 - }  
1312 - // 1099 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method - 4");
  1100 +
  1101 + // Test Cholesky solver (fifth test)
1313 e0[3] = 0.0; 1102 e0[3] = 0.0;
1314 e0[4] = 1.0; 1103 e0[4] = 1.0;
1315 a0[3] = 0.0; 1104 a0[3] = 0.0;
1316 a0[4] = 1.0; 1105 a0[4] = 1.0;
1317 s0 = cd.cholesky_solver(a0); 1106 s0 = cd.cholesky_solver(a0);
1318 res = max(abs(s0-e0)); 1107 res = max(abs(s0-e0));
1319 - if (res < 1.0e-15)  
1320 - #if defined(DUMP_RESIDUALS)  
1321 - cout << "Res(S4)=" << res << " ";  
1322 - #else  
1323 - cout << ".";  
1324 - #endif  
1325 - else {  
1326 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1327 - cout << "Residual vector (all elements should be zero):" << endl;  
1328 - cout << s0 << endl;  
1329 - throw;  
1330 - }  
1331 - // 1108 + test_value(res, 0.0, 1.0e-15, "Test cholesky_solver() method - 5");
  1109 +
1332 // Setup matrix for Cholesky decomposition with zero row/col 1110 // Setup matrix for Cholesky decomposition with zero row/col
1333 - GSparseMatrix m_chol_test_zero(6,6);  
1334 - m_chol_test_zero(0,0) = 1.0;  
1335 - m_chol_test_zero(0,1) = 0.2;  
1336 - m_chol_test_zero(0,2) = 0.2;  
1337 - m_chol_test_zero(0,4) = 0.2;  
1338 - m_chol_test_zero(0,5) = 0.2;  
1339 - m_chol_test_zero(1,0) = 0.2;  
1340 - m_chol_test_zero(2,0) = 0.2;  
1341 - m_chol_test_zero(4,0) = 0.2;  
1342 - m_chol_test_zero(5,0) = 0.2;  
1343 - m_chol_test_zero(1,1) = 1.0;  
1344 - m_chol_test_zero(2,2) = 1.0;  
1345 - m_chol_test_zero(4,4) = 1.0;  
1346 - m_chol_test_zero(5,5) = 1.0;  
1347 - // 1111 + GSparseMatrix chol_test_zero(6,6);
  1112 + chol_test_zero(0,0) = 1.0;
  1113 + chol_test_zero(0,1) = 0.2;
  1114 + chol_test_zero(0,2) = 0.2;
  1115 + chol_test_zero(0,4) = 0.2;
  1116 + chol_test_zero(0,5) = 0.2;
  1117 + chol_test_zero(1,0) = 0.2;
  1118 + chol_test_zero(2,0) = 0.2;
  1119 + chol_test_zero(4,0) = 0.2;
  1120 + chol_test_zero(5,0) = 0.2;
  1121 + chol_test_zero(1,1) = 1.0;
  1122 + chol_test_zero(2,2) = 1.0;
  1123 + chol_test_zero(4,4) = 1.0;
  1124 + chol_test_zero(5,5) = 1.0;
  1125 +
1348 // Test compressed Cholesky decomposition 1126 // Test compressed Cholesky decomposition
1349 - GSparseMatrix cd_zero = m_chol_test_zero; 1127 + GSparseMatrix cd_zero = chol_test_zero;
1350 cd_zero.cholesky_decompose(); 1128 cd_zero.cholesky_decompose();
1351 - //  
1352 - // Test compressed Cholesky solver 1129 +
  1130 + // Test compressed Cholesky solver (first test)
1353 e0 = GVector(6); 1131 e0 = GVector(6);
1354 a0 = GVector(6); 1132 a0 = GVector(6);
1355 e0[0] = 1.0; 1133 e0[0] = 1.0;
@@ -1360,19 +1138,9 @@ void test_cholesky(void) @@ -1360,19 +1138,9 @@ void test_cholesky(void)
1360 a0[5] = 0.2; 1138 a0[5] = 0.2;
1361 s0 = cd_zero.cholesky_solver(a0); 1139 s0 = cd_zero.cholesky_solver(a0);
1362 res = max(abs(s0-e0)); 1140 res = max(abs(s0-e0));
1363 - if (res < 1.0e-15)  
1364 - #if defined(DUMP_RESIDUALS)  
1365 - cout << "Res(S0Z)=" << res << " ";  
1366 - #else  
1367 - cout << ".";  
1368 - #endif  
1369 - else {  
1370 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1371 - cout << "Residual vector (all elements should be zero):" << endl;  
1372 - cout << s0 << endl;  
1373 - throw;  
1374 - }  
1375 - // 1141 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method - 1");
  1142 +
  1143 + // Test compressed Cholesky solver (second test)
1376 e0[0] = 0.0; 1144 e0[0] = 0.0;
1377 e0[1] = 1.0; 1145 e0[1] = 1.0;
1378 a0[0] = 0.2; 1146 a0[0] = 0.2;
@@ -1382,96 +1150,55 @@ void test_cholesky(void) @@ -1382,96 +1150,55 @@ void test_cholesky(void)
1382 a0[5] = 0.0; 1150 a0[5] = 0.0;
1383 s0 = cd_zero.cholesky_solver(a0); 1151 s0 = cd_zero.cholesky_solver(a0);
1384 res = max(abs(s0-e0)); 1152 res = max(abs(s0-e0));
1385 - if (res < 1.0e-15)  
1386 - #if defined(DUMP_RESIDUALS)  
1387 - cout << "Res(S1Z)=" << res << " ";  
1388 - #else  
1389 - cout << ".";  
1390 - #endif  
1391 - else {  
1392 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1393 - cout << "Residual vector (all elements should be zero):" << endl;  
1394 - cout << s0 << endl;  
1395 - throw;  
1396 - }  
1397 - // 1153 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method - 2");
  1154 +
  1155 + // Test compressed Cholesky solver (third test)
1398 e0[1] = 0.0; 1156 e0[1] = 0.0;
1399 e0[2] = 1.0; 1157 e0[2] = 1.0;
1400 a0[1] = 0.0; 1158 a0[1] = 0.0;
1401 a0[2] = 1.0; 1159 a0[2] = 1.0;
1402 s0 = cd_zero.cholesky_solver(a0); 1160 s0 = cd_zero.cholesky_solver(a0);
1403 res = max(abs(s0-e0)); 1161 res = max(abs(s0-e0));
1404 - if (res < 1.0e-15)  
1405 - #if defined(DUMP_RESIDUALS)  
1406 - cout << "Res(S2Z)=" << res << " ";  
1407 - #else  
1408 - cout << ".";  
1409 - #endif  
1410 - else {  
1411 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1412 - cout << "Residual vector (all elements should be zero):" << endl;  
1413 - cout << s0 << endl;  
1414 - throw;  
1415 - }  
1416 - // 1162 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method - 3");
  1163 +
  1164 + // Test compressed Cholesky solver (forth test)
1417 e0[2] = 0.0; 1165 e0[2] = 0.0;
1418 e0[4] = 1.0; 1166 e0[4] = 1.0;
1419 a0[2] = 0.0; 1167 a0[2] = 0.0;
1420 a0[4] = 1.0; 1168 a0[4] = 1.0;
1421 s0 = cd_zero.cholesky_solver(a0); 1169 s0 = cd_zero.cholesky_solver(a0);
1422 res = max(abs(s0-e0)); 1170 res = max(abs(s0-e0));
1423 - if (res < 1.0e-15)  
1424 - #if defined(DUMP_RESIDUALS)  
1425 - cout << "Res(S3Z)=" << res << " ";  
1426 - #else  
1427 - cout << ".";  
1428 - #endif  
1429 - else {  
1430 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1431 - cout << "Residual vector (all elements should be zero):" << endl;  
1432 - cout << s0 << endl;  
1433 - throw;  
1434 - }  
1435 - // 1171 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method - 4");
  1172 +
  1173 + // Test compressed Cholesky solver (fifth test)
1436 e0[4] = 0.0; 1174 e0[4] = 0.0;
1437 e0[5] = 1.0; 1175 e0[5] = 1.0;
1438 a0[4] = 0.0; 1176 a0[4] = 0.0;
1439 a0[5] = 1.0; 1177 a0[5] = 1.0;
1440 s0 = cd_zero.cholesky_solver(a0); 1178 s0 = cd_zero.cholesky_solver(a0);
1441 res = max(abs(s0-e0)); 1179 res = max(abs(s0-e0));
1442 - if (res < 1.0e-15)  
1443 - #if defined(DUMP_RESIDUALS)  
1444 - cout << "Res(S4Z)=" << res << " ";  
1445 - #else  
1446 - cout << ".";  
1447 - #endif  
1448 - else {  
1449 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1450 - cout << "Residual vector (all elements should be zero):" << endl;  
1451 - cout << s0 << endl;  
1452 - throw;  
1453 - }  
1454 - // 1180 + test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method - 5");
  1181 +
1455 // Setup matrix for Cholesky decomposition with zero row/col (unsymmetric case) 1182 // Setup matrix for Cholesky decomposition with zero row/col (unsymmetric case)
1456 - GSparseMatrix m_chol_test_zero2(6,5);  
1457 - m_chol_test_zero2(0,0) = 1.0;  
1458 - m_chol_test_zero2(0,1) = 0.2;  
1459 - m_chol_test_zero2(0,2) = 0.2;  
1460 - m_chol_test_zero2(0,3) = 0.2;  
1461 - m_chol_test_zero2(0,4) = 0.2;  
1462 - m_chol_test_zero2(1,0) = 0.2;  
1463 - m_chol_test_zero2(2,0) = 0.2;  
1464 - m_chol_test_zero2(4,0) = 0.2;  
1465 - m_chol_test_zero2(5,0) = 0.2;  
1466 - m_chol_test_zero2(1,1) = 1.0;  
1467 - m_chol_test_zero2(2,2) = 1.0;  
1468 - m_chol_test_zero2(4,3) = 1.0;  
1469 - m_chol_test_zero2(5,4) = 1.0;  
1470 - // 1183 + GSparseMatrix chol_test_zero2(6,5);
  1184 + chol_test_zero2(0,0) = 1.0;
  1185 + chol_test_zero2(0,1) = 0.2;
  1186 + chol_test_zero2(0,2) = 0.2;
  1187 + chol_test_zero2(0,3) = 0.2;
  1188 + chol_test_zero2(0,4) = 0.2;
  1189 + chol_test_zero2(1,0) = 0.2;
  1190 + chol_test_zero2(2,0) = 0.2;
  1191 + chol_test_zero2(4,0) = 0.2;
  1192 + chol_test_zero2(5,0) = 0.2;
  1193 + chol_test_zero2(1,1) = 1.0;
  1194 + chol_test_zero2(2,2) = 1.0;
  1195 + chol_test_zero2(4,3) = 1.0;
  1196 + chol_test_zero2(5,4) = 1.0;
  1197 +
1471 // Test compressed Cholesky decomposition (unsymmetric case) 1198 // Test compressed Cholesky decomposition (unsymmetric case)
1472 - GSparseMatrix cd_zero2 = m_chol_test_zero2; 1199 + GSparseMatrix cd_zero2 = chol_test_zero2;
1473 cd_zero2.cholesky_decompose(); 1200 cd_zero2.cholesky_decompose();
1474 - // 1201 +
1475 // Test compressed Cholesky solver (unsymmetric case) 1202 // Test compressed Cholesky solver (unsymmetric case)
1476 e0 = GVector(5); 1203 e0 = GVector(5);
1477 a0 = GVector(6); 1204 a0 = GVector(6);
@@ -1483,19 +1210,10 @@ void test_cholesky(void) @@ -1483,19 +1210,10 @@ void test_cholesky(void)
1483 a0[5] = 0.2; 1210 a0[5] = 0.2;
1484 s0 = cd_zero2.cholesky_solver(a0); 1211 s0 = cd_zero2.cholesky_solver(a0);
1485 res = max(abs(s0-e0)); 1212 res = max(abs(s0-e0));
1486 - if (res < 1.0e-15)  
1487 - #if defined(DUMP_RESIDUALS)  
1488 - cout << "Res(S0ZU)=" << res << " ";  
1489 - #else  
1490 - cout << ".";  
1491 - #endif  
1492 - else {  
1493 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1494 - cout << "Residual vector (all elements should be zero):" << endl;  
1495 - cout << s0 << endl;  
1496 - throw;  
1497 - }  
1498 - // 1213 + test_value(res, 0.0, 1.0e-15,
  1214 + "Test unsymmetric compressed cholesky_solver() method - 1");
  1215 +
  1216 + // Test compressed Cholesky solver (unsymmetric case)
1499 e0[0] = 0.0; 1217 e0[0] = 0.0;
1500 e0[1] = 1.0; 1218 e0[1] = 1.0;
1501 a0[0] = 0.2; 1219 a0[0] = 0.2;
@@ -1505,76 +1223,39 @@ void test_cholesky(void) @@ -1505,76 +1223,39 @@ void test_cholesky(void)
1505 a0[5] = 0.0; 1223 a0[5] = 0.0;
1506 s0 = cd_zero2.cholesky_solver(a0); 1224 s0 = cd_zero2.cholesky_solver(a0);
1507 res = max(abs(s0-e0)); 1225 res = max(abs(s0-e0));
1508 - if (res < 1.0e-15)  
1509 - #if defined(DUMP_RESIDUALS)  
1510 - cout << "Res(S1ZU)=" << res << " ";  
1511 - #else  
1512 - cout << ".";  
1513 - #endif  
1514 - else {  
1515 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1516 - cout << "Residual vector (all elements should be zero):" << endl;  
1517 - cout << s0 << endl;  
1518 - throw;  
1519 - }  
1520 - // 1226 + test_value(res, 0.0, 1.0e-15,
  1227 + "Test unsymmetric compressed cholesky_solver() method - 2");
  1228 +
  1229 + // Test compressed Cholesky solver (unsymmetric case)
1521 e0[1] = 0.0; 1230 e0[1] = 0.0;
1522 e0[2] = 1.0; 1231 e0[2] = 1.0;
1523 a0[1] = 0.0; 1232 a0[1] = 0.0;
1524 a0[2] = 1.0; 1233 a0[2] = 1.0;
1525 s0 = cd_zero2.cholesky_solver(a0); 1234 s0 = cd_zero2.cholesky_solver(a0);
1526 res = max(abs(s0-e0)); 1235 res = max(abs(s0-e0));
1527 - if (res < 1.0e-15)  
1528 - #if defined(DUMP_RESIDUALS)  
1529 - cout << "Res(S2ZU)=" << res << " ";  
1530 - #else  
1531 - cout << ".";  
1532 - #endif  
1533 - else {  
1534 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1535 - cout << "Residual vector (all elements should be zero):" << endl;  
1536 - cout << s0 << endl;  
1537 - throw;  
1538 - }  
1539 - // 1236 + test_value(res, 0.0, 1.0e-15,
  1237 + "Test unsymmetric compressed cholesky_solver() method - 3");
  1238 +
  1239 + // Test compressed Cholesky solver (unsymmetric case)
1540 e0[2] = 0.0; 1240 e0[2] = 0.0;
1541 e0[3] = 1.0; 1241 e0[3] = 1.0;
1542 a0[2] = 0.0; 1242 a0[2] = 0.0;
1543 a0[4] = 1.0; 1243 a0[4] = 1.0;
1544 s0 = cd_zero2.cholesky_solver(a0); 1244 s0 = cd_zero2.cholesky_solver(a0);
1545 res = max(abs(s0-e0)); 1245 res = max(abs(s0-e0));
1546 - if (res < 1.0e-15)  
1547 - #if defined(DUMP_RESIDUALS)  
1548 - cout << "Res(S3ZU)=" << res << " ";  
1549 - #else  
1550 - cout << ".";  
1551 - #endif  
1552 - else {  
1553 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1554 - cout << "Residual vector (all elements should be zero):" << endl;  
1555 - cout << s0 << endl;  
1556 - throw;  
1557 - }  
1558 - // 1246 + test_value(res, 0.0, 1.0e-15,
  1247 + "Test unsymmetric compressed cholesky_solver() method - 4");
  1248 +
  1249 + // Test compressed Cholesky solver (unsymmetric case)
1559 e0[3] = 0.0; 1250 e0[3] = 0.0;
1560 e0[4] = 1.0; 1251 e0[4] = 1.0;
1561 a0[4] = 0.0; 1252 a0[4] = 0.0;
1562 a0[5] = 1.0; 1253 a0[5] = 1.0;
1563 s0 = cd_zero2.cholesky_solver(a0); 1254 s0 = cd_zero2.cholesky_solver(a0);
1564 res = max(abs(s0-e0)); 1255 res = max(abs(s0-e0));
1565 - if (res < 1.0e-15)  
1566 - #if defined(DUMP_RESIDUALS)  
1567 - cout << "Res(S4ZU)=" << res << " ";  
1568 - #else  
1569 - cout << ".";  
1570 - #endif  
1571 - else {  
1572 - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;  
1573 - cout << "Residual vector (all elements should be zero):" << endl;  
1574 - cout << s0 << endl;  
1575 - throw;  
1576 - }  
1577 - // 1256 + test_value(res, 0.0, 1.0e-15,
  1257 + "Test unsymmetric compressed cholesky_solver() method - 5");
  1258 +
1578 // Test Cholesky inverter (inplace) 1259 // Test Cholesky inverter (inplace)
1579 GSparseMatrix unit(5,5); 1260 GSparseMatrix unit(5,5);
1580 unit(0,0) = 1.0; 1261 unit(0,0) = 1.0;
@@ -1582,77 +1263,33 @@ void test_cholesky(void) @@ -1582,77 +1263,33 @@ void test_cholesky(void)
1582 unit(2,2) = 1.0; 1263 unit(2,2) = 1.0;
1583 unit(3,3) = 1.0; 1264 unit(3,3) = 1.0;
1584 unit(4,4) = 1.0; 1265 unit(4,4) = 1.0;
1585 - GSparseMatrix m_chol_test_inv = m_chol_test;  
1586 - m_chol_test_inv.cholesky_invert();  
1587 - GSparseMatrix ci_product = m_chol_test * m_chol_test_inv; 1266 + GSparseMatrix chol_test_inv = chol_test;
  1267 + chol_test_inv.cholesky_invert();
  1268 + GSparseMatrix ci_product = chol_test * chol_test_inv;
1588 GSparseMatrix ci_residuals = ci_product - unit; 1269 GSparseMatrix ci_residuals = ci_product - unit;
1589 res = (abs(ci_residuals)).max(); 1270 res = (abs(ci_residuals)).max();
1590 - if (res < 1.0e-15)  
1591 - #if defined(DUMP_RESIDUALS)  
1592 - cout << "Res(CI)=" << res << " ";  
1593 - #else  
1594 - cout << ".";  
1595 - #endif  
1596 - else {  
1597 - cout << endl << "TEST ERROR: Corrupt GSparseMatrix.cholesky_invert() function." << endl;  
1598 - cout << "Residual matrix (all elements should be zero):" << endl;  
1599 - cout << ci_residuals << endl;  
1600 - throw;  
1601 - }  
1602 - // 1271 + test_value(res, 0.0, 1.0e-15, "Test inplace Cholesky inverter");
  1272 +
1603 // Test Cholesky inverter 1273 // Test Cholesky inverter
1604 - m_chol_test_inv = cholesky_invert(m_chol_test);  
1605 - ci_product = m_chol_test * m_chol_test_inv;  
1606 - ci_residuals = ci_product - unit;  
1607 - res = (abs(ci_residuals)).max();  
1608 - if (res < 1.0e-15)  
1609 - #if defined(DUMP_RESIDUALS)  
1610 - cout << "Res(CI2)=" << res << " ";  
1611 - #else  
1612 - cout << ".";  
1613 - #endif  
1614 - else {  
1615 - cout << endl << "TEST ERROR: Corrupt cholesky_invert(GSparseMatrix) function." << endl;  
1616 - cout << "Residual matrix (all elements should be zero):" << endl;  
1617 - cout << ci_residuals << endl;  
1618 - throw;  
1619 - }  
1620 - //  
1621 - // Test Cholesky inverter for compressed matrix  
1622 - unit = GSparseMatrix(6,6);  
1623 - unit(0,0) = 1.0;  
1624 - unit(1,1) = 1.0;  
1625 - unit(2,2) = 1.0;  
1626 - unit(4,4) = 1.0;  
1627 - unit(5,5) = 1.0;  
1628 - GSparseMatrix m_chol_test_zero_inv = m_chol_test_zero;  
1629 - m_chol_test_zero_inv.cholesky_invert();  
1630 - GSparseMatrix ciz_product = m_chol_test_zero * m_chol_test_zero_inv;  
1631 - GSparseMatrix ciz_residuals = ciz_product - unit;  
1632 - res = (abs(ciz_residuals)).max();  
1633 - if (res < 1.0e-15) {  
1634 - #if defined(DUMP_RESIDUALS)  
1635 - std::cout << "Res(CIZ)=" << res << " ";  
1636 - #else  
1637 - std::cout << ".";  
1638 - #endif  
1639 - }  
1640 - else {  
1641 - std::cout << std::endl  
1642 - << "TEST ERROR: Corrupt compressed GSparseMatrix.cholesky_invert()"  
1643 - " method." << std::endl;  
1644 - std::cout << "Residual matrix (all elements should be zero):" << std::endl;  
1645 - std::cout << ciz_residuals << std::endl;  
1646 - throw;  
1647 - }  
1648 - }  
1649 - catch (exception &e) {  
1650 - std::cout << e.what() << std::endl;  
1651 - throw;  
1652 - }  
1653 -  
1654 - // Dump final test ok  
1655 - std::cout << "ok." << std::endl; 1274 + chol_test_inv = cholesky_invert(chol_test);
  1275 + ci_product = chol_test * chol_test_inv;
  1276 + ci_residuals = ci_product - unit;
  1277 + res = (abs(ci_residuals)).max();
  1278 + test_value(res, 0.0, 1.0e-15, "Test Cholesky inverter");
  1279 +
  1280 + // Test Cholesky inverter for compressed matrix
  1281 + unit = GSparseMatrix(6,6);
  1282 + unit(0,0) = 1.0;
  1283 + unit(1,1) = 1.0;
  1284 + unit(2,2) = 1.0;
  1285 + unit(4,4) = 1.0;
  1286 + unit(5,5) = 1.0;
  1287 + GSparseMatrix chol_test_zero_inv = chol_test_zero;
  1288 + chol_test_zero_inv.cholesky_invert();
  1289 + GSparseMatrix ciz_product = chol_test_zero * chol_test_zero_inv;
  1290 + GSparseMatrix ciz_residuals = ciz_product - unit;
  1291 + res = (abs(ciz_residuals)).max();
  1292 + test_value(res, 0.0, 1.0e-15, "Test compressed matrix Cholesky inverter");
1656 1293
1657 // Return 1294 // Return
1658 return; 1295 return;
@@ -1660,227 +1297,97 @@ void test_cholesky(void) @@ -1660,227 +1297,97 @@ void test_cholesky(void)
1660 1297
1661 1298
1662 /*************************************************************************** 1299 /***************************************************************************
1663 - * Test: Heavy calculus * 1300 + * @brief Test matrix printing
1664 ***************************************************************************/ 1301 ***************************************************************************/
1665 -void test_heavy(int block = 512) 1302 +void TestGSparseMatrix::matrix_print(void)
1666 { 1303 {
1667 - // Allocate some variables  
1668 - clock_t t_start;  
1669 - double t_elapse;  
1670 - int number;  
1671 - int num_cols;  
1672 - int i_min, i_max;  
1673 -  
1674 - // Perform test  
1675 - cout << "Test GSparseMatrix: Heavy calculus: ";  
1676 - try {  
1677 - //  
1678 - // Fill 10000 x 10000 matrix with 10000 values  
1679 - t_start = clock();  
1680 - number = 10000;  
1681 - GSparseMatrix large(number,number);  
1682 - large.set_mem_block(block);  
1683 - for (int i = 0; i < number; ++i)  
1684 - large(i,i) = (i+1)*3.14;  
1685 - t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);  
1686 - #if defined(DUMP_TIMING)  
1687 - cout << endl << " - Fill of 10000 values needed " << t_elapse <<  
1688 - " sec (reference ~ 0.2 sec)";  
1689 - #endif  
1690 - //  
1691 - // Modify 10000 x 10000 matrix with 10000 (exisiting) values  
1692 - t_start = clock();  
1693 - for (int i = 0; i < number; ++i)  
1694 - large(i,i) = (i+1)*1.0;  
1695 - t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);  
1696 - #if defined(DUMP_TIMING)  
1697 - cout << endl << " - Modification of 10000 values needed " << t_elapse <<  
1698 - " sec (reference ~ 0 sec)";  
1699 - #endif  
1700 - //  
1701 - // Insert columns into 10000 x 10000 matrix  
1702 - number = 10000;  
1703 - num_cols = 10000;  
1704 - GSparseMatrix diag(number,number);  
1705 - for (int j = 10; j < number-10; ++j)  
1706 - diag(j,j) = 3.14;  
1707 - large = diag;  
1708 - large.set_mem_block(block);  
1709 - GVector column(number);  
1710 - t_start = clock();  
1711 - for (int j = 0; j < num_cols; ++j) {  
1712 - i_min = (j < 2) ? 0 : j-2;  
1713 - i_max = (j > num_cols-2) ? num_cols : j+2;  
1714 - column = 0.0;  
1715 - for (int i = i_min; i < i_max; ++i)  
1716 - column[i] = (i+1)*0.01;  
1717 - large.add_col(column, j);  
1718 - }  
1719 - t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);  
1720 - #if defined(DUMP_TIMING)  
1721 - cout << endl << " - " << num_cols << " columns adding needed " << t_elapse <<  
1722 - " sec (reference ~ 4.6 sec)";  
1723 - #endif  
1724 - //  
1725 - // Insert columns into 10000 x 10000 matrix using matrix stack  
1726 - GSparseMatrix stack(number,number);  
1727 - stack = diag;  
1728 - stack.set_mem_block(block);  
1729 - for (int j = 200; j < 400; ++j)  
1730 - stack(j,j) = 3.14;  
1731 - t_start = clock();  
1732 - column = GVector(number);  
1733 - stack.stack_init(10000, 1000);  
1734 - for (int j = 0; j < num_cols; ++j) {  
1735 - i_min = (j < 2) ? 0 : j-2;  
1736 - i_max = (j > num_cols-2) ? num_cols : j+2;  
1737 - column = 0.0;  
1738 - for (int i = i_min; i < i_max; ++i)  
1739 - column[i] = (i+1)*0.01;  
1740 - stack.add_col(column, j);  
1741 - }  
1742 - stack.stack_destroy();  
1743 - t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);  
1744 - #if defined(DUMP_TIMING)  
1745 - cout << endl << " - " << num_cols << " columns stack-adding needed " << t_elapse <<  
1746 - " sec (reference ~ 1 sec)";  
1747 - #endif  
1748 - //  
1749 - // Check that both are identical  
1750 - if (large != stack) {  
1751 - cout << endl << "TEST ERROR: Stack-based insertion corrupted." << endl;  
1752 - for (int i = 0; i < number; ++i) {  
1753 - for (int j = 0; j < number; ++j) {  
1754 - if (large(i,j) != stack(i,j)) {  
1755 - std::cout << "(" << i << "," << j << "): ";  
1756 - std::cout << large(i,j) << " is not stack " << stack(i,j) << std::endl;  
1757 - } 1304 + // Set reference
  1305 + std::string reference;
  1306 + reference.append("=== GSparseMatrix ===\n");
  1307 + reference.append(" Number of rows ............: 4\n");
  1308 + reference.append(" Number of columns .........: 5\n");
  1309 + reference.append(" Number of nonzero elements : 9\n");
  1310 + reference.append(" Pending element ...........: (3,4)=9\n");
  1311 + reference.append(" Number of allocated cells .: 512\n");
  1312 + reference.append(" Memory block size .........: 512\n");
  1313 + reference.append(" Sparse matrix fill ........: 0.45\n");
  1314 + reference.append(" Pending element ...........: 9\n");
  1315 + reference.append(" Fill stack size ...........: 0 (none)\n");
  1316 + reference.append(" 1, 0, 0, 0, 7\n");
  1317 + reference.append(" 0, 3, 0, 0, 0\n");
  1318 + reference.append(" 2, 0, 4, 0, 8\n");
  1319 + reference.append(" 0, 0, 5, 6, 9");
  1320 +
  1321 + // Test the print method
  1322 + std::string output = m_test.print();
  1323 + test_assert((output == reference), "Test matrix printing",
  1324 + "Unexpected print() output:\n"+output+"\nExpected:\n"+reference);
  1325 +
  1326 + // Allocate big matrix
  1327 + GSparseMatrix big(30, 40);
  1328 + for (int i = 0; i < 30; ++i) {
  1329 + for (int j = 0; j < 40; ++j) {
  1330 + big(i,j) = i+j;
1758 } 1331 }
1759 - }  
1760 - throw;  
1761 - }  
1762 - //  
1763 - // Subsequent addition and subtractions  
1764 - number = 300;  
1765 - t_start = clock();  
1766 - GSparseMatrix sum(number,number);  
1767 - GSparseMatrix add(number,number);  
1768 - GSparseMatrix sub(number,number);  
1769 - sum.set_mem_block(block);  
1770 - add.set_mem_block(block);  
1771 - sub.set_mem_block(block);  
1772 - for (int i = 0; i < number; ++i) {  
1773 - add(i,number-i-1) = (i*i+1)*0.01;  
1774 - sum += add;  
1775 - sub(i,i) = (i+1)*1.0;  
1776 - sum -= sub;  
1777 - }  
1778 - t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);  
1779 - #if defined(DUMP_TIMING)  
1780 - cout << endl << " - Matrix adding and subtraction 300 times " << t_elapse <<  
1781 - " sec (reference ~ 1.4 sec)";  
1782 - #endif  
1783 - //  
1784 - // Subsequent multiplications  
1785 - GSparseMatrix fac(number,number);  
1786 - fac.set_mem_block(block);  
1787 - t_start = clock();  
1788 - for (int i = 0; i < 1; ++i) {  
1789 - for (int j = 0; j < number; ++j)  
1790 - fac(j,j) = (i+1)*1.0;  
1791 - sum *= fac;  
1792 - }  
1793 - t_elapse = double(clock() - t_start)/double(CLOCKS_PER_SEC);  
1794 - #if defined(DUMP_TIMING)  
1795 - cout << endl << " - Matrix multiplication 1 times " << t_elapse <<  
1796 - " sec (reference ~ 3.1 sec)";  
1797 - #endif  
1798 - }  
1799 - catch (exception &e) {  
1800 - cout << e.what() << endl;  
1801 - throw;  
1802 - }  
1803 - #if defined(DUMP_TIMING)  
1804 - cout << endl;  
1805 - #endif  
1806 - cout << "ok." << endl;  
1807 -}  
1808 -  
1809 -  
1810 -void test_stack(int block = 512)  
1811 -{  
1812 - // Set parameters  
1813 - int nrows = 10;  
1814 - int ncols = 10;  
1815 -  
1816 - // Initialise column vector  
1817 - GVector column(nrows);  
1818 - for (int i = 0; i < nrows; i+=2) {  
1819 - column[i] = 5.0;  
1820 } 1332 }
1821 1333
1822 - // Initialise reference matrix  
1823 - GSparseMatrix ref(nrows, ncols);  
1824 - for (int i = 0; i < ncols; i+=2) {  
1825 - ref.add_col(column, i);  
1826 - } 1334 + // Set reference for big matrix
  1335 + std::string ref_big;
  1336 + ref_big.append("=== GSparseMatrix ===\n");
  1337 + ref_big.append(" Number of rows ............: 30\n");
  1338 + ref_big.append(" Number of columns .........: 40\n");
  1339 + ref_big.append(" Number of nonzero elements : 1199\n");
  1340 + ref_big.append(" Pending element ...........: (29,39)=68\n");
  1341 + ref_big.append(" Number of allocated cells .: 1536\n");
  1342 + ref_big.append(" Memory block size .........: 512\n");
  1343 + ref_big.append(" Sparse matrix fill ........: 0.999167\n");
  1344 + ref_big.append(" Pending element ...........: 68\n");
  1345 + ref_big.append(" Fill stack size ...........: 0 (none)\n");
  1346 + ref_big.append(" 0, 1, 2, 3, 4, ... 35, 36, 37, 38, 39\n");
  1347 + ref_big.append(" 1, 2, 3, 4, 5, ... 36, 37, 38, 39, 40\n");
  1348 + ref_big.append(" 2, 3, 4, 5, 6, ... 37, 38, 39, 40, 41\n");
  1349 + ref_big.append(" 3, 4, 5, 6, 7, ... 38, 39, 40, 41, 42\n");
  1350 + ref_big.append(" 4, 5, 6, 7, 8, ... 39, 40, 41, 42, 43\n");
  1351 + ref_big.append(" ... ... ... ... ... ... ... ... ... ... ... \n");
  1352 + ref_big.append(" 25, 26, 27, 28, 29, ... 60, 61, 62, 63, 64\n");
  1353 + ref_big.append(" 26, 27, 28, 29, 30, ... 61, 62, 63, 64, 65\n");
  1354 + ref_big.append(" 27, 28, 29, 30, 31, ... 62, 63, 64, 65, 66\n");
  1355 + ref_big.append(" 28, 29, 30, 31, 32, ... 63, 64, 65, 66, 67\n");
  1356 + ref_big.append(" 29, 30, 31, 32, 33, ... 64, 65, 66, 67, 68");
  1357 +
  1358 + // Test the print method
  1359 + output = big.print();
  1360 + test_assert((output == ref_big), "Test big matrix printing",
  1361 + "Unexpected print() output:\n"+output+"\nExpected:\n"+ref_big);
1827 1362
1828 - // Initialise stack matrix  
1829 - GSparseMatrix stack(nrows, ncols);  
1830 - stack.stack_init(10000, 4);  
1831 - for (int i = 0; i < ncols; i+=2) {  
1832 - stack.add_col(column, i);  
1833 - }  
1834 - stack.stack_destroy();  
1835 -  
1836 - // Compare  
1837 - if (ref != stack) {  
1838 - std::cout << endl << "TEST ERROR: Stack-based insertion corrupted." << std::endl;  
1839 - for (int i = 0; i < nrows; ++i) {  
1840 - for (int j = 0; j < ncols; ++j) {  
1841 - if (ref(i,j) != stack(i,j)) {  
1842 - std::cout << "(" << i << "," << j << "): ";  
1843 - std::cout << "ref=" << ref(i,j) << " stack=" << stack(i,j) << std::endl;  
1844 - }  
1845 - }  
1846 - }  
1847 - throw;  
1848 - } 1363 + // Return
  1364 + return;
1849 } 1365 }
1850 1366
1851 1367
1852 /*************************************************************************** 1368 /***************************************************************************
1853 - * Main test function * 1369 + * @brief Main entry point for test executable
1854 ***************************************************************************/ 1370 ***************************************************************************/
1855 int main(void) 1371 int main(void)
1856 { 1372 {
1857 - // Dump header  
1858 - std::cout << std::endl;  
1859 - std::cout << "*******************************" << std::endl;  
1860 - std::cout << "* GSparseMatrix class testing *" << std::endl;  
1861 - std::cout << "*******************************" << std::endl;  
1862 -  
1863 - // Set test matrix and vector  
1864 - GSparseMatrix m_test = set_matrix();  
1865 - GVector v_test = set_vector();  
1866 -  
1867 - // Execute the tests  
1868 - test_output(m_test);  
1869 - test_allocate();  
1870 - test_assign_values(m_test);  
1871 - test_copy_constructor(m_test);  
1872 - test_assign_matrix(m_test);  
1873 - test_transpose(m_test);  
1874 - test_matrix_vector(m_test, v_test);  
1875 - test_matrix_matrix(m_test);  
1876 - test_arithmetics(m_test);  
1877 - test_functions(m_test);  
1878 - test_comparison(m_test);  
1879 - test_conversion();  
1880 - test_cholesky();  
1881 - test_heavy();  
1882 - //test_stack(); 1373 + // Allocate test suit container
  1374 + GTestSuites testsuites("GSparseMatrix class testing");
1883 1375
1884 - // Return  
1885 - return 0; 1376 + // Initially assume that we pass all tests
  1377 + bool success = true;
  1378 +
  1379 + // Create a test suite
  1380 + TestGSparseMatrix test;
  1381 +
  1382 + // Append test suite to the container
  1383 + testsuites.append(test);
  1384 +
  1385 + // Run the testsuites
  1386 + success = testsuites.run();
  1387 +
  1388 + // Save test report
  1389 + testsuites.save("reports/GSparseMatrix.xml");
  1390 +
  1391 + // Return success status
  1392 + return (success ? 0 : 1);
1886 } 1393 }
test/test_GSparseMatrix.hpp 0 → 100644
  View file @973bbf0
  1 +/***************************************************************************
  2 + * test_GSparseMatrix.hpp - Test sparse matrix class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file test_GSparseMatrix.hpp
  23 + * @brief Definition of unit tests for sparse matrices
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef TEST_GSPARSEMATRIX_HPP
  28 +#define TEST_GSPARSEMATRIX_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +
  33 +
  34 +/***********************************************************************//**
  35 + * @class TestGSparseMatrix
  36 + *
  37 + * @brief Test suite for sparse matrix class testing
  38 + *
  39 + * This test suite performs a unit test for the sparse 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 TestGSparseMatrix : public GTestSuite {
  44 +
  45 +public:
  46 + // Constructors and destructors
  47 + TestGSparseMatrix(void) : GTestSuite() {}
  48 + virtual ~TestGSparseMatrix(void) {}
  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);
  61 +
  62 +private:
  63 + // Private methods
  64 + GSparseMatrix set_matrix(void) const;
  65 + //GSparseMatrix set_matrix_zero(void) const;
  66 + GVector set_vector(void) const;
  67 + bool check_matrix(const GSparseMatrix& matrix,
  68 + const double& scale = 1.0,
  69 + const double& offset = 0.0) const;
  70 + bool check_matrix_trans(const GSparseMatrix& matrix,
  71 + const double& scale,
  72 + const double& offset) const;
  73 + bool check_matrix_lt(const GSparseMatrix& matrix,
  74 + const GSparseMatrix& ref) const;
  75 + bool check_matrix_ut(const GSparseMatrix& matrix,
  76 + const GSparseMatrix& ref) const;
  77 +
  78 + // Private members;
  79 + GSparseMatrix m_test; //!< Test matrix
  80 + GSparseMatrix m_bigger; //!< Bigger test matrix (for collisions)
  81 + GVector v_test; //!< Test vector
  82 +};
  83 +
  84 +#endif /* TEST_GSPARSEMATRIX_HPP */