1 /**@(#)SVD
2 **/
3
4 #include <stdio.h>
5 #include <math.h>
6 #include <tina/sys.h>
7 #include <tina/sysfuncs.h>
8 #include <tina/math.h>
9
10 /* smc 2/2/93 */
11 #ifdef Print
12 #undef Print
13 #endif
14 /* smc 2/2/93 */
15
16 #define Print(x) { if(out_text!=NULL) out_text(x); }
17
18 /* This macro is now redundant as function exists but kept incase
19 * function is bad re-implementation of macro */
20 #define Radius(u,v) ( (Au=Abs((u))) > (Av=Abs((v))) ? \
21 (Aw=Av/Au, Au*sqrt(1.+Aw*Aw)) : \
22 (Av ? (Aw=Au/Av, Av*sqrt(1.+Aw*Aw)) : 0.0) )
23
24 #define Sign(u,v) ( (v)>=0.0 ? Abs(u) : -Abs(u) )
25
26 #define Max(u,v) ( (u)>(v)? (u) : (v) )
27
28 #define Abs(x) ( (x)>0.0? (x) : (-(x)) )
29
30 static void (*out_text) () = NULL;
31
32 static double radius(double u, double v)
33 {
34 double Au, Av, Aw;
35
36 Au = Abs(u);
37 Av = Abs(v);
38
39 if (Au > Av)
40 {
41 Aw = Av / Au;
42 return (Au * sqrt(1. + Aw * Aw));
43 } else
44 {
45 if (Av)
46 {
47 Aw = Au / Av;
48 return (Av * sqrt(1. + Aw * Aw));
49 } else
50 return 0.0;
51 }
52 }
53
54
55
56 void svd_err_func(void (*text) ( /* ??? */ ))
57 {
58 out_text = text;
59 }
60
61 /*************************** SVDcmp *****************************************
62 * Given matrix A[m][n], m>=n, using svd decomposition A = U W V' to get *
63 * U[m][n], W[n][n] and V[n][n], where U occupies the position of A. *
64 * NOTE: if m<n, A should be filled up to square with zero rows. *
65 * A[m][n] has been destroyed by U[m][n] after the decomposition. *
66 ****************************************************************************/
67 void SVDcmp(double **a, int m, int n, double *w, double **v)
68 {
69 /* BUG `nm' may be used uninitialized in this function */
70 int flag, i, its, j, jj, k, l, nm, nm1 = n - 1, mm1 = m - 1;
71 double c, f, h, s, x, y, z;
72 double anorm = 0.0, g = 0.0, scale = 0.0;
73 double *rv1;
74
75 if (m < n)
76 Print("SVDCMP: You must augment A with extra zero rows");
77 rv1 = (double *) ralloc((unsigned) n * sizeof(double));
78
79 /* Householder reduction to bidigonal form */
80 for (i = 0; i < n; i++)
81 {
82 l = i + 1;
83 rv1[i] = scale * g;
84 g = s = scale = 0.0;
85 if (i < m)
86 {
87 for (k = i; k < m; k++)
88 scale += Abs(a[k][i]);
89 if (scale)
90 {
91 for (k = i; k < m; k++)
92 {
93 a[k][i] /= scale;
94 s += a[k][i] * a[k][i];
95 }
96 f = a[i][i];
97 g = -Sign(sqrt(s), f);
98 h = f * g - s;
99 a[i][i] = f - g;
100 if (i != nm1)
101 {
102 for (j = l; j < n; j++)
103 {
104 for (s = 0.0, k = i; k < m; k++)
105 s += a[k][i] * a[k][j];
106 f = s / h;
107 for (k = i; k < m; k++)
108 a[k][j] += f * a[k][i];
109 }
110 }
111 for (k = i; k < m; k++)
112 a[k][i] *= scale;
113 }
114 }
115 w[i] = scale * g;
116 g = s = scale = 0.0;
117 if (i < m && i != nm1)
118 {
119 for (k = l; k < n; k++)
120 scale += Abs(a[i][k]);
121 if (scale)
122 {
123 for (k = l; k < n; k++)
124 {
125 a[i][k] /= scale;
126 s += a[i][k] * a[i][k];
127 }
128 f = a[i][l];
129 g = -Sign(sqrt(s), f);
130 h = f * g - s;
131 a[i][l] = f - g;
132 for (k = l; k < n; k++)
133 rv1[k] = a[i][k] / h;
134 if (i != mm1)
135 {
136 for (j = l; j < m; j++)
137 {
138 for (s = 0.0, k = l; k < n; k++)
139 s += a[j][k] * a[i][k];
140 for (k = l; k < n; k++)
141 a[j][k] += s * rv1[k];
142 }
143 }
144 for (k = l; k < n; k++)
145 a[i][k] *= scale;
146 }
147 }
148 anorm = Max(anorm, (Abs(w[i]) + Abs(rv1[i])));
149 }
150
151 /* Accumulation of right-hand transformations */
152 for (i = n - 1; i >= 0; i--)
153 {
154 if (i < nm1)
155 {
156 if (g)
157 {
158 /* double division to avoid possible underflow */
159 for (j = l; j < n; j++)
160 v[j][i] = (a[i][j] / a[i][l]) / g;
161 for (j = l; j < n; j++)
162 {
163 for (s = 0.0, k = l; k < n; k++)
164 s += a[i][k] * v[k][j];
165 for (k = l; k < n; k++)
166 v[k][j] += s * v[k][i];
167 }
168 }
169 for (j = l; j < n; j++)
170 v[i][j] = v[j][i] = 0.0;
171 }
172 v[i][i] = 1.0;
173 g = rv1[i];
174 l = i;
175 }
176 /* Accumulation of left-hand transformations */
177 for (i = n - 1; i >= 0; i--)
178 {
179 l = i + 1;
180 g = w[i];
181 if (i < nm1)
182 for (j = l; j < n; j++)
183 a[i][j] = 0.0;
184 if (g)
185 {
186 g = 1.0 / g;
187 if (i != nm1)
188 {
189 for (j = l; j < n; j++)
190 {
191 for (s = 0.0, k = l; k < m; k++)
192 s += a[k][i] * a[k][j];
193 f = (s / a[i][i]) * g;
194 for (k = i; k < m; k++)
195 a[k][j] += f * a[k][i];
196 }
197 }
198 for (j = i; j < m; j++)
199 a[j][i] *= g;
200 } else
201 for (j = i; j < m; j++)
202 a[j][i] = 0.0;
203 ++a[i][i];
204 }
205 /* diagonalization of the bidigonal form */
206 for (k = n - 1; k >= 0; k--)
207 { /* loop over singlar values */
208 for (its = 0; its < 30; its++)
209 { /* loop over allowed iterations */
210 flag = 1;
211 for (l = k; l >= 0; l--)
212 { /* test for splitting */
213 nm = l - 1; /* note that rv1[l] is always zero */
214 if (Abs(rv1[l]) + anorm == anorm)
215 {
216 flag = 0;
217 break;
218 }
219 if (Abs(w[nm]) + anorm == anorm)
220 break;
221 }
222 if (flag)
223 {
224 c = 0.0; /* cancellation of rv1[l], if l>1 */
225 s = 1.0;
226 for (i = l; i <= k; i++)
227 {
228 f = s * rv1[i];
229 if (Abs(f) + anorm != anorm)
230 {
231 g = w[i];
232 h = radius(f, g);
233 w[i] = h;
234 h = 1.0 / h;
235 c = g * h;
236 s = (-f * h);
237 for (j = 0; j < m; j++)
238 {
239 y = a[j][nm];
240 z = a[j][i];
241 a[j][nm] = y * c + z * s;
242 a[j][i] = z * c - y * s;
243 }
244 }
245 }
246 }
247 z = w[k];
248 if (l == k)
249 { /* convergence */
250 if (z < 0.0)
251 {
252 w[k] = -z;
253 for (j = 0; j < n; j++)
254 v[j][k] = (-v[j][k]);
255 }
256 break;
257 }
258 if (its == 30)
259 Print("No convergence in 30 SVDCMP iterations");
260 x = w[l]; /* shift from bottom 2-by-2 minor */
261 nm = k - 1;
262 y = w[nm];
263 g = rv1[nm];
264 h = rv1[k];
265 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
266 g = radius(f, 1.0);
267 /* next QR transformation */
268 f = ((x - z) * (x + z) + h * ((y / (f + Sign(g, f))) - h)) / x;
269 c = s = 1.0;
270 for (j = l; j <= nm; j++)
271 {
272 i = j + 1;
273 g = rv1[i];
274 y = w[i];
275 h = s * g;
276 g = c * g;
277 z = radius(f, h);
278 rv1[j] = z;
279 c = f / z;
280 s = h / z;
281 f = x * c + g * s;
282 g = g * c - x * s;
283 h = y * s;
284 y = y * c;
285 for (jj = 0; jj < n; jj++)
286 {
287 x = v[jj][j];
288 z = v[jj][i];
289 v[jj][j] = x * c + z * s;
290 v[jj][i] = z * c - x * s;
291 }
292 z = radius(f, h);
293 w[j] = z; /* rotation can be arbitrary id z=0 */
294 if (z)
295 {
296 z = 1.0 / z;
297 c = f * z;
298 s = h * z;
299 }
300 f = (c * g) + (s * y);
301 x = (c * y) - (s * g);
302 for (jj = 0; jj < m; jj++)
303 {
304 y = a[jj][j];
305 z = a[jj][i];
306 a[jj][j] = y * c + z * s;
307 a[jj][i] = z * c - y * s;
308 }
309 }
310 rv1[l] = 0.0;
311 rv1[k] = f;
312 w[k] = x;
313 }
314 }
315 rfree((void *) rv1);
316 }
317
318 /******************************* SVDbksb ****************************
319 * Given A[m][n], b[m], solves A x = b in the svd form U W V'x = b *
320 * so x = V U'b/W *
321 * No input quantities are destroyed, so sequential calling is OK *
322 ********************************************************************/
323 void SVDbksb(double **u, double *w, double **v, int m, int n, double *b, double *x)
324 {
325 int jj, j, i;
326 double s, *tmp;
327
328 tmp = (double *) ralloc((unsigned) n * sizeof(double));
329 for (j = 0; j < n; j++)
330 { /* calculate <U,b> */
331 s = 0.0;
332 if (w[j])
333 { /* nonzero result only if nonzero w[j] */
334 for (i = 0; i < m; i++)
335 s += u[i][j] * b[i];
336 s /= w[j]; /* premultiplied by inverse W */
337 }
338 tmp[j] = s;
339 }
340 for (j = 0; j < n; j++)
341 { /* V(U'b/W) */
342 s = 0.0;
343 for (jj = 0; jj < n; jj++)
344 s += v[j][jj] * tmp[jj];
345 x[j] = s;
346 }
347 rfree((void *) tmp);
348 }
349
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.