GRASS Programmer's Manual  6.4.4(2014)-r
svd.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <grass/gis.h>
3 #include <grass/gmath.h>
4 
5 static double at, bt, ct;
6 
7 #define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
8  (ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
9 
10 static double maxarg1, maxarg2;
11 
12 #define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
13  (maxarg1) : (maxarg2))
14 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
15 
16 int G_svdcmp(double **a, int m, int n, double *w, double **v)
17 {
18  int flag, i, its, j, jj, k, ii = 0, nm = 0;
19  double c, f, h, s, x, y, z;
20  double anorm = 0.0, g = 0.0, scale = 0.0;
21  double *rv1, *G_alloc_vector();
22 
23 
24  if (m < n)
25  return -1; /* must augment A with extra zero rows */
26  rv1 = G_alloc_vector(n);
27 
28  n--;
29  m--;
30 
31  for (i = 0; i <= n; i++) {
32  ii = i + 1;
33  rv1[i] = scale * g;
34  g = s = scale = 0.0;
35  if (i <= m) {
36  for (k = i; k <= m; k++)
37  scale += fabs(a[k][i]);
38  if (scale) {
39  for (k = i; k <= m; k++) {
40  a[k][i] /= scale;
41  s += a[k][i] * a[k][i];
42  }
43  f = a[i][i];
44  g = -SIGN(sqrt(s), f);
45  h = f * g - s;
46  a[i][i] = f - g;
47  if (i != n) {
48  for (j = ii; j <= n; j++) {
49  for (s = 0.0, k = i; k <= m; k++)
50  s += a[k][i] * a[k][j];
51  f = s / h;
52  for (k = i; k <= m; k++)
53  a[k][j] += f * a[k][i];
54  }
55  }
56  for (k = i; k <= m; k++)
57  a[k][i] *= scale;
58  }
59  }
60  w[i] = scale * g;
61  g = s = scale = 0.0;
62  if (i <= m && i != n) {
63  for (k = ii; k <= n; k++)
64  scale += fabs(a[i][k]);
65  if (scale) {
66  for (k = ii; k <= n; k++) {
67  a[i][k] /= scale;
68  s += a[i][k] * a[i][k];
69  }
70  f = a[i][ii];
71  g = -SIGN(sqrt(s), f);
72  h = f * g - s;
73  a[i][ii] = f - g;
74  for (k = ii; k <= n; k++)
75  rv1[k] = a[i][k] / h;
76  if (i != m) {
77  for (j = ii; j <= m; j++) {
78  for (s = 0.0, k = ii; k <= n; k++)
79  s += a[j][k] * a[i][k];
80  for (k = ii; k <= n; k++)
81  a[j][k] += s * rv1[k];
82  }
83  }
84  for (k = ii; k <= n; k++)
85  a[i][k] *= scale;
86  }
87  }
88  anorm = MAX(anorm, (fabs(w[i]) + fabs(rv1[i])));
89  }
90  for (i = n; i >= 0; i--) {
91  if (i < n) {
92  if (g) {
93  for (j = ii; j <= n; j++)
94  v[j][i] = (a[i][j] / a[i][ii]) / g;
95  for (j = ii; j <= n; j++) {
96  for (s = 0.0, k = ii; k <= n; k++)
97  s += a[i][k] * v[k][j];
98  for (k = ii; k <= n; k++)
99  v[k][j] += s * v[k][i];
100  }
101  }
102  for (j = ii; j <= n; j++)
103  v[i][j] = v[j][i] = 0.0;
104  }
105  v[i][i] = 1.0;
106  g = rv1[i];
107  ii = i;
108  }
109  for (i = n; i >= 0; i--) {
110  ii = i + 1;
111  g = w[i];
112  if (i < n)
113  for (j = ii; j <= n; j++)
114  a[i][j] = 0.0;
115  if (g) {
116  g = 1.0 / g;
117  if (i != n) {
118  for (j = ii; j <= n; j++) {
119  for (s = 0.0, k = ii; k <= m; k++)
120  s += a[k][i] * a[k][j];
121  f = (s / a[i][i]) * g;
122  for (k = i; k <= m; k++)
123  a[k][j] += f * a[k][i];
124  }
125  }
126  for (j = i; j <= m; j++)
127  a[j][i] *= g;
128  }
129  else {
130  for (j = i; j <= m; j++)
131  a[j][i] = 0.0;
132  }
133  ++a[i][i];
134  }
135  for (k = n; k >= 0; k--) {
136  for (its = 1; its <= 30; its++) {
137  flag = 1;
138  for (ii = k; ii >= 0; ii--) {
139  nm = ii - 1;
140  if (fabs(rv1[ii]) + anorm == anorm) {
141  flag = 0;
142  break;
143  }
144  if (fabs(w[nm]) + anorm == anorm)
145  break;
146  }
147  if (flag) {
148  c = 0.0;
149  s = 1.0;
150  for (i = ii; i <= k; i++) {
151  f = s * rv1[i];
152  if (fabs(f) + anorm != anorm) {
153  g = w[i];
154  h = PYTHAG(f, g);
155  w[i] = h;
156  h = 1.0 / h;
157  c = g * h;
158  s = (-f * h);
159  for (j = 0; j <= m; j++) {
160  y = a[j][nm];
161  z = a[j][i];
162  a[j][nm] = y * c + z * s;
163  a[j][i] = z * c - y * s;
164  }
165  }
166  }
167  }
168  z = w[k];
169  if (ii == k) {
170  if (z < 0.0) {
171  w[k] = -z;
172  for (j = 0; j <= n; j++)
173  v[j][k] = (-v[j][k]);
174  }
175  break;
176  }
177  if (its == 30)
178  return -2; /*No convergence in 30 SVDCMP iterations */
179  x = w[ii];
180  nm = k - 1;
181  y = w[nm];
182  g = rv1[nm];
183  h = rv1[k];
184  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
185  g = PYTHAG(f, 1.0);
186  f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
187  c = s = 1.0;
188  for (j = ii; j <= nm; j++) {
189  i = j + 1;
190  g = rv1[i];
191  y = w[i];
192  h = s * g;
193  g = c * g;
194  z = PYTHAG(f, h);
195  rv1[j] = z;
196  c = f / z;
197  s = h / z;
198  f = x * c + g * s;
199  g = g * c - x * s;
200  h = y * s;
201  y = y * c;
202  for (jj = 0; jj <= n; jj++) {
203  x = v[jj][j];
204  z = v[jj][i];
205  v[jj][j] = x * c + z * s;
206  v[jj][i] = z * c - x * s;
207  }
208  z = PYTHAG(f, h);
209  w[j] = z;
210  if (z) {
211  z = 1.0 / z;
212  c = f * z;
213  s = h * z;
214  }
215  f = (c * g) + (s * y);
216  x = (c * y) - (s * g);
217  for (jj = 0; jj <= m; jj++) {
218  y = a[jj][j];
219  z = a[jj][i];
220  a[jj][j] = y * c + z * s;
221  a[jj][i] = z * c - y * s;
222  }
223  }
224  rv1[ii] = 0.0;
225  rv1[k] = f;
226  w[k] = x;
227  }
228  }
229  G_free_vector(rv1);
230  return 0;
231 }
232 
233 #undef SIGN
234 #undef MAX
235 #undef PYTHAG
236 
237 int G_svbksb(double **u, double w[], double **v,
238  int m, int n, double b[], double x[])
239 {
240  int j, i;
241  double s, *tmp, *G_alloc_vector();
242 
243  tmp = G_alloc_vector(n);
244  for (j = 0; j < n; j++) {
245  s = 0.0;
246  if (w[j]) {
247  for (i = 0; i < m; i++)
248  s += u[i][j] * b[i];
249  s /= w[j];
250  }
251  tmp[j] = s;
252  }
253  for (j = 0; j < n; j++) {
254  s = 0.0;
255  for (i = 0; i < n; i++)
256  s += v[j][i] * tmp[i];
257  x[j] = s;
258  }
259  G_free_vector(tmp);
260 
261  return 0;
262 }
263 
264 #define TOL 1e-8
265 
266 int G_svelim(double *w, int n)
267 {
268  int i;
269  double thresh;
270 
271  thresh = 0.0; /* remove singularity */
272  for (i = 0; i < n; i++)
273  if (w[i] > thresh)
274  thresh = w[i];
275  thresh *= TOL;
276  for (i = 0; i < n; i++)
277  if (w[i] < thresh)
278  w[i] = 0.0;
279 
280  return 0;
281 }
282 
283 #undef TOL
b
float b
Definition: named_colr.c:8
G_svbksb
int G_svbksb(double **u, double w[], double **v, int m, int n, double b[], double x[])
Definition: svd.c:237
MAX
#define MAX(a, b)
Definition: svd.c:12
y
int y
Definition: plot.c:34
G_svdcmp
int G_svdcmp(double **a, int m, int n, double *w, double **v)
Definition: svd.c:16
G_free_vector
void G_free_vector(double *v)
Vector memory deallocation.
Definition: dalloc.c:129
g
float g
Definition: named_colr.c:8
tools.flag
flag
Definition: tools.py:1403
SIGN
#define SIGN(a, b)
Definition: svd.c:14
TOL
#define TOL
Definition: svd.c:264
G_alloc_vector
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
G_svelim
int G_svelim(double *w, int n)
Definition: svd.c:266
dialogs.h
tuple h
panel.defaultSize = wx.CheckBox(panel, id = wx.ID_ANY, label = _("Use default size")) panel...
Definition: psmap/dialogs.py:2736
PYTHAG
#define PYTHAG(a, b)
Definition: svd.c:7