GRASS GIS 7 Programmer's Manual  7.7.svn(2018)-r73786
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
solvers_classic_iter.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass numerical math interface
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> googlemail <dot> com
7 *
8 * PURPOSE: linear equation system solvers
9 * part of the gmath library
10 *
11 * COPYRIGHT: (C) 2010 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 
28 /*!
29  * \brief The iterative jacobi solver for sparse matrices
30  *
31  * The Jacobi solver solves the linear equation system Ax = b
32  * The result is written to the vector x.
33  *
34  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
35  * solver will abort the calculation and writes the current result into the vector x.
36  * The parameter <i>err</i> defines the error break criteria for the solver.
37  *
38  * \param Asp G_math_spvector ** -- the sparse matrix
39  * \param x double * -- the vector of unknowns
40  * \param b double * -- the right side vector
41  * \param rows int -- number of rows
42  * \param maxit int -- the maximum number of iterations
43  * \param sor double -- defines the successive overrelaxion parameter [0:1]
44  * \param error double -- defines the error break criteria
45  * \return int -- 1=success, -1=could not solve the les
46  *
47  * */
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 
98 /*!
99  * \brief The iterative gauss seidel solver for sparse matrices
100  *
101  * The Jacobi solver solves the linear equation system Ax = b
102  * The result is written to the vector x.
103  *
104  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
105  * solver will abort the calculation and writes the current result into the vector x.
106  * The parameter <i>err</i> defines the error break criteria for the solver.
107  *
108  * \param Asp G_math_spvector ** -- the sparse matrix
109  * \param x double * -- the vector of unknowns
110  * \param b double * -- the right side vector
111  * \param rows int -- number of rows
112  * \param maxit int -- the maximum number of iterations
113  * \param sor double -- defines the successive overrelaxion parameter [0:2]
114  * \param error double -- defines the error break criteria
115  * \return int -- 1=success, -1=could not solve the les
116  *
117  * */
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 
170 /*!
171  * \brief The iterative jacobi solver for quadratic matrices
172  *
173  * The Jacobi solver solves the linear equation system Ax = b
174  * The result is written to the vector x.
175  *
176  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
177  * solver will abort the calculation and writes the current result into the vector x.
178  * The parameter <i>err</i> defines the error break criteria for the solver.
179  *
180  * \param A double ** -- the dense matrix
181  * \param x double * -- the vector of unknowns
182  * \param b double * -- the right side vector
183  * \param rows int -- number of rows
184  * \param maxit int -- the maximum number of iterations
185  * \param sor double -- defines the successive overrelaxion parameter [0:1]
186  * \param error double -- defines the error break criteria
187  * \return int -- 1=success, -1=could not solve the les
188  *
189  * */
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 
227 /*!
228  * \brief The iterative gauss seidel solver for quadratic matrices
229  *
230  * The Jacobi solver solves the linear equation system Ax = b
231  * The result is written to the vector x.
232  *
233  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
234  * solver will abort the calculation and writes the current result into the vector x.
235  * The parameter <i>err</i> defines the error break criteria for the solver.
236  *
237  * \param A double ** -- the dense matrix
238  * \param x double * -- the vector of unknowns
239  * \param b double * -- the right side vector
240  * \param rows int -- number of rows
241  * \param maxit int -- the maximum number of iterations
242  * \param sor double -- defines the successive overrelaxion parameter [0:2]
243  * \param error double -- defines the error break criteria
244  * \return int -- 1=success, -1=could not solve the les
245  *
246  * */
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:149
The row vector of the sparse matrix.
Definition: gmath.h:54
#define x
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
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.
unsigned int * index
Definition: gmath.h:58
double b
Definition: r_raster.c:39
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.
double * values
Definition: gmath.h:56
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.
unsigned int cols
Definition: gmath.h:57
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.
#define _(str)
Definition: glocale.h:13
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41