GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 int G_ludcmp(double **a, int n, int *indx, double *d)
10 {
11  int i, imax = 0, j, k;
12  double big, dum, sum, temp;
13  double *vv;
14 
15  vv = G_alloc_vector(n);
16  *d = 1.0;
17  for (i = 0; i < n; i++) {
18  big = 0.0;
19  for (j = 0; j < n; j++)
20  if ((temp = fabs(a[i][j])) > big)
21  big = temp;
22  if (big == 0.0) {
23  *d = 0.0;
24  return 0; /* Singular matrix */
25  }
26  vv[i] = 1.0 / big;
27  }
28  for (j = 0; j < n; j++) {
29  for (i = 0; i < j; i++) {
30  sum = a[i][j];
31  for (k = 0; k < i; k++)
32  sum -= a[i][k] * a[k][j];
33  a[i][j] = sum;
34  }
35  big = 0.0;
36  for (i = j; i < n; i++) {
37  sum = a[i][j];
38  for (k = 0; k < j; k++)
39  sum -= a[i][k] * a[k][j];
40  a[i][j] = sum;
41  if ((dum = vv[i] * fabs(sum)) >= big) {
42  big = dum;
43  imax = i;
44  }
45  }
46  if (j != imax) {
47  for (k = 0; k < n; k++) {
48  dum = a[imax][k];
49  a[imax][k] = a[j][k];
50  a[j][k] = dum;
51  }
52  *d = -(*d);
53  vv[imax] = vv[j];
54  }
55  indx[j] = imax;
56  if (a[j][j] == 0.0)
57  a[j][j] = TINY;
58  if (j != n) {
59  dum = 1.0 / (a[j][j]);
60  for (i = j + 1; i < n; i++)
61  a[i][j] *= dum;
62  }
63  }
64  G_free_vector(vv);
65 
66  return 1;
67 }
68 
69 #undef TINY
70 
71 void G_lubksb(double **a, int n, int *indx, double b[])
72 {
73  int i, ii, ip, j;
74  double sum;
75 
76  ii = -1;
77  for (i = 0; i < n; i++) {
78  ip = indx[i];
79  sum = b[ip];
80  b[ip] = b[i];
81  if (ii >= 0)
82  for (j = ii; j < i; j++)
83  sum -= a[i][j] * b[j];
84  else if (sum)
85  ii = i;
86  b[i] = sum;
87  }
88  for (i = n - 1; i >= 0; i--) {
89  sum = b[i];
90  for (j = i + 1; j < n; j++)
91  sum -= a[i][j] * b[j];
92  b[i] = sum / a[i][i];
93  }
94 }
float b
Definition: named_colr.c:8
int G_ludcmp(double **a, int n, int *indx, double *d)
Definition: lu.c:9
void G_lubksb(double **a, int n, int *indx, double b[])
Definition: lu.c:71
void G_free_vector(double *v)
Vector memory deallocation.
Definition: dalloc.c:129
#define TINY
Definition: lu.c:6
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
int n
Definition: dataquad.c:291