GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72721
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  * */
43 int G_math_solver_gauss(double **A, double *x, double *b, int rows)
44 {
45  G_message(_("Starting direct gauss elimination solver"));
46
47  G_math_gauss_elimination(A, b, rows);
48  G_math_backward_substitution(A, x, b, rows);
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  * */
67 int 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  {
91  G_math_forward_substitution(A, b, b, rows);
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  {
101  G_math_backward_substitution(A, x, b, rows);
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  * */
127 int 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
138  G_math_forward_substitution(A, b, b, rows);
139  G_math_backward_substitution(A, x, b, rows);
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  * */
156 void 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  * */
189 void 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  * */
220 int 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  * */
289 void 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  * */
313 void 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_math_backward_substitution(double **A, double *x, double *b, int rows)
backward substitution
int G_math_solver_lu(double **A, double *x, double *b, int rows)
The LU solver for quardatic matrices.
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
int G_math_cholesky_decomposition(double **A, int rows, int bandwidth)
cholesky decomposition for symmetric, positiv definite matrices with bandwidth optimization ...
#define x
int G_math_solver_gauss(double **A, double *x, double *b, int rows)
The gauss elimination solver for quardatic matrices.
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_lu_decomposition(double **A, double *b, int rows)
lu decomposition
void G_math_forward_substitution(double **A, double *x, double *b, int rows)
forward substitution
double b
Definition: r_raster.c:39
#define _(str)
Definition: glocale.h:13
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
void G_math_gauss_elimination(double **A, double *b, int rows)
Gauss elimination.
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204