GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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 */
17int 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 */
103void 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}
void G_free_vector(double *v)
Vector memory deallocation.
Definition dalloc.c:123
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:39
#define TRUE
Definition dbfopen.c:75
#define FALSE
Definition dbfopen.c:74
double b
#define TINY
Definition findzc.c:30
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