GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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 */
18int 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 */
104void 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}
void G_free_vector(double *v)
Vector memory deallocation.
Definition: dalloc.c:129
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
#define TRUE
Definition: dbfopen.c:183
#define FALSE
Definition: dbfopen.c:182
double b
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
#define TINY
Definition: lu.c:6