GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-8828094037
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 #define TINY 1.0e-20;
6 
7 /*!
8  * \brief LU decomposition
9  *
10  * \param a double **
11  * \param n int
12  * \param indx int *
13  * \param d double *
14  *
15  * \return 0 on singular matrix, 1 on success
16  */
17 int G_ludcmp(double **a, int n, int *indx, double *d)
18 {
19  int i, imax = 0, j, k;
20  double big, dum, sum, temp;
21  double *vv;
22  int is_singular = FALSE;
23 
24  vv = G_alloc_vector(n);
25  *d = 1.0;
26  /* this pragma works, but doesn't really help speed things up */
27  /* #pragma omp parallel for private(i, j, big, temp) shared(n, a, vv,
28  * 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  * \brief LU backward substitution
95  *
96  * \param a double **
97  * \param n int
98  * \param indx int *
99  * \param b double []
100  *
101  * \return void
102  */
103 void G_lubksb(double **a, int n, int *indx, double b[])
104 {
105  int i, ii, ip, j;
106  double sum;
107 
108  ii = -1;
109  for (i = 0; i < n; i++) {
110  ip = indx[i];
111  sum = b[ip];
112  b[ip] = b[i];
113  if (ii >= 0)
114  for (j = ii; j < i; j++)
115  sum -= a[i][j] * b[j];
116  else if (sum)
117  ii = i;
118  b[i] = sum;
119  }
120  for (i = n - 1; i >= 0; i--) {
121  sum = b[i];
122  for (j = i + 1; j < n; j++)
123  sum -= a[i][j] * b[j];
124  b[i] = sum / a[i][i];
125  }
126 }
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:39
void G_free_vector(double *)
Vector memory deallocation.
Definition: dalloc.c:123
#define TRUE
Definition: gis.h:79
#define FALSE
Definition: gis.h:83
int G_ludcmp(double **a, int n, int *indx, double *d)
LU decomposition.
Definition: lu.c:17
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition: lu.c:103
#define TINY
Definition: lu.c:5
double b
Definition: r_raster.c:39