GRASS Programmer's Manual  6.4.4(2014)-r
N_solvers.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: direkt linear equation system solvers
9 * part of the gpde library
10 *
11 * COPYRIGHT: (C) 2000 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/N_pde.h"
24 #include "solvers_local_proto.h"
25 
26 /*prototypes */
27 static void gauss_elimination(double **A, double *b, int rows);
28 static void lu_decomposition(double **A, int rows);
29 static int cholesky_decomposition(double **A, int rows);
30 static void backward_solving(double **A, double *x, double *b, int rows);
31 static void forward_solving(double **A, double *x, double *b, int rows);
32 
33 /***********************************************************
34  * GAUSS elimination solver for Ax = b *********************
35  * ********************************************************/
47 {
48 
49  if (les->type != N_NORMAL_LES) {
50  G_warning(_("The gauss elimination solver does not work with sparse matrices"));
51  return 0;
52  }
53 
54  if (les->quad != 1) {
55  G_fatal_error(_("The linear equation system is not quadratic"));
56  return 0;
57  }
58 
59 
60  G_message(_("Starting direct gauss elimination solver"));
61 
62  N_les_pivot_create(les);
63  gauss_elimination(les->A, les->b, les->rows);
64  backward_solving(les->A, les->x, les->b, les->rows);
65 
66  return 1;
67 }
68 
69 /***********************************************************
70  * LU solver for Ax = b ************************************
71  * ********************************************************/
82 int N_solver_lu(N_les * les)
83 {
84  int i;
85  double *c, *tmpv;
86 
87  if (les->type != N_NORMAL_LES) {
88  G_warning(_("The lu solver does not work with sparse matrices"));
89  return 0;
90  }
91 
92  if (les->quad != 1) {
93  G_warning(_("The linear equation system is not quadratic"));
94  return -1;
95  }
96 
97 
98  G_message(_("Starting direct lu decomposition solver"));
99 
100  tmpv = vectmem(les->rows);
101  c = vectmem(les->rows);
102 
103  N_les_pivot_create(les);
104  lu_decomposition(les->A, les->rows);
105 
106 #pragma omp parallel
107  {
108 
109 #pragma omp for schedule (static) private(i)
110  for (i = 0; i < les->rows; i++) {
111  tmpv[i] = les->A[i][i];
112  les->A[i][i] = 1;
113  }
114 
115 #pragma omp single
116  {
117  forward_solving(les->A, les->b, les->b, les->rows);
118  }
119 
120 #pragma omp for schedule (static) private(i)
121  for (i = 0; i < les->rows; i++) {
122  les->A[i][i] = tmpv[i];
123  }
124 
125 #pragma omp single
126  {
127  backward_solving(les->A, les->x, les->b, les->rows);
128  }
129  }
130 
131  G_free(c);
132  G_free(tmpv);
133 
134 
135  return 1;
136 }
137 
138 /***********************************************************
139  * cholesky solver for Ax = b ******************************
140  * ********************************************************/
153 {
154  if (les->type != N_NORMAL_LES) {
155  G_warning(_("The cholesky solver does not work with sparse matrices"));
156  return 0;
157  }
158 
159  if (les->quad != 1) {
160  G_warning(_("The linear equation system is not quadratic"));
161  return -1;
162  }
163 
164  /* check for symmetry */
165  if (check_symmetry(les) != 1) {
166  G_warning(_("Matrix is not symmetric!"));
167  return -3;
168  }
169 
170  G_message(_("Starting cholesky decomposition solver"));
171 
172  if (cholesky_decomposition(les->A, les->rows) != 1) {
173  G_warning(_("Unable to solve the linear equation system"));
174  return -2;
175  }
176 
177  forward_solving(les->A, les->b, les->b, les->rows);
178  backward_solving(les->A, les->x, les->b, les->rows);
179 
180  return 1;
181 }
182 
183 
184 /***********************************************************
185  * gauss elimination ***************************************
186  * ********************************************************/
198 void gauss_elimination(double **A, double *b, int rows)
199 {
200  int i, j, k;
201  double tmpval = 0.0;
202 
203  for (k = 0; k < rows - 1; k++) {
204 #pragma omp parallel for schedule (static) private(i, j, tmpval) shared(k, A, b, rows)
205  for (i = k + 1; i < rows; i++) {
206  tmpval = A[i][k] / A[k][k];
207  b[i] = b[i] - tmpval * b[k];
208  for (j = k + 1; j < rows; j++) {
209  A[i][j] = A[i][j] - tmpval * A[k][j];
210  }
211  }
212  }
213 
214  return;
215 }
216 
217 /***********************************************************
218  * lu decomposition ****************************************
219  * ********************************************************/
230 void lu_decomposition(double **A, int rows)
231 {
232 
233  int i, j, k;
234 
235  for (k = 0; k < rows - 1; k++) {
236 #pragma omp parallel for schedule (static) private(i, j) shared(k, A, rows)
237  for (i = k + 1; i < rows; i++) {
238  A[i][k] = A[i][k] / A[k][k];
239  for (j = k + 1; j < rows; j++) {
240  A[i][j] = A[i][j] - A[i][k] * A[k][j];
241  }
242  }
243  }
244 
245  return;
246 }
247 
248 /***********************************************************
249  * cholesky decomposition **********************************
250  * ********************************************************/
262 int cholesky_decomposition(double **A, int rows)
263 {
264 
265  int i, j, k;
266  double sum_1 = 0.0;
267  double sum_2 = 0.0;
268  int error = 0;
269 
270 
271  for (k = 0; k < rows; k++) {
272 #pragma omp parallel private(i, j, sum_2) shared(A, rows, sum_1)
273  {
274 #pragma omp for schedule (static) private(j) reduction(+:sum_1)
275  for (j = 0; j < k; j++) {
276  sum_1 += A[k][j] * A[k][j];
277  }
278 #pragma omp single
279  {
280  if ((A[k][k] - sum_1) < 0) {
281  error++; /*not allowed to exit the OpenMP region */
282  }
283  A[k][k] = sqrt(A[k][k] - sum_1);
284  sum_1 = 0.0;
285  }
286 #pragma omp for schedule (static) private(i, j, sum_2)
287  for (i = k + 1; i < rows; i++) {
288  sum_2 = 0.0;
289  for (j = 0; j < k; j++) {
290  sum_2 += A[i][j] * A[k][j];
291  }
292  A[i][k] = (A[i][k] - sum_2) / A[k][k];
293  }
294  }
295  }
296 
297  /*we need to copy the lower triangle matrix to the upper trianle */
298 #pragma omp parallel for schedule (static) private(i, k) shared(A, rows)
299  for (k = 0; k < rows; k++) {
300  for (i = k + 1; i < rows; i++) {
301  A[k][i] = A[i][k];
302  }
303  }
304 
305  if (error > 0) {
306  G_warning("Matrix is not positive definite");
307  return -1;
308  }
309 
310  return 1;
311 }
312 
313 
314 /***********************************************************
315  * backward solving ****************************************
316  * ********************************************************/
327 void backward_solving(double **A, double *x, double *b, int rows)
328 {
329  int i, j;
330  double tmpval = 0.0;
331 
332  for (i = rows - 1; i >= 0; i--) {
333  tmpval = 0;
334  for (j = i + 1; j < rows; j++) {
335  /*tmpval += A[i][j] * x[j]; */
336  b[i] = b[i] - A[i][j] * x[j];
337  }
338  /*x[i] = (b[i] - tmpval) / A[i][i]; */
339  x[i] = (b[i]) / A[i][i];
340  }
341 
342  return;
343 }
344 
345 
346 /***********************************************************
347  * forward solving *****************************************
348  * ********************************************************/
359 void forward_solving(double **A, double *x, double *b, int rows)
360 {
361  int i, j;
362  double tmpval = 0.0;
363 
364  for (i = 0; i < rows; i++) {
365  tmpval = 0;
366  for (j = 0; j < i; j++) {
367  tmpval += A[i][j] * x[j];
368  }
369  x[i] = (b[i] - tmpval) / A[i][i];
370  }
371 
372  return;
373 }
374 
375 
376 /* ******************************************************* *
377  * ***** solving a tridiagonal eq system ***************** *
378  * ******************************************************* */
379 void thomalg(double **M, double *V, int rows)
380 {
381  double *Vtmp;
382  double *g;
383  double b;
384  int i;
385 
386  Vtmp = vectmem(rows);
387  g = vectmem(rows);
388 
389  for (i = 0; i < rows; i++) {
390  if (i == 0) {
391  b = M[i][i];
392  Vtmp[i] = V[i] / b;
393  }
394  else {
395  b = M[i][i] - M[i][i - 1] * g[i - 1];
396  Vtmp[i] = (V[i] - Vtmp[i - 1] * M[i][i - 1]) / b;
397  }
398  if (i < rows - 1) {
399  g[i] = M[i][i + 1] / b;
400  }
401  }
402 
403  V[rows - 1] = Vtmp[rows - 1];
404  for (i = rows - 2; i >= 0; i--) {
405  V[i] = Vtmp[i] - g[i] * V[i + 1];
406  }
407 
408  G_free(Vtmp);
409  G_free(g);
410 }
411 
412 
413 /***********************************************************
414  * vectmem *************************************************
415  * ********************************************************/
423 double *vectmem(int rows)
424 {
425  double *vector;
426 
427  vector = (double *)(G_calloc(rows, (sizeof(double))));
428  return vector;
429 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
float b
Definition: named_colr.c:8
def error(msg)
Display an error message using g.message -e
Definition: core.py:370
int quad
Definition: N_pde.h:104
int N_les_pivot_create(N_les *les)
Optimize the structure of the linear equation system with a common pivoting strategy.
Definition: N_les_pivot.c:47
double * x
Definition: N_pde.h:98
int G_warning(const char *msg,...)
Print a warning message to stderr.
double ** A
Definition: N_pde.h:100
void G_message(const char *msg,...)
Print a message to stderr.
Definition: lib/gis/error.c:74
double * vectmem(int rows)
Allocate vector memory.
Definition: N_solvers.c:423
int N_solver_gauss(N_les *les)
The gauss elimination solver for quardatic matrices.
Definition: N_solvers.c:46
float g
Definition: named_colr.c:8
int N_solver_cholesky(N_les *les)
The choleksy decomposition solver for quardatic, symmetric positiv definite matrices.
Definition: N_solvers.c:152
double * b
Definition: N_pde.h:99
#define N_NORMAL_LES
Definition: N_pde.h:41
int check_symmetry(N_les *L)
Check if the matrix in les is symmetric.
void thomalg(double **M, double *V, int rows)
Definition: N_solvers.c:379
int G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
int rows
Definition: N_pde.h:102
The linear equation system (les) structure.
Definition: N_pde.h:96
int type
Definition: N_pde.h:105
int N_solver_lu(N_les *les)
The LU solver for quardatic matrices.
Definition: N_solvers.c:82