GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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 * */
48int 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 * */
118int 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 * */
190int 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 * */
247int 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: alloc.c:149
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
double b
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
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.
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_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_gs(double **A, double *x, double *b, int rows, int maxit, double sor, double error)
The iterative gauss seidel solver for quadratic matrices.
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
#define x