GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
test_solvers.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: Unit tests for les solving
9 *
10 * COPYRIGHT: (C) 2000 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General Public
13 * License (>=v2). Read the file COPYING that comes with GRASS
14 * for details.
15 *
16 *****************************************************************************/
17 
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21 #include <grass/glocale.h>
22 #include <grass/N_pde.h>
23 #include "test_gpde_lib.h"
24 
25 /* prototypes */
26 static int test_solvers(void);
27 static N_les *create_normal_les(int rows);
28 static N_les *create_sparse_les(int rows);
29 
30 /* ************************************************************************* */
31 /* Performe the solver unit tests ****************************************** */
32 /* ************************************************************************* */
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 /* create a normal matrix with values ** Hilbert matrix ********** */
51 /* *************************************************************** */
52 N_les *create_normal_les(int rows)
53 {
54  N_les *les;
55  int i, j;
56  int size = rows;
57  double val;
58 
59  les = N_alloc_les(rows, N_NORMAL_LES);
60  for (i = 0; i < size; i++) {
61  val = 0.0;
62  for (j = 0; j < size; j++) {
63  les->A[i][j] = (double)(1.0 / (((double)i + 1.0) +
64  ((double)j + 1.0) - 1.0));
65  val += les->A[i][j];
66  }
67  les->b[i] = val;
68  }
69 
70  return les;
71 }
72 
73 /* *************************************************************** */
74 /* create a sparse matrix with values ** Hilbert matrix ********** */
75 /* *************************************************************** */
76 N_les *create_sparse_les(int rows)
77 {
78  N_les *les;
79  N_spvector *spvector;
80  int i, j;
81  double val;
82 
83  les = N_alloc_les(rows, N_SPARSE_LES);
84 
85  for (i = 0; i < rows; i++) {
86  spvector = N_alloc_spvector(rows);
87  val = 0;
88 
89  for (j = 0; j < rows; j++) {
90  spvector->values[j] =
91  (double)(1.0 / (((double)i + 1.0) + ((double)j + 1.0) - 1.0));
92  spvector->index[j] = j;
93  val += spvector->values[j];
94  }
95 
96  N_add_spvector_to_les(les, spvector, i);
97  les->b[i] = val;
98  }
99 
100 
101  return les;
102 }
103 
104 
105 /* *************************************************************** */
106 /* Test all implemented solvers for sparse and normal matrices *** */
107 /* *************************************************************** */
108 int test_solvers(void)
109 {
110  N_les *les;
111  N_les *sples;
112 
113  G_message("\t * testing jacobi solver\n");
114 
115  les = create_normal_les(TEST_N_NUM_ROWS);
116  sples = create_sparse_les(TEST_N_NUM_ROWS);
117 
118  N_solver_jacobi(les, 100, 1, 0.1e-4);
119  /*N_print_les(les); */
120  N_solver_jacobi(sples, 100, 1, 0.1e-4);
121  /*N_print_les(sples); */
122 
123  N_free_les(les);
124  N_free_les(sples);
125 
126 
127  G_message("\t * testing SOR solver\n");
128 
129  les = create_normal_les(TEST_N_NUM_ROWS);
130  sples = create_sparse_les(TEST_N_NUM_ROWS);
131 
132  N_solver_SOR(les, 100, 1, 0.1e-4);
133  /*N_print_les(les); */
134  N_solver_SOR(sples, 100, 1, 0.1e-4);
135  /*N_print_les(sples); */
136 
137  N_free_les(les);
138  N_free_les(sples);
139 
140  G_message("\t * testing cg solver\n");
141 
142  les = create_normal_les(TEST_N_NUM_ROWS);
143  sples = create_sparse_les(TEST_N_NUM_ROWS);
144 
145  N_solver_cg(les, 100, 0.1e-8);
146  /*N_print_les(les); */
147  N_solver_cg(sples, 100, 0.1e-8);
148  /*N_print_les(sples); */
149 
150  N_free_les(les);
151  N_free_les(sples);
152 
153  G_message("\t * testing pcg solver with N_DIAGONAL_PRECONDITION\n");
154 
155  les = create_normal_les(TEST_N_NUM_ROWS);
156  sples = create_sparse_les(TEST_N_NUM_ROWS);
157 
158  N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
159  N_print_les(les);
160  N_solver_pcg(sples, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
161  N_print_les(sples);
162 
163  N_free_les(les);
164  N_free_les(sples);
165 
166  G_message
167  ("\t * testing pcg solver with N_ROWSCALE_EUKLIDNORM_PRECONDITION\n");
168 
169  les = create_normal_les(TEST_N_NUM_ROWS);
170  sples = create_sparse_les(TEST_N_NUM_ROWS);
171 
173  N_print_les(les);
175  N_print_les(sples);
176 
177  N_free_les(les);
178  N_free_les(sples);
179 
180  G_message
181  ("\t * testing pcg solver with N_ROWSCALE_ABSSUMNORM_PRECONDITION\n");
182 
183  les = create_normal_les(TEST_N_NUM_ROWS);
184  sples = create_sparse_les(TEST_N_NUM_ROWS);
185 
187  N_print_les(les);
189  N_print_les(sples);
190 
191  N_free_les(les);
192  N_free_les(sples);
193 
194 
195  G_message("\t * testing bicgstab solver\n");
196 
197  les = create_normal_les(TEST_N_NUM_ROWS);
198  sples = create_sparse_les(TEST_N_NUM_ROWS);
199 
200  N_solver_bicgstab(les, 100, 0.1e-8);
201  /*N_print_les(les); */
202  N_solver_bicgstab(sples, 100, 0.1e-8);
203  /*N_print_les(sples); */
204 
205  N_free_les(les);
206  N_free_les(sples);
207 
208  G_message("\t * testing gauss elimination solver\n");
209 
210  les = create_normal_les(TEST_N_NUM_ROWS);
211 
212  /*GAUSS*/ N_solver_gauss(les);
213  N_print_les(les);
214 
215  N_free_les(les);
216 
217  G_message("\t * testing lu decomposition solver\n");
218 
219  les = create_normal_les(TEST_N_NUM_ROWS);
220 
221  /*LU*/ N_solver_lu(les);
222  N_print_les(les);
223 
224  N_free_les(les);
225 
226  les = create_normal_les(TEST_N_NUM_ROWS);
227 
228  /*cholesky */ N_solver_cholesky(les);
229  N_print_les(les);
230 
231  N_free_les(les);
232 
233 
234  return 0;
235 }