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