GRASS Programmer's Manual  6.4.4(2014)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
eigen_tools.c
Go to the documentation of this file.
1 #include <grass/gis.h>
2 #include <math.h>
3 
4 
5 #define MAX_ITERS 30
6 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
7 
8 
9 int G_tqli(double d[], double e[], int n, double **z)
10 {
11  int m, l, iter, i, k;
12  double s, r, p, g, f, dd, c, b;
13 
14  for (i = 1; i < n; i++)
15  e[i - 1] = e[i];
16  e[n - 1] = 0.0;
17  for (l = 0; l < n; l++) {
18  iter = 0;
19 
20  do {
21  for (m = l; m < n - 1; m++) {
22  dd = fabs(d[m]) + fabs(d[m + 1]);
23  if (fabs(e[m]) + dd == dd)
24  break;
25  }
26 
27  if (m != l) {
28  if (iter++ == MAX_ITERS)
29  return 0; /* Too many iterations in TQLI */
30  g = (d[l + 1] - d[l]) / (2.0 * e[l]);
31  r = sqrt((g * g) + 1.0);
32  g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
33  s = c = 1.0;
34  p = 0.0;
35 
36  for (i = m - 1; i >= l; i--) {
37  f = s * e[i];
38  b = c * e[i];
39 
40  if (fabs(f) >= fabs(g)) {
41  c = g / f;
42  r = sqrt((c * c) + 1.0);
43  e[i + 1] = f * r;
44  c *= (s = 1.0 / r);
45  }
46  else {
47  s = f / g;
48  r = sqrt((s * s) + 1.0);
49  e[i + 1] = g * r;
50  s *= (c = 1.0 / r);
51  }
52 
53  g = d[i + 1] - p;
54  r = (d[i] - g) * s + 2.0 * c * b;
55  p = s * r;
56  d[i + 1] = g + p;
57  g = c * r - b;
58 
59  /* Next loop can be omitted if eigenvectors not wanted */
60  for (k = 0; k < n; k++) {
61  f = z[k][i + 1];
62  z[k][i + 1] = s * z[k][i] + c * f;
63  z[k][i] = c * z[k][i] - s * f;
64  }
65  }
66  d[l] = d[l] - p;
67  e[l] = g;
68  e[m] = 0.0;
69  }
70  } while (m != l);
71  }
72 
73  return 1;
74 }
75 
76 
77 void G_tred2(double **a, int n, double d[], double e[])
78 {
79  int l, k, j, i;
80  double scale, hh, h, g, f;
81 
82  for (i = n - 1; i >= 1; i--) {
83  l = i - 1;
84  h = scale = 0.0;
85 
86  if (l > 0) {
87  for (k = 0; k <= l; k++)
88  scale += fabs(a[i][k]);
89 
90  if (scale == 0.0)
91  e[i] = a[i][l];
92  else {
93  for (k = 0; k <= l; k++) {
94  a[i][k] /= scale;
95  h += a[i][k] * a[i][k];
96  }
97 
98  f = a[i][l];
99  g = f > 0 ? -sqrt(h) : sqrt(h);
100  e[i] = scale * g;
101  h -= f * g;
102  a[i][l] = f - g;
103  f = 0.0;
104 
105  for (j = 0; j <= l; j++) {
106  /* Next statement can be omitted if eigenvectors not wanted */
107  a[j][i] = a[i][j] / h;
108  g = 0.0;
109  for (k = 0; k <= j; k++)
110  g += a[j][k] * a[i][k];
111  for (k = j + 1; k <= l; k++)
112  g += a[k][j] * a[i][k];
113  e[j] = g / h;
114  f += e[j] * a[i][j];
115  }
116 
117  hh = f / (h + h);
118  for (j = 0; j <= l; j++) {
119  f = a[i][j];
120  e[j] = g = e[j] - hh * f;
121 
122  for (k = 0; k <= j; k++)
123  a[j][k] -= (f * e[k] + g * a[i][k]);
124  }
125  }
126  }
127  else
128  e[i] = a[i][l];
129  d[i] = h;
130  }
131 
132  /* Next statement can be omitted if eigenvectors not wanted */
133  d[0] = 0.0;
134  e[0] = 0.0;
135 
136  /* Contents of this loop can be omitted if eigenvectors not
137  wanted except for statement d[i]=a[i][i]; */
138  for (i = 0; i < n; i++) {
139  l = i - 1;
140 
141  if (d[i]) {
142  for (j = 0; j <= l; j++) {
143  g = 0.0;
144  for (k = 0; k <= l; k++)
145  g += a[i][k] * a[k][j];
146  for (k = 0; k <= l; k++)
147  a[k][j] -= g * a[k][i];
148  }
149  }
150 
151  d[i] = a[i][i];
152  a[i][i] = 1.0;
153  for (j = 0; j <= l; j++)
154  a[j][i] = a[i][j] = 0.0;
155  }
156 }
#define SIGN(a, b)
Definition: eigen_tools.c:6
float b
Definition: named_colr.c:8
float r
Definition: named_colr.c:8
float g
Definition: named_colr.c:8
#define MAX_ITERS
Definition: eigen_tools.c:5
int G_tqli(double d[], double e[], int n, double **z)
Definition: eigen_tools.c:9
void G_tred2(double **a, int n, double d[], double e[])
Definition: eigen_tools.c:77
tuple h
panel.defaultSize = wx.CheckBox(panel, id = wx.ID_ANY, label = _("Use default size")) panel...