GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-d6dec75dd4
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  * */
42 int G_math_solver_gauss(double **A, double *x, double *b, int rows)
43 {
44  G_message(_("Starting direct gauss elimination solver"));
45 
46  G_math_gauss_elimination(A, b, rows);
47  G_math_backward_substitution(A, x, b, rows);
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  * */
66 int 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  {
90  G_math_forward_substitution(A, b, b, rows);
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  {
100  G_math_backward_substitution(A, x, b, rows);
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  * */
124 int 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 
135  G_math_forward_substitution(A, b, b, rows);
136  G_math_backward_substitution(A, x, b, rows);
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  * */
153 void 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  * */
186 void 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  * */
217 int 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) \
235  reduction(+ : sum_1)
236  for (j = 0; j < k; j++) {
237  sum_1 += A[k][j] * A[k][j];
238  }
239 
240  if (0 > (A[k][k] - sum_1)) {
241  G_warning("Matrix is not positive definite. break.");
242  return -1;
243  }
244  A[k][k] = sqrt(A[k][k] - sum_1);
245  sum_1 = 0.0;
246 
247  if ((k + bandwidth) > rows) {
248  colsize = rows;
249  }
250  else {
251  colsize = k + bandwidth;
252  }
253 
254 #pragma omp parallel for schedule(static) private(i, j, sum_2) \
255  shared(A, k, sum_1, colsize)
256 
257  for (i = k + 1; i < colsize; i++) {
258  sum_2 = 0.0;
259  for (j = 0; j < k; j++) {
260  sum_2 += A[i][j] * A[k][j];
261  }
262  A[i][k] = (A[i][k] - sum_2) / A[k][k];
263  }
264  }
265  /* we need to copy the lower triangle matrix to the upper triangle */
266 #pragma omp parallel for schedule(static) private(i, k) shared(A, rows)
267  for (k = 0; k < rows; k++) {
268  for (i = k + 1; i < rows; i++) {
269  A[k][i] = A[i][k];
270  }
271  }
272 
273  return 1;
274 }
275 
276 /*!
277  * \brief backward substitution
278  *
279  * \param A double **
280  * \param x double *
281  * \param b double *
282  * \param rows int
283  * \return void
284  *
285  * */
286 void G_math_backward_substitution(double **A, double *x, double *b, int rows)
287 {
288  int i, j;
289 
290  for (i = rows - 1; i >= 0; i--) {
291  for (j = i + 1; j < rows; j++) {
292  b[i] = b[i] - A[i][j] * x[j];
293  }
294  x[i] = (b[i]) / A[i][i];
295  }
296 
297  return;
298 }
299 
300 /*!
301  * \brief forward substitution
302  *
303  * \param A double **
304  * \param x double *
305  * \param b double *
306  * \param rows int
307  * \return void
308  *
309  * */
310 void G_math_forward_substitution(double **A, double *x, double *b, int rows)
311 {
312  int i, j;
313 
314  double tmpval = 0.0;
315 
316  for (i = 0; i < rows; i++) {
317  tmpval = 0;
318  for (j = 0; j < i; j++) {
319  tmpval += A[i][j] * x[j];
320  }
321  x[i] = (b[i] - tmpval) / A[i][i];
322  }
323 
324  return;
325 }
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
void G_warning(const char *,...) __attribute__((format(printf
void G_message(const char *,...) __attribute__((format(printf
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:39
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition: gis.h:47
#define _(str)
Definition: glocale.h:10
double b
Definition: r_raster.c:39
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