GRASS Programmer's Manual  6.5.svn(2014)-r66266
solvers_classic_iter.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: linear equation system solvers
9 * part of the gpde library
10 *
11 * COPYRIGHT: (C) 2007 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18
19 #include <math.h>
20 #include <unistd.h>
21 #include <stdio.h>
22 #include <string.h>
23 #include <grass/gis.h>
24 #include <grass/glocale.h>
25 #include <grass/gmath.h>
26
27
48 int G_math_solver_sparse_jacobi(G_math_spvector ** Asp, double *x, double *b,
49  int rows, int maxit, double sor, double error)
50 {
51  int i, j, k, center, finished = 0;
52
53  double *Enew;
54
55  double E, err = 0;
56
57  Enew = G_alloc_vector(rows);
58
59  for (k = 0; k < maxit; k++) {
60  err = 0;
61  {
62  if (k == 0) {
63  for (j = 0; j < rows; j++) {
64  Enew[j] = x[j];
65  }
66  }
67  for (i = 0; i < rows; i++) {
68  E = 0;
69  center = 0;
70  for (j = 0; j < Asp[i]->cols; j++) {
71  E += Asp[i]->values[j] * x[Asp[i]->index[j]];
72  if (Asp[i]->index[j] == i)
73  center = j;
74  }
75  Enew[i] = x[i] - sor * (E - b[i]) / Asp[i]->values[center];
76  }
77  for (j = 0; j < rows; j++) {
78  err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
79
80  x[j] = Enew[j];
81  }
82  }
83
84  G_message(_("sparse Jacobi -- iteration %5i error %g\n"), k, err);
85
86  if (err < error) {
87  finished = 1;
88  break;
89  }
90  }
91
92  G_free(Enew);
93
94  return finished;
95 }
96
97
118 int G_math_solver_sparse_gs(G_math_spvector ** Asp, double *x, double *b,
119  int rows, int maxit, double sor, double error)
120 {
121  int i, j, k, finished = 0;
122
123  double *Enew;
124
125  double E, err = 0;
126
127  int center;
128
129  Enew = G_alloc_vector(rows);
130
131  for (k = 0; k < maxit; k++) {
132  err = 0;
133  {
134  if (k == 0) {
135  for (j = 0; j < rows; j++) {
136  Enew[j] = x[j];
137  }
138  }
139  for (i = 0; i < rows; i++) {
140  E = 0;
141  center = 0;
142  for (j = 0; j < Asp[i]->cols; j++) {
143  E += Asp[i]->values[j] * Enew[Asp[i]->index[j]];
144  if (Asp[i]->index[j] == i)
145  center = j;
146  }
147  Enew[i] = x[i] - sor * (E - b[i]) / Asp[i]->values[center];
148  }
149  for (j = 0; j < rows; j++) {
150  err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
151
152  x[j] = Enew[j];
153  }
154  }
155
156  G_message(_("sparse SOR -- iteration %5i error %g\n"), k, err);
157
158  if (err < error) {
159  finished = 1;
160  break;
161  }
162  }
163
164  G_free(Enew);
165
166  return finished;
167 }
168
169
190 int G_math_solver_jacobi(double **A, double *x, double *b, int rows,
191  int maxit, double sor, double error)
192 {
193  int i, j, k;
194
195  double *Enew;
196
197  double E, err = 0;
198
199  Enew = G_alloc_vector(rows);
200
201  for (j = 0; j < rows; j++) {
202  Enew[j] = x[j];
203  }
204
205  for (k = 0; k < maxit; k++) {
206  for (i = 0; i < rows; i++) {
207  E = 0;
208  for (j = 0; j < rows; j++) {
209  E += A[i][j] * x[j];
210  }
211  Enew[i] = x[i] - sor * (E - b[i]) / A[i][i];
212  }
213  err = 0;
214  for (j = 0; j < rows; j++) {
215  err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
216  x[j] = Enew[j];
217  }
218  G_message(_("Jacobi -- iteration %5i error %g\n"), k, err);
219  if (err < error)
220  break;
221  }
222
223  return 1;
224 }
225
226
247 int G_math_solver_gs(double **A, double *x, double *b, int rows, int maxit,
248  double sor, double error)
249 {
250  int i, j, k;
251
252  double *Enew;
253
254  double E, err = 0;
255
256  Enew = G_alloc_vector(rows);
257
258  for (j = 0; j < rows; j++) {
259  Enew[j] = x[j];
260  }
261
262  for (k = 0; k < maxit; k++) {
263  for (i = 0; i < rows; i++) {
264  E = 0;
265  for (j = 0; j < rows; j++) {
266  E += A[i][j] * Enew[j];
267  }
268  Enew[i] = x[i] - sor * (E - b[i]) / A[i][i];
269  }
270  err = 0;
271  for (j = 0; j < rows; j++) {
272  err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
273  x[j] = Enew[j];
274  }
275  G_message(_("SOR -- iteration %5i error %g\n"), k, err);
276  if (err < error)
277  break;
278  }
279
280  return 1;
281 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
float b
Definition: named_colr.c:8
const char * err
Definition: g3dcolor.c:50
int G_math_solver_sparse_jacobi(G_math_spvector **Asp, double *x, double *b, int rows, int maxit, double sor, double error)
The iterative jacobi solver for sparse matrices.
void G_message(const char *msg,...)
Print a message to stderr.
Definition: lib/gis/error.c:74
int G_math_solver_jacobi(double **A, double *x, double *b, int rows, int maxit, double sor, double error)
The iterative jacobi solver for quadratic matrices.
int G_math_solver_sparse_gs(G_math_spvector **Asp, double *x, double *b, int rows, int maxit, double sor, double error)
The iterative gauss seidel solver for sparse matrices.
int G_math_solver_gs(double **A, double *x, double *b, int rows, int maxit, double sor, double error)
The iterative gauss seidel solver for quadratic matrices.
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41