GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
solvers_classic_iter.c
Go to the documentation of this file.
1/*****************************************************************************
2 *
3 * MODULE: Grass numerical math interface
4 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
5 * soerengebbert <at> googlemail <dot> com
6 *
7 * PURPOSE: linear equation system solvers
8 * part of the gmath library
9 *
10 * COPYRIGHT: (C) 2010 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 <assert.h>
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 * \brief The iterative jacobi solver for sparse matrices
29 *
30 * The Jacobi solver solves the linear equation system Ax = b
31 * The result is written to the vector x.
32 *
33 * The parameter <i>maxit</i> specifies the maximum number of iterations. If the
34 * maximum is reached, the solver will abort the calculation and writes the
35 * current result into the vector x. The parameter <i>err</i> defines the error
36 * 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 unsigned int i, j, center, finished = 0;
52
53 int k;
54
55 double *Enew;
56
57 double E, err = 0;
58
59 assert(rows >= 0);
60
61 Enew = G_alloc_vector(rows);
62
63 for (k = 0; k < maxit; k++) {
64 err = 0;
65 {
66 if (k == 0) {
67 for (j = 0; j < (unsigned int)rows; j++) {
68 Enew[j] = x[j];
69 }
70 }
71 for (i = 0; i < (unsigned int)rows; i++) {
72 E = 0;
73 center = 0;
74 for (j = 0; j < Asp[i]->cols; j++) {
75 E += Asp[i]->values[j] * x[Asp[i]->index[j]];
76 if (Asp[i]->index[j] == i)
77 center = j;
78 }
79 Enew[i] = x[i] - sor * (E - b[i]) / Asp[i]->values[center];
80 }
81 for (j = 0; j < (unsigned int)rows; j++) {
82 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
83
84 x[j] = Enew[j];
85 }
86 }
87
88 G_message(_("sparse Jacobi -- iteration %5i error %g\n"), k, err);
89
90 if (err < error) {
91 finished = 1;
92 break;
93 }
94 }
95
96 G_free(Enew);
97
98 return finished;
99}
100
101/*!
102 * \brief The iterative gauss seidel solver for sparse matrices
103 *
104 * The Jacobi solver solves the linear equation system Ax = b
105 * The result is written to the vector x.
106 *
107 * The parameter <i>maxit</i> specifies the maximum number of iterations. If the
108 * maximum is reached, the solver will abort the calculation and writes the
109 * current result into the vector x. The parameter <i>err</i> defines the error
110 * break criteria for the solver.
111 *
112 * \param Asp G_math_spvector ** -- the sparse matrix
113 * \param x double * -- the vector of unknowns
114 * \param b double * -- the right side vector
115 * \param rows int -- number of rows
116 * \param maxit int -- the maximum number of iterations
117 * \param sor double -- defines the successive overrelaxion parameter [0:2]
118 * \param error double -- defines the error break criteria
119 * \return int -- 1=success, -1=could not solve the les
120 *
121 * */
122int G_math_solver_sparse_gs(G_math_spvector **Asp, double *x, double *b,
123 int rows, int maxit, double sor, double error)
124{
125 unsigned int i, j, finished = 0;
126
127 int k;
128
129 double *Enew;
130
131 double E, err = 0;
132
133 int center;
134
135 assert(rows >= 0);
136
137 Enew = G_alloc_vector(rows);
138
139 for (k = 0; k < maxit; k++) {
140 err = 0;
141 {
142 if (k == 0) {
143 for (j = 0; j < (unsigned int)rows; j++) {
144 Enew[j] = x[j];
145 }
146 }
147 for (i = 0; i < (unsigned int)rows; i++) {
148 E = 0;
149 center = 0;
150 for (j = 0; j < Asp[i]->cols; j++) {
151 E += Asp[i]->values[j] * Enew[Asp[i]->index[j]];
152 if (Asp[i]->index[j] == i)
153 center = j;
154 }
155 Enew[i] = x[i] - sor * (E - b[i]) / Asp[i]->values[center];
156 }
157 for (j = 0; j < (unsigned int)rows; j++) {
158 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
159
160 x[j] = Enew[j];
161 }
162 }
163
164 G_message(_("sparse SOR -- iteration %5i error %g\n"), k, err);
165
166 if (err < error) {
167 finished = 1;
168 break;
169 }
170 }
171
172 G_free(Enew);
173
174 return finished;
175}
176
177/*!
178 * \brief The iterative jacobi solver for quadratic matrices
179 *
180 * The Jacobi solver solves the linear equation system Ax = b
181 * The result is written to the vector x.
182 *
183 * The parameter <i>maxit</i> specifies the maximum number of iterations. If the
184 * maximum is reached, the solver will abort the calculation and writes the
185 * current result into the vector x. The parameter <i>err</i> defines the error
186 * break criteria for the solver.
187 *
188 * \param A double ** -- the dense matrix
189 * \param x double * -- the vector of unknowns
190 * \param b double * -- the right side vector
191 * \param rows int -- number of rows
192 * \param maxit int -- the maximum number of iterations
193 * \param sor double -- defines the successive overrelaxion parameter [0:1]
194 * \param error double -- defines the error break criteria
195 * \return int -- 1=success, -1=could not solve the les
196 *
197 * */
198int G_math_solver_jacobi(double **A, double *x, double *b, int rows, int maxit,
199 double sor, double error)
200{
201 int i, j, k;
202
203 double *Enew;
204
205 double E, err = 0;
206
207 Enew = G_alloc_vector(rows);
208
209 for (j = 0; j < rows; j++) {
210 Enew[j] = x[j];
211 }
212
213 for (k = 0; k < maxit; k++) {
214 for (i = 0; i < rows; i++) {
215 E = 0;
216 for (j = 0; j < rows; j++) {
217 E += A[i][j] * x[j];
218 }
219 Enew[i] = x[i] - sor * (E - b[i]) / A[i][i];
220 }
221 err = 0;
222 for (j = 0; j < rows; j++) {
223 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
224 x[j] = Enew[j];
225 }
226 G_message(_("Jacobi -- iteration %5i error %g\n"), k, err);
227 if (err < error)
228 break;
229 }
230
231 return 1;
232}
233
234/*!
235 * \brief The iterative gauss seidel solver for quadratic matrices
236 *
237 * The Jacobi solver solves the linear equation system Ax = b
238 * The result is written to the vector x.
239 *
240 * The parameter <i>maxit</i> specifies the maximum number of iterations. If the
241 * maximum is reached, the solver will abort the calculation and writes the
242 * current result into the vector x. The parameter <i>err</i> defines the error
243 * break criteria for the solver.
244 *
245 * \param A double ** -- the dense matrix
246 * \param x double * -- the vector of unknowns
247 * \param b double * -- the right side vector
248 * \param rows int -- number of rows
249 * \param maxit int -- the maximum number of iterations
250 * \param sor double -- defines the successive overrelaxion parameter [0:2]
251 * \param error double -- defines the error break criteria
252 * \return int -- 1=success, -1=could not solve the les
253 *
254 * */
255int G_math_solver_gs(double **A, double *x, double *b, int rows, int maxit,
256 double sor, double error)
257{
258 int i, j, k;
259
260 double *Enew;
261
262 double E, err = 0;
263
264 Enew = G_alloc_vector(rows);
265
266 for (j = 0; j < rows; j++) {
267 Enew[j] = x[j];
268 }
269
270 for (k = 0; k < maxit; k++) {
271 for (i = 0; i < rows; i++) {
272 E = 0;
273 for (j = 0; j < rows; j++) {
274 E += A[i][j] * Enew[j];
275 }
276 Enew[i] = x[i] - sor * (E - b[i]) / A[i][i];
277 }
278 err = 0;
279 for (j = 0; j < rows; j++) {
280 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
281 x[j] = Enew[j];
282 }
283 G_message(_("SOR -- iteration %5i error %g\n"), k, err);
284 if (err < error)
285 break;
286 }
287
288 return 1;
289}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:39
double b
void G_message(const char *msg,...)
Print a message to stderr.
Definition gis/error.c:89
#define assert(condition)
Definition lz4.c:393
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)
#define x