GRASS GIS 7 Programmer's Manual
7.8.4(2020)-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
*/
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
}
G_lubksb
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition:
lu.c:104
FALSE
#define FALSE
Definition:
dbfopen.c:182
TINY
#define TINY
Definition:
lu.c:6
G_free_vector
void G_free_vector(double *v)
Vector memory deallocation.
Definition:
dalloc.c:129
TRUE
#define TRUE
Definition:
dbfopen.c:183
b
double b
Definition:
driver/set_window.c:5
G_ludcmp
int G_ludcmp(double **a, int n, int *indx, double *d)
LU decomposition.
Definition:
lu.c:18
G_alloc_vector
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition:
dalloc.c:41
gmath
lu.c
Generated on Mon Oct 5 2020 08:56:03 for GRASS GIS 7 Programmer's Manual by
1.8.20