GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
solvers_direct.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/gmath.h>
25#include <grass/glocale.h>
26
27#define TINY 1.0e-20
28#define COMP_PIVOT 100
29
30/*!
31 * \brief The gauss elimination solver for quardatic matrices
32 *
33 * This solver does not support sparse matrices
34 * The matrix A will be overwritten.
35 * The result is written to the vector x
36 *
37 * \param A double **
38 * \param x double *
39 * \param b double *
40 * \param rows int
41 * \return int -- 1 success
42 * */
43int G_math_solver_gauss(double **A, double *x, double *b, int rows)
44{
45 G_message(_("Starting direct gauss elimination solver"));
46
49
50 return 1;
51}
52
53/*!
54 * \brief The LU solver for quardatic matrices
55 *
56 * This solver does not support sparse matrices
57 * The matrix A will be overwritten.
58 * The result is written to the vector x in the G_math_les structure
59 *
60 *
61 * \param A double **
62 * \param x double *
63 * \param b double *
64 * \param rows int
65 * \return int -- 1 success
66 * */
67int G_math_solver_lu(double **A, double *x, double *b, int rows)
68{
69 int i;
70
71 double *c, *tmpv;
72
73 G_message(_("Starting direct lu decomposition solver"));
74
75 tmpv = G_alloc_vector(rows);
76 c = G_alloc_vector(rows);
77
78 G_math_lu_decomposition(A, b, rows);
79
80#pragma omp parallel
81 {
82
83#pragma omp for schedule (static) private(i)
84 for (i = 0; i < rows; i++) {
85 tmpv[i] = A[i][i];
86 A[i][i] = 1;
87 }
88
89#pragma omp single
90 {
92 }
93
94#pragma omp for schedule (static) private(i)
95 for (i = 0; i < rows; i++) {
96 A[i][i] = tmpv[i];
97 }
98
99#pragma omp single
100 {
102 }
103 }
104
105 G_free(c);
106 G_free(tmpv);
107
108
109 return 1;
110}
111
112/*!
113 * \brief The choleksy decomposition solver for quardatic, symmetric
114 * positiv definite matrices
115 *
116 * This solver does not support sparse matrices
117 * The matrix A will be overwritten.
118 * The result is written to the vector x
119 *
120 * \param A double **
121 * \param x double *
122 * \param b double *
123 * \param bandwidth int -- the bandwidth of the band matrix, if unsure set to rows
124 * \param rows int
125 * \return int -- 1 success
126 * */
127int G_math_solver_cholesky(double **A, double *x, double *b, int bandwidth,
128 int rows)
129{
130
131 G_message(_("Starting cholesky decomposition solver"));
132
133 if (G_math_cholesky_decomposition(A, rows, bandwidth) != 1) {
134 G_warning(_("Unable to solve the linear equation system"));
135 return -2;
136 }
137
140
141 return 1;
142}
143
144/*!
145 * \brief Gauss elimination
146 *
147 * To run this solver efficiently,
148 * no pivoting is supported.
149 * The matrix will be overwritten with the decomposite form
150 * \param A double **
151 * \param b double *
152 * \param rows int
153 * \return void
154 *
155 * */
156void G_math_gauss_elimination(double **A, double *b, int rows)
157{
158 int i, j, k;
159
160 double tmpval = 0.0;
161
162 for (k = 0; k < rows - 1; k++) {
163#pragma omp parallel for schedule (static) private(i, j, tmpval) shared(k, A, b, rows)
164 for (i = k + 1; i < rows; i++) {
165 tmpval = A[i][k] / A[k][k];
166 b[i] = b[i] - tmpval * b[k];
167 for (j = k + 1; j < rows; j++) {
168 A[i][j] = A[i][j] - tmpval * A[k][j];
169 }
170 }
171 }
172
173 return;
174}
175
176/*!
177 * \brief lu decomposition
178 *
179 * To run this solver efficiently,
180 * no pivoting is supported.
181 * The matrix will be overwritten with the decomposite form
182 *
183 * \param A double **
184 * \param b double * -- this vector is needed if its part of the linear equation system, otherwise set it to NULL
185 * \param rows int
186 * \return void
187 *
188 * */
189void G_math_lu_decomposition(double **A, double *b, int rows)
190{
191
192 int i, j, k;
193
194 for (k = 0; k < rows - 1; k++) {
195#pragma omp parallel for schedule (static) private(i, j) shared(k, A, rows)
196 for (i = k + 1; i < rows; i++) {
197 A[i][k] = A[i][k] / A[k][k];
198 for (j = k + 1; j < rows; j++) {
199 A[i][j] = A[i][j] - A[i][k] * A[k][j];
200 }
201 }
202 }
203
204 return;
205}
206
207/*!
208 * \brief cholesky decomposition for symmetric, positiv definite matrices
209 * with bandwidth optimization
210 *
211 * The provided matrix will be overwritten with the lower and
212 * upper triangle matrix A = LL^T
213 *
214 * \param A double **
215 * \param rows int
216 * \param bandwidth int -- the bandwidth of the matrix (0 > bandwidth <= cols)
217 * \return void
218 *
219 * */
220int G_math_cholesky_decomposition(double **A, int rows, int bandwidth)
221{
222
223 int i = 0, j = 0, k = 0;
224
225 double sum_1 = 0.0;
226
227 double sum_2 = 0.0;
228
229 int colsize;
230
231 if (bandwidth <= 0)
232 bandwidth = rows;
233
234 colsize = bandwidth;
235
236 for (k = 0; k < rows; k++) {
237#pragma omp parallel for schedule (static) private(i, j, sum_2) shared(A, k) reduction(+:sum_1)
238 for (j = 0; j < k; j++) {
239 sum_1 += A[k][j] * A[k][j];
240 }
241
242 if (0 > (A[k][k] - sum_1)) {
243 G_warning("Matrix is not positive definite. break.");
244 return -1;
245 }
246 A[k][k] = sqrt(A[k][k] - sum_1);
247 sum_1 = 0.0;
248
249 if ((k + bandwidth) > rows) {
250 colsize = rows;
251 }
252 else {
253 colsize = k + bandwidth;
254 }
255
256#pragma omp parallel for schedule (static) private(i, j, sum_2) shared(A, k, sum_1, colsize)
257
258 for (i = k + 1; i < colsize; i++) {
259 sum_2 = 0.0;
260 for (j = 0; j < k; j++) {
261 sum_2 += A[i][j] * A[k][j];
262 }
263 A[i][k] = (A[i][k] - sum_2) / A[k][k];
264 }
265
266 }
267 /* we need to copy the lower triangle matrix to the upper triangle */
268#pragma omp parallel for schedule (static) private(i, k) shared(A, rows)
269 for (k = 0; k < rows; k++) {
270 for (i = k + 1; i < rows; i++) {
271 A[k][i] = A[i][k];
272 }
273 }
274
275
276 return 1;
277}
278
279/*!
280 * \brief backward substitution
281 *
282 * \param A double **
283 * \param x double *
284 * \param b double *
285 * \param rows int
286 * \return void
287 *
288 * */
289void G_math_backward_substitution(double **A, double *x, double *b, int rows)
290{
291 int i, j;
292
293 for (i = rows - 1; i >= 0; i--) {
294 for (j = i + 1; j < rows; j++) {
295 b[i] = b[i] - A[i][j] * x[j];
296 }
297 x[i] = (b[i]) / A[i][i];
298 }
299
300 return;
301}
302
303/*!
304 * \brief forward substitution
305 *
306 * \param A double **
307 * \param x double *
308 * \param b double *
309 * \param rows int
310 * \return void
311 *
312 * */
313void G_math_forward_substitution(double **A, double *x, double *b, int rows)
314{
315 int i, j;
316
317 double tmpval = 0.0;
318
319 for (i = 0; i < rows; i++) {
320 tmpval = 0;
321 for (j = 0; j < i; j++) {
322 tmpval += A[i][j] * x[j];
323 }
324 x[i] = (b[i] - tmpval) / A[i][i];
325 }
326
327 return;
328}
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
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int G_math_solver_lu(double **A, double *x, double *b, int rows)
The LU solver for quardatic matrices.
void G_math_lu_decomposition(double **A, double *b, int rows)
lu decomposition
int G_math_solver_gauss(double **A, double *x, double *b, int rows)
The gauss elimination solver for quardatic matrices.
void G_math_forward_substitution(double **A, double *x, double *b, int rows)
forward substitution
int G_math_solver_cholesky(double **A, double *x, double *b, int bandwidth, int rows)
The choleksy decomposition solver for quardatic, symmetric positiv definite matrices.
void G_math_backward_substitution(double **A, double *x, double *b, int rows)
backward substitution
void G_math_gauss_elimination(double **A, double *b, int rows)
Gauss elimination.
int G_math_cholesky_decomposition(double **A, int rows, int bandwidth)
cholesky decomposition for symmetric, positiv definite matrices with bandwidth optimization
#define x