GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
lu.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <grass/gis.h>
3 #include <grass/gmath.h>
4 
5 
6 #define TINY 1.0e-20;
7 
8 /*!
9  * \brief LU decomposition
10  *
11  * \param a double **
12  * \param n int
13  * \param indx int *
14  * \param d double *
15  *
16  * \return 0 on singular matrix, 1 on success
17  */
18 int G_ludcmp(double **a, int n, int *indx, double *d)
19 {
20  int i, imax = 0, j, k;
21  double big, dum, sum, temp;
22  double *vv;
23  int is_singular = FALSE;
24 
25  vv = G_alloc_vector(n);
26  *d = 1.0;
27 /* this pragma works, but doesn't really help speed things up */
28 /* #pragma omp parallel for private(i, j, big, temp) shared(n, a, vv, is_singular) */
29  for (i = 0; i < n; i++) {
30  big = 0.0;
31  for (j = 0; j < n; j++)
32  if ((temp = fabs(a[i][j])) > big)
33  big = temp;
34 
35  if (big == 0.0) {
36  is_singular = TRUE;
37  break;
38  }
39 
40  vv[i] = 1.0 / big;
41  }
42  if (is_singular) {
43  *d = 0.0;
44  return 0; /* Singular matrix */
45  }
46 
47  for (j = 0; j < n; j++) {
48  for (i = 0; i < j; i++) {
49  sum = a[i][j];
50  for (k = 0; k < i; k++)
51  sum -= a[i][k] * a[k][j];
52  a[i][j] = sum;
53  }
54 
55  big = 0.0;
56 /* not very efficient, but this pragma helps speed things up a bit */
57 #pragma omp parallel for private(i, k, sum, dum) shared(j, n, a, vv, big, imax)
58  for (i = j; i < n; i++) {
59  sum = a[i][j];
60  for (k = 0; k < j; k++)
61  sum -= a[i][k] * a[k][j];
62  a[i][j] = sum;
63  if ((dum = vv[i] * fabs(sum)) >= big) {
64  big = dum;
65  imax = i;
66  }
67  }
68  if (j != imax) {
69  for (k = 0; k < n; k++) {
70  dum = a[imax][k];
71  a[imax][k] = a[j][k];
72  a[j][k] = dum;
73  }
74  *d = -(*d);
75  vv[imax] = vv[j];
76  }
77  indx[j] = imax;
78  if (a[j][j] == 0.0)
79  a[j][j] = TINY;
80  if (j != n) {
81  dum = 1.0 / (a[j][j]);
82  for (i = j + 1; i < n; i++)
83  a[i][j] *= dum;
84  }
85  }
86  G_free_vector(vv);
87 
88  return 1;
89 }
90 
91 #undef TINY
92 
93 
94 /*!
95  * \brief LU backward substitution
96  *
97  * \param a double **
98  * \param n int
99  * \param indx int *
100  * \param b double []
101  *
102  * \return void
103  */
104 void G_lubksb(double **a, int n, int *indx, double b[])
105 {
106  int i, ii, ip, j;
107  double sum;
108 
109  ii = -1;
110  for (i = 0; i < n; i++) {
111  ip = indx[i];
112  sum = b[ip];
113  b[ip] = b[i];
114  if (ii >= 0)
115  for (j = ii; j < i; j++)
116  sum -= a[i][j] * b[j];
117  else if (sum)
118  ii = i;
119  b[i] = sum;
120  }
121  for (i = n - 1; i >= 0; i--) {
122  sum = b[i];
123  for (j = i + 1; j < n; j++)
124  sum -= a[i][j] * b[j];
125  b[i] = sum / a[i][i];
126  }
127 }
#define TRUE
Definition: gis.h:59
void G_free_vector(double *)
Vector memory deallocation.
Definition: dalloc.c:129
int G_ludcmp(double **a, int n, int *indx, double *d)
LU decomposition.
Definition: lu.c:18
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition: lu.c:104
double b
Definition: r_raster.c:39
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:41
#define FALSE
Definition: gis.h:63
#define TINY
Definition: lu.c:6