GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
N_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) 2000 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/N_pde.h"
24 #include "solvers_local_proto.h"
25 
26 static int sparse_jacobi_gauss(N_les * L, int maxit, double sor, double error,
27  const char *type);
28 static int jacobi(double **M, double *b, double *x, int rows, int maxit,
29  double sor, double error);
30 static int gauss_seidel(double **M, double *b, double *x, int rows, int maxit,
31  double sor, double error);
32 
33 /* ******************************************************* *
34  * ******** overrelaxed jacobian ************************* *
35  * ******************************************************* */
53 int N_solver_jacobi(N_les * L, int maxit, double sor, double error)
54 {
55 
56  if (L->quad != 1) {
57  G_warning(_("The linear equation system is not quadratic"));
58  return -1;
59  }
60 
61  if (L->type == N_NORMAL_LES)
62  return jacobi(L->A, L->b, L->x, L->rows, maxit, sor, error);
63  else
64  return sparse_jacobi_gauss(L, maxit, sor, error,
66 }
67 
68 
69 /* ******************************************************* *
70  * ********* overrelaxed gauss seidel ******************** *
71  * ******************************************************* */
90 int N_solver_SOR(N_les * L, int maxit, double sor, double error)
91 {
92 
93  if (L->quad != 1) {
94  G_warning(_("The linear equation system is not quadratic"));
95  return -1;
96  }
97 
98  if (L->type == N_NORMAL_LES)
99  return gauss_seidel(L->A, L->b, L->x, L->rows, maxit, sor, error);
100  else
101  return sparse_jacobi_gauss(L, maxit, sor, error,
103 }
104 
105 /* ******************************************************* *
106  * ****** sparse jacobi and SOR algorithm **************** *
107  * ******************************************************* */
108 int
109 sparse_jacobi_gauss(N_les * L, int maxit, double sor, double error,
110  const char *type)
111 {
112  int i, j, k, rows, finished = 0;
113  double *Enew, *x, *b;
114  double E, err = 0;
115 
116  x = L->x;
117  b = L->b;
118  rows = L->rows;
119 
120  Enew = vectmem(rows);
121 
122  for (k = 0; k < maxit; k++) {
123  err = 0;
124  {
125  if (k == 0) {
126  for (j = 0; j < rows; j++) {
127  Enew[j] = x[j];
128  }
129  }
130  for (i = 0; i < rows; i++) {
131  E = 0;
132  if (strcmp(type, N_SOLVER_ITERATIVE_JACOBI) == 0) {
133  for (j = 0; j < L->Asp[i]->cols; j++) {
134  E += L->Asp[i]->values[j] * x[L->Asp[i]->index[j]];
135  }
136  }
137  else {
138  for (j = 0; j < L->Asp[i]->cols; j++) {
139  E += L->Asp[i]->values[j] * Enew[L->Asp[i]->index[j]];
140  }
141  }
142  Enew[i] = x[i] - sor * (E - b[i]) / L->Asp[i]->values[0];
143  }
144  for (j = 0; j < rows; j++) {
145  err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
146 
147  x[j] = Enew[j];
148  }
149  }
150 
151  if (strcmp(type, N_SOLVER_ITERATIVE_JACOBI) == 0)
152  G_message(_("sparse Jacobi -- iteration %5i error %g\n"), k, err);
153  else if (strcmp(type, N_SOLVER_ITERATIVE_SOR) == 0)
154  G_message(_("sparse SOR -- iteration %5i error %g\n"), k, err);
155 
156  if (err < error) {
157  finished = 1;
158  break;
159  }
160  }
161 
162  G_free(Enew);
163 
164  return finished;
165 }
166 
167 /* ******************************************************* *
168  * ******* direct jacobi ********************************* *
169  * ******************************************************* */
170 int jacobi(double **M, double *b, double *x, int rows, int maxit, double sor,
171  double error)
172 {
173  int i, j, k;
174  double *Enew;
175  double E, err = 0;
176 
177  Enew = vectmem(rows);
178 
179  for (j = 0; j < rows; j++) {
180  Enew[j] = x[j];
181  }
182 
183  for (k = 0; k < maxit; k++) {
184  for (i = 0; i < rows; i++) {
185  E = 0;
186  for (j = 0; j < rows; j++) {
187  E += M[i][j] * x[j];
188  }
189  Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
190  }
191  err = 0;
192  for (j = 0; j < rows; j++) {
193  err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
194 
195  x[j] = Enew[j];
196  }
197  G_message(_("Jacobi -- iteration %5i error %g\n"), k, err);
198  if (err < error)
199  break;
200  }
201 
202  return 1;
203 }
204 
205 /* ******************************************************* *
206  * ******* direct gauss seidel *************************** *
207  * ******************************************************* */
208 int gauss_seidel(double **M, double *b, double *x, int rows, int maxit,
209  double sor, double error)
210 {
211  int i, j, k;
212  double *Enew;
213  double E, err = 0;
214 
215  Enew = vectmem(rows);
216 
217  for (j = 0; j < rows; j++) {
218  Enew[j] = x[j];
219  }
220 
221  for (k = 0; k < maxit; k++) {
222  for (i = 0; i < rows; i++) {
223  E = 0;
224  for (j = 0; j < rows; j++) {
225  E += M[i][j] * Enew[j];
226  }
227  Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
228  }
229  err = 0;
230  for (j = 0; j < rows; j++) {
231  err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
232 
233  x[j] = Enew[j];
234  }
235  G_message(_("SOR -- iteration %5i error %g\n"), k, err);
236  if (err < error)
237  break;
238  }
239 
240  return 1;
241 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
float b
Definition: named_colr.c:8
int quad
Definition: N_pde.h:104
#define N_SOLVER_ITERATIVE_JACOBI
Definition: N_pde.h:29
int cols
Definition: N_pde.h:81
const char * err
Definition: g3dcolor.c:50
double * x
Definition: N_pde.h:98
int jacobi(double a[MX][MX], long n, double d[MX], double v[MX][MX])
Definition: jacobi.c:11
#define N_SOLVER_ITERATIVE_SOR
Definition: N_pde.h:30
double ** A
Definition: N_pde.h:100
int * index
Definition: N_pde.h:83
void G_message(const char *msg,...)
Print a message to stderr.
Definition: lib/gis/error.c:74
double * vectmem(int rows)
Allocate vector memory.
Definition: N_solvers.c:423
int N_solver_jacobi(N_les *L, int maxit, double sor, double error)
The iterative jacobian solver for regular matrices.
double * values
Definition: N_pde.h:82
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
double * b
Definition: N_pde.h:99
#define N_NORMAL_LES
Definition: N_pde.h:41
int rows
Definition: N_pde.h:102
int N_solver_SOR(N_les *L, int maxit, double sor, double error)
The iterative overrelaxed gauss seidel solver for regular matrices.
The linear equation system (les) structure.
Definition: N_pde.h:96
N_spvector ** Asp
Definition: N_pde.h:101
int type
Definition: N_pde.h:105