source: grass/trunk/lib/gmath/test/test_solvers.c

Last change on this file was 53679, checked in by neteler, 12 years ago

svn propset
  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-csrc
File size: 15.1 KB
Line 
1/*****************************************************************************
2 *
3 * MODULE: Grass PDE Numerical Library
4 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
5 * soerengebbert <at> gmx <dot> de
6 *
7 * PURPOSE: Unit tests for les solving
8 *
9 * COPYRIGHT: (C) 2000 by the GRASS Development Team
10 *
11 * This program is free software under the GNU General Public
12 * License (>=v2). Read the file COPYING that comes with GRASS
13 * for details.
14 *
15 *****************************************************************************/
16
17#include <stdio.h>
18#include <stdlib.h>
19#include <string.h>
20#include <grass/glocale.h>
21#include <grass/gmath.h>
22#include "test_gmath_lib.h"
23
24#define EPSILON_DIRECT 1.0E-10
25#define EPSILON_ITER 1.0E-4
26
27/* prototypes */
28static int test_solvers(void);
29
30/* ************************************************************************* */
31/* Performe the solver unit tests ****************************************** */
32/* ************************************************************************* */
33int unit_test_solvers(void)
34{
35 int sum = 0;
36
37 G_message(_("\n++ Running solver unit tests ++"));
38
39 sum += test_solvers();
40
41 if (sum > 0)
42 G_warning(_("\n-- Solver unit tests failure --"));
43 else
44 G_message(_("\n-- Solver unit tests finished successfully --"));
45
46 return sum;
47}
48
49/* *************************************************************** */
50/* Test all implemented solvers for sparse and normal matrix *** */
51/* *************************************************************** */
52int test_solvers(void)
53{
54 G_math_les *les;
55 G_math_les *sples;
56 int sum = 0;
57 double val = 0.0;
58
59 G_message("\t * testing jacobi solver with symmetric matrix\n");
60
61 les = create_normal_symmetric_les(TEST_NUM_ROWS);
62 sples = create_sparse_symmetric_les(TEST_NUM_ROWS);
63
64 G_math_solver_jacobi(les->A, les->x, les->b, les->rows, 250, 1, 0.1e-10);
65 G_math_d_asum_norm(les->x, &val, les->rows);
66 if ((val - (double)les->rows) > EPSILON_ITER)
67 {
68 G_warning("Error in G_math_solver_jacobi abs %2.20f != %i", val,
69 les->rows);
70 sum++;
71 }
72 G_math_solver_sparse_jacobi(sples->Asp, sples->x, sples->b, les->rows, 250,
73 1, 0.1e-10);
74 G_math_d_asum_norm(sples->x, &val, sples->rows);
75 if ((val - (double)sples->rows) > EPSILON_ITER)
76 {
77 G_warning("Error in G_math_solver_sparse_jacobi abs %2.20f != %i", val,
78 sples->rows);
79 sum++;
80 }
81
82 G_math_free_les(les);
83 G_math_free_les(sples);
84
85 G_message("\t * testing jacobi solver with unsymmetric matrix\n");
86
87 les = create_normal_unsymmetric_les(TEST_NUM_ROWS);
88 sples = create_sparse_unsymmetric_les(TEST_NUM_ROWS);
89
90 G_math_solver_jacobi(les->A, les->x, les->b, les->rows, 250, 1, 0.1e-10);
91 G_math_d_asum_norm(les->x, &val, les->rows);
92 if ((val - (double)les->rows) > EPSILON_ITER)
93 {
94 G_warning("Error in G_math_solver_jacobi abs %2.20f != %i", val,
95 les->rows);
96 sum++;
97 }
98
99 G_math_solver_sparse_jacobi(sples->Asp, sples->x, sples->b, les->rows, 250,
100 1, 0.1e-10);
101 G_math_d_asum_norm(sples->x, &val, sples->rows);
102 if ((val - (double)sples->rows) > EPSILON_ITER)
103 {
104 G_warning("Error in G_math_solver_sparse_jacobi abs %2.20f != %i", val,
105 sples->rows);
106 sum++;
107 }
108
109 G_math_free_les(les);
110 G_math_free_les(sples);
111
112 G_message("\t * testing gauss seidel solver with symmetric matrix\n");
113
114 les = create_normal_symmetric_les(TEST_NUM_ROWS);
115 sples = create_sparse_symmetric_les(TEST_NUM_ROWS);
116
117 G_math_solver_gs(les->A, les->x, les->b, les->rows, 150, 1, 0.1e-9);
118 G_math_d_asum_norm(les->x, &val, les->rows);
119 if ((val - (double)les->rows) > EPSILON_ITER)
120 {
121 G_warning("Error in G_math_solver_gs abs %2.20f != %i", val, les->rows);
122 sum++;
123 }
124
125 G_math_solver_sparse_gs(sples->Asp, sples->x, sples->b, les->rows, 150, 1,
126 0.1e-9);
127 G_math_d_asum_norm(sples->x, &val, sples->rows);
128 if ((val - (double)sples->rows) > EPSILON_ITER)
129 {
130 G_warning("Error in G_math_solver_sparse_gs abs %2.20f != %i", val,
131 sples->rows);
132 sum++;
133 }
134
135 G_math_free_les(les);
136 G_math_free_les(sples);
137
138 G_message("\t * testing gauss seidel solver with unsymmetric matrix\n");
139
140 les = create_normal_unsymmetric_les(TEST_NUM_ROWS);
141 sples = create_sparse_unsymmetric_les(TEST_NUM_ROWS);
142
143 G_math_solver_gs(les->A, les->x, les->b, les->rows, 150, 1, 0.1e-9);
144 G_math_d_asum_norm(les->x, &val, les->rows);
145 if ((val - (double)les->rows) > EPSILON_ITER)
146 {
147 G_warning("Error in G_math_solver_gs abs %2.20f != %i", val, les->rows);
148 sum++;
149 }
150
151 G_math_solver_sparse_gs(sples->Asp, sples->x, sples->b, les->rows, 150, 1,
152 0.1e-9);
153 G_math_d_asum_norm(sples->x, &val, sples->rows);
154 if ((val - (double)sples->rows) > EPSILON_ITER)
155 {
156 G_warning("Error in G_math_solver_sparse_gs abs %2.20f != %i", val,
157 sples->rows);
158 sum++;
159 }
160
161 G_math_free_les(les);
162 G_math_free_les(sples);
163
164 G_message("\t * testing pcg solver with symmetric bad conditioned matrix and preconditioner 3\n");
165
166 les = create_normal_symmetric_pivot_les(TEST_NUM_ROWS);
167
168 G_math_solver_pcg(les->A, les->x, les->b, les->rows, 250, 0.1e-9, 3);
169 G_math_d_asum_norm(les->x, &val, les->rows);
170 if ((val - (double)les->rows) > EPSILON_ITER)
171 {
172 G_warning("Error in G_math_solver_pcg abs %2.20f != %i", val, les->rows);
173 sum++;
174 }
175 G_math_print_les(les);
176 G_math_free_les(les);
177
178 G_message("\t * testing pcg solver with symmetric matrix and preconditioner 1\n");
179
180 les = create_normal_symmetric_les(TEST_NUM_ROWS);
181 sples = create_sparse_symmetric_les(TEST_NUM_ROWS);
182
183 G_math_solver_pcg(les->A, les->x, les->b, les->rows, 250, 0.1e-9, 1);
184 G_math_d_asum_norm(les->x, &val, les->rows);
185 if ((val - (double)les->rows) > EPSILON_ITER)
186 {
187 G_warning("Error in G_math_solver_pcg abs %2.20f != %i", val, les->rows);
188 sum++;
189 }
190 G_math_print_les(les);
191
192 G_math_solver_sparse_pcg(sples->Asp, sples->x, sples->b, les->rows, 250,
193 0.1e-9, 1);
194 G_math_d_asum_norm(sples->x, &val, sples->rows);
195 if ((val - (double)sples->rows) > EPSILON_ITER)
196 {
197 G_warning("Error in G_math_solver_sparse_pcg abs %2.20f != %i", val,
198 sples->rows);
199 sum++;
200 }
201 G_math_print_les(sples);
202
203 G_math_free_les(les);
204 G_math_free_les(sples);
205
206 G_message("\t * testing pcg solver with symmetric matrix and preconditioner 2\n");
207
208 les = create_normal_symmetric_les(TEST_NUM_ROWS);
209 sples = create_sparse_symmetric_les(TEST_NUM_ROWS);
210
211 G_math_solver_pcg(les->A, les->x, les->b, les->rows, 250, 0.1e-9, 2);
212 G_math_d_asum_norm(les->x, &val, les->rows);
213 if ((val - (double)les->rows) > EPSILON_ITER)
214 {
215 G_warning("Error in G_math_solver_pcg abs %2.20f != %i", val, les->rows);
216 sum++;
217 }
218 G_math_print_les(les);
219
220 G_math_solver_sparse_pcg(sples->Asp, sples->x, sples->b, les->rows, 250,
221 0.1e-9, 2);
222 G_math_d_asum_norm(sples->x, &val, sples->rows);
223 if ((val - (double)sples->rows) > EPSILON_ITER)
224 {
225 G_warning("Error in G_math_solver_sparse_pcg abs %2.20f != %i", val,
226 sples->rows);
227 sum++;
228 }
229 G_math_print_les(sples);
230
231 G_math_free_les(les);
232 G_math_free_les(sples);
233
234 G_message("\t * testing pcg solver with symmetric matrix and preconditioner 3\n");
235
236 les = create_normal_symmetric_les(TEST_NUM_ROWS);
237 sples = create_sparse_symmetric_les(TEST_NUM_ROWS);
238
239 G_math_solver_pcg(les->A, les->x, les->b, les->rows, 250, 0.1e-9, 3);
240 G_math_d_asum_norm(les->x, &val, les->rows);
241 if ((val - (double)les->rows) > EPSILON_ITER)
242 {
243 G_warning("Error in G_math_solver_pcg abs %2.20f != %i", val, les->rows);
244 sum++;
245 }
246 G_math_print_les(les);
247
248 G_math_solver_sparse_pcg(sples->Asp, sples->x, sples->b, les->rows, 250,
249 0.1e-9, 3);
250 G_math_d_asum_norm(sples->x, &val, sples->rows);
251 if ((val - (double)sples->rows) > EPSILON_ITER)
252 {
253 G_warning("Error in G_math_solver_sparse_pcg abs %2.20f != %i", val,
254 sples->rows);
255 sum++;
256 }
257 G_math_print_les(sples);
258
259 G_math_free_les(les);
260 G_math_free_les(sples);
261
262 G_message("\t * testing cg solver with symmetric matrix\n");
263
264 les = create_normal_symmetric_les(TEST_NUM_ROWS);
265 sples = create_sparse_symmetric_les(TEST_NUM_ROWS);
266
267 G_math_solver_cg(les->A, les->x, les->b, les->rows, 250, 0.1e-9);
268 G_math_d_asum_norm(les->x, &val, les->rows);
269 if ((val - (double)les->rows) > EPSILON_ITER)
270 {
271 G_warning("Error in G_math_solver_cg abs %2.20f != %i", val, les->rows);
272 sum++;
273 }
274 G_math_print_les(les);
275
276 G_math_solver_sparse_cg(sples->Asp, sples->x, sples->b, les->rows, 250,
277 0.1e-9);
278 G_math_d_asum_norm(sples->x, &val, sples->rows);
279 if ((val - (double)sples->rows) > EPSILON_ITER)
280 {
281 G_warning("Error in G_math_solver_sparse_cg abs %2.20f != %i", val,
282 sples->rows);
283 sum++;
284 }
285 G_math_print_les(sples);
286
287 G_math_free_les(les);
288 G_math_free_les(sples);
289
290
291 G_message("\t * testing cg solver with symmetric bad conditioned matrix\n");
292
293 les = create_normal_symmetric_pivot_les(TEST_NUM_ROWS);
294
295 G_math_solver_cg(les->A, les->x, les->b, les->rows, 250, 0.1e-9);
296 G_math_d_asum_norm(les->x, &val, les->rows);
297 if ((val - (double)les->rows) > EPSILON_ITER)
298 {
299 G_warning("Error in G_math_solver_cg abs %2.20f != %i", val, les->rows);
300 sum++;
301 }
302 G_math_print_les(les);
303 G_math_free_les(les);
304
305 G_message("\t * testing bicgstab solver with symmetric matrix\n");
306
307 les = create_normal_symmetric_les(TEST_NUM_ROWS);
308 sples = create_sparse_symmetric_les(TEST_NUM_ROWS);
309
310 G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 250, 0.1e-9);
311 G_math_d_asum_norm(les->x, &val, les->rows);
312 if ((val - (double)les->rows) > EPSILON_ITER)
313 {
314 G_warning("Error in G_math_solver_bicgstab abs %2.20f != %i", val,
315 les->rows);
316 sum++;
317 }
318 G_math_print_les(les);
319
320 G_math_solver_sparse_bicgstab(sples->Asp, sples->x, sples->b, les->rows,
321 250, 0.1e-9);
322 G_math_d_asum_norm(sples->x, &val, sples->rows);
323 if ((val - (double)sples->rows) > EPSILON_ITER)
324 {
325 G_warning("Error in G_math_solver_sparse_bicgstab abs %2.20f != %i",
326 val, sples->rows);
327 sum++;
328 }
329 G_math_print_les(sples);
330
331 G_math_free_les(les);
332 G_math_free_les(sples);
333
334 G_message("\t * testing bicgstab solver with unsymmetric matrix\n");
335
336 les = create_normal_unsymmetric_les(TEST_NUM_ROWS);
337 sples = create_sparse_unsymmetric_les(TEST_NUM_ROWS);
338
339 G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 250, 0.1e-9);
340 G_math_d_asum_norm(les->x, &val, les->rows);
341 if ((val - (double)les->rows) > EPSILON_ITER)
342 {
343 G_warning("Error in G_math_solver_bicgstab abs %2.20f != %i", val,
344 les->rows);
345 sum++;
346 }
347 G_math_print_les(les);
348
349 G_math_solver_sparse_bicgstab(sples->Asp, sples->x, sples->b, les->rows,
350 250, 0.1e-9);
351 G_math_d_asum_norm(sples->x, &val, sples->rows);
352 if ((val - (double)les->rows) > EPSILON_ITER)
353 {
354 G_warning("Error in G_math_solver_sparse_bicgstab abs %2.20f != %i",
355 val, sples->rows);
356 sum++;
357 }
358
359 G_math_print_les(sples);
360 G_math_free_les(les);
361 G_math_free_les(sples);
362
363 G_message("\t * testing gauss elimination solver with symmetric matrix\n");
364 les = create_normal_symmetric_les(TEST_NUM_ROWS);
365 G_math_solver_gauss(les->A, les->x, les->b, les->rows);
366 G_math_d_asum_norm(les->x, &val, les->rows);
367 if ((val - (double)les->rows) > EPSILON_DIRECT)
368 {
369 G_warning("Error in G_math_solver_gauss abs %2.20f != %i", val,
370 les->rows);
371 sum++;
372 }
373 G_math_print_les(les);
374 G_math_free_les(les);
375
376 G_message("\t * testing lu decomposition solver with symmetric matrix\n");
377 les = create_normal_symmetric_les(TEST_NUM_ROWS);
378 G_math_solver_lu(les->A, les->x, les->b, les->rows);
379 G_math_d_asum_norm(les->x, &val, les->rows);
380 if ((val - (double)les->rows) > EPSILON_DIRECT)
381 {
382 G_warning("Error in G_math_solver_gauss abs %2.20f != %i", val,
383 les->rows);
384 sum++;
385 }
386 G_math_print_les(les);
387 G_math_free_les(les);
388
389 G_message("\t * testing gauss elimination solver with unsymmetric matrix\n");
390 les = create_normal_unsymmetric_les(TEST_NUM_ROWS);
391 G_math_solver_gauss(les->A, les->x, les->b, les->rows);
392 G_math_d_asum_norm(les->x, &val, les->rows);
393 if ((val - (double)les->rows) > EPSILON_DIRECT)
394 {
395 G_warning("Error in G_math_solver_gauss abs %2.20f != %i", val,
396 les->rows);
397 sum++;
398 }
399 G_math_print_les(les);
400 G_math_free_les(les);
401
402 G_message("\t * testing lu decomposition solver with unsymmetric matrix\n");
403 les = create_normal_unsymmetric_les(TEST_NUM_ROWS);
404 G_math_solver_lu(les->A, les->x, les->b, les->rows);
405
406 G_math_d_asum_norm(les->x, &val, les->rows);
407 if ((val - (double)les->rows) > EPSILON_DIRECT)
408 {
409 G_warning("Error in G_math_solver_gauss abs %2.20f != %i", val,
410 les->rows);
411 sum++;
412 }
413 G_math_print_les(les);
414 G_math_free_les(les);
415
416 G_message("\t * testing gauss elimination solver with symmetric bad conditioned matrix\n");
417 les = create_normal_symmetric_pivot_les(TEST_NUM_ROWS);
418 G_math_solver_gauss(les->A, les->x, les->b, les->rows);
419 G_math_d_asum_norm(les->x, &val, les->rows);
420 if ((val - (double)les->rows) > EPSILON_DIRECT)
421 {
422 G_warning("Error in G_math_solver_gauss abs %2.20f != %i", val,
423 les->rows);
424 sum++;
425 }
426 G_math_print_les(les);
427 G_math_free_les(les);
428
429 G_message("\t * testing lu decomposition solver with symmetric bad conditioned matrix\n");
430 les = create_normal_symmetric_pivot_les(TEST_NUM_ROWS);
431 G_math_solver_lu(les->A, les->x, les->b, les->rows);
432 G_math_d_asum_norm(les->x, &val, les->rows);
433 if ((val - (double)les->rows) > EPSILON_DIRECT)
434 {
435 G_warning("Error in G_math_solver_lu abs %2.20f != %i", val,
436 les->rows);
437 sum++;
438 }
439 G_math_print_les(les);
440 G_math_free_les(les);
441
442 G_message("\t * testing cholesky decomposition solver with symmetric matrix\n");
443 les = create_normal_symmetric_les(TEST_NUM_ROWS);
444 /*cholesky*/G_math_solver_cholesky(les->A, les->x, les->b, les->rows,
445 les->rows);
446 G_math_d_asum_norm(les->x, &val, les->rows);
447 if ((val - (double)les->rows) > EPSILON_DIRECT)
448 {
449 G_warning("Error in G_math_solver_solver_cholesky abs %2.20f != %i", val,
450 les->rows);
451 sum++;
452 }
453 G_math_print_les(les);
454 G_math_free_les(les);
455
456 G_message("\t * testing cholesky band decomposition solver with symmetric band matrix 1\n");
457 les = create_normal_symmetric_les(TEST_NUM_ROWS);
458 G_math_print_les(les);
459 /* Create a band matrix*/
460 G_message("\t * Creating symmetric band matrix\n");
461 les->A = G_math_matrix_to_sband_matrix(les->A, les->rows, les->rows);
462 G_math_print_les(les);
463
464 /*cholesky*/G_math_solver_cholesky_sband(les->A, les->x, les->b, les->rows,les->rows);
465 G_math_d_asum_norm(les->x, &val, les->rows);
466 if ((val - (double)les->rows) > EPSILON_DIRECT)
467 {
468 G_warning("Error in G_math_solver_solver_cholesky_band abs %2.20f != %i", val,
469 les->rows);
470 sum++;
471 }
472 G_math_print_les(les);
473 G_math_free_les(les);
474
475 G_message("\t * testing cholesky band decomposition solver with symmetric band matrix 2\n");
476 les = create_symmetric_band_les(TEST_NUM_ROWS);
477 G_math_print_les(les);
478
479 /*cholesky*/G_math_solver_cholesky_sband(les->A, les->x, les->b, les->rows,les->rows);
480 G_math_d_asum_norm(les->x, &val, les->rows);
481 if ((val - (double)les->rows) > EPSILON_DIRECT)
482 {
483 G_warning("Error in G_math_solver_solver_cholesky_band abs %2.20f != %i", val,
484 les->rows);
485 sum++;
486 }
487 G_math_print_les(les);
488 G_math_free_les(les);
489
490 G_message("\t * testing cg solver with symmetric band matrix\n");
491
492 les = create_symmetric_band_les(TEST_NUM_ROWS);
493
494 G_math_solver_cg_sband(les->A, les->x, les->b, les->rows, les->rows, 250, 0.1e-9);
495 G_math_d_asum_norm(les->x, &val, les->rows);
496 if ((val - (double)les->rows) > EPSILON_ITER)
497 {
498 G_warning("Error in G_math_solver_cg_sband abs %2.20f != %i", val, les->rows);
499 sum++;
500 }
501 G_math_print_les(les);
502 G_math_free_les(les);
503
504 return sum;
505}
506
Note: See TracBrowser for help on using the repository browser.