source: trunk/MagicSoft/slalib/svd.c@ 7525

Last change on this file since 7525 was 732, checked in by tbretz, 25 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 11.6 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3double rms ( double a, double b );
4void slaSvd ( int m, int n, int mp, int np, double *a, double *w,
5 double *v, double *work, int *jstat )
6/*
7** - - - - - - -
8** s l a S v d
9** - - - - - - -
10**
11** Singular value decomposition.
12**
13** (double precision)
14**
15** This routine expresses a given matrix a as the product of
16** three matrices u, w, v:
17**
18** a = u x w x vt
19**
20** where:
21**
22** a is any m (rows) x n (columns) matrix, where m >= n
23** u is an m x n column-orthogonal matrix
24** w is an n x n diagonal matrix with w(i,i) >= 0
25** vt is the transpose of an n x n orthogonal matrix
26**
27** Note that m and n, above, are the logical dimensions of the
28** matrices and vectors concerned, which can be located in
29** arrays of larger physical dimensions, given by mp and np.
30**
31** Given:
32** m,n int numbers of rows and columns in matrix a
33** mp,np int physical dimensions of the array containing a
34** a double[mp][np] array containing m x n matrix a
35**
36** Returned:
37** *a double[mp][np] array containing m x n column-orthogonal matrix u
38** *w double[n] n x n diagonal matrix w (diagonal elements only)
39** *v double[np][np] array containing n x n orthogonal matrix v
40** *work double[n] workspace
41** *jstat int 0 = OK
42** -1 = the a array is the wrong shape
43** >0 = 1 + index of w for which convergence failed.
44**
45** (n.b. v contains matrix v, not the transpose of matrix v)
46**
47** References:
48** The algorithm is an adaptation of the routine SVD in the EISPACK
49** library (Garbow et al 1977, Eispack guide extension, Springer
50** Verlag), which is a Fortran 66 implementation of the Algol
51** routine SVD of Wilkinson & Reinsch 1971 (Handbook for Automatic
52** Computation, vol 2, Ed Bauer et al, Springer Verlag). For the
53** non-specialist, probably the clearest general account of the use
54** of SVD in least squares problems is given in Numerical Recipes
55** (Press et al 1986, Cambridge University Press).
56**
57** From slamac.h: TRUE, FALSE
58**
59** Example call (note handling of "adjustable dimension" 2D arrays):
60**
61** double a[MP][NP], w[NP], v[NP][NP], work[NP];
62** int m, n, j;
63** :
64** slaSvd ( m, n, MP, NP, (double *) a, w, (double *) v, work, &j );
65**
66** Last revision: 24 June 1997
67**
68** Copyright P.T.Wallace. All rights reserved.
69*/
70
71/* Maximum number of iterations in QR phase */
72#define ITMAX 30
73{
74
75 int i, k, l, j, k1, its, l1, i1, cancel;
76 double g, scale, an, s, x, f, h, cn, c, y, z;
77 double *ai, *aj, *ak;
78 double *vi, *vj, *vk;
79
80/* Check that the matrix is the right size and shape */
81 if ( m < n || m > mp || n > np ) {
82 *jstat = -1;
83 } else {
84 *jstat = 0;
85
86 /* Householder reduction to bidiagonal form */
87 g = 0.0;
88 scale = 0.0;
89 an = 0.0;
90 for ( i = 0, ai = a; i < n; i++, ai += np ) {
91 l = i + 1;
92 work[i] = scale * g;
93 g = 0.0;
94 s = 0.0;
95 scale = 0.0;
96 if ( i < m ) {
97 for ( k = i, ak = ai; k < m; k++, ak += np ) {
98 scale += fabs ( ak[i] );
99 }
100 if ( scale != 0.0 ) {
101 for ( k = i, ak = ai; k < m; k++, ak += np ) {
102 x = ak[i] / scale;
103 ak[i] = x;
104 s += x * x;
105 }
106 f = ai[i];
107 g = - dsign ( sqrt ( s ), f );
108 h = f * g - s;
109 ai[i] = f - g;
110 if ( i != n - 1 ) {
111 for ( j = l; j < n; j++ ) {
112 s = 0.0;
113 for ( k = i, ak = ai; k < m; k++, ak += np ) {
114 s += ak[i] * ak[j];
115 }
116 f = s / h;
117 for ( k = i, ak = ai; k < m; k++, ak += np ) {
118 ak[j] += f * ak[i];
119 }
120 }
121 }
122 for ( k = i, ak = ai; k < m; k++, ak += np ) {
123 ak[i] *= scale;
124 }
125 }
126 }
127 w[i] = scale * g;
128 g = 0.0;
129 s = 0.0;
130 scale = 0.0;
131 if ( i < m && i != n - 1 ) {
132 for ( k = l; k < n; k++ ) {
133 scale += fabs ( ai[k] );
134 }
135 if ( scale != 0.0 ) {
136 for ( k = l; k < n; k++ ) {
137 x = ai[k] / scale;
138 ai[k] = x;
139 s += x * x;
140 }
141 f = ai[l];
142 g = - dsign ( sqrt ( s ), f );
143 h = f * g - s;
144 ai[l] = f - g;
145 for ( k = l; k < n; k++ ) {
146 work[k] = ai[k] / h;
147 }
148 if ( i != m-1 ) {
149 for ( j = l, aj = a + l*np; j < m; j++, aj += np ) {
150 s = 0.0;
151 for ( k = l; k < n; k++ ) {
152 s += aj[k] * ai[k];
153 }
154 for ( k = l; k < n; k++ ) {
155 aj[k] += s * work[k];
156 }
157 }
158 }
159 for ( k = l; k < n; k++ ) {
160 ai[k] *= scale;
161 }
162 }
163 }
164
165 /* Overestimate of largest column norm for convergence test */
166 cn = fabs ( w[i] ) + fabs ( work[i] );
167 an = gmax ( an, cn );
168 }
169
170 /* Accumulation of right-hand transformations */
171 for ( i = n - 1, ai = a + ( n - 1 ) * np, vi = v + ( n - 1 ) * np;
172 i >= 0;
173 i--, ai -= np, vi -= np ) {
174 if ( i != n - 1 ) {
175 if ( g != 0.0 ) {
176 for ( j = l, vj = v + l * np; j < n; j++, vj += np ) {
177 vj[i] = ( ai[j] / ai[l] ) / g;
178 }
179 for ( j = l; j < n; j++ ) {
180 s = 0.0;
181 for ( k = l, vk = v + l*np; k < n; k++, vk += np ) {
182 s += ai[k] * vk[j];
183 }
184 for ( k = l, vk = v + l*np; k < n; k++, vk += np ) {
185 vk[j] += s * vk[i];
186 }
187 }
188 }
189 for ( j = l, vj = v + l*np; j < n; j++, vj += np ) {
190 vi[j] = 0.0;
191 vj[i] = 0.0;
192 }
193 }
194 vi[i] = 1.0;
195 g = work[i];
196 l = i;
197 }
198
199 /* Accumulation of left-hand transformations */
200 for ( i = n - 1, ai = a + i*np; i >= 0; i--, ai -= np ) {
201 l = i + 1;
202 g = w[i];
203 if ( i != n - 1 ) {
204 for ( j = l; j < n; j++ ) {
205 ai[j] = 0.0;
206 }
207 }
208 if ( g != 0.0 ) {
209 if ( i != n - 1 ) {
210 for ( j = l; j < n; j++ ) {
211 s = 0.0;
212 for ( k = l, ak = a + l * np; k < m; k++, ak += np ) {
213 s += ak[i] * ak[j];
214 }
215 f = ( s / ai[i] ) / g;
216 for ( k = i, ak = a + i * np; k < m; k++, ak += np ) {
217 ak[j] += f * ak[i];
218 }
219 }
220 }
221 for ( j = i, aj = ai; j < m; j++, aj += np ) {
222 aj[i] /= g;
223 }
224 } else {
225 for ( j = i, aj = ai; j < m; j++, aj += np ) {
226 aj[i] = 0.0;
227 }
228 }
229 ai[i] += 1.0;
230 }
231
232 /* Diagonalization of the bidiagonal form */
233 for ( k = n - 1; k >= 0; k-- ) {
234 k1 = k - 1;
235
236 /* Iterate until converged */
237 for ( its = 1; its <= ITMAX; its++ ) {
238
239 /* Test for splitting into submatrices */
240 cancel = TRUE;
241 for ( l = k; l >= 0; l-- ) {
242 l1 = l - 1;
243 if ( an + fabs ( work[l] ) == an ) {
244 cancel = FALSE;
245 break;
246 }
247 /* (Following never attempted for l=0 because work[0] is zero) */
248 if ( an + fabs ( w[l1] ) == an ) {
249 break;
250 }
251 }
252
253 /* Cancellation of work[l] if l>0 */
254 if ( cancel ) {
255 c = 0.0;
256 s = 1.0;
257 for ( i = l; i <= k; i++ ) {
258 f = s * work[i];
259 if ( an + fabs ( f ) == an ) {
260 break;
261 }
262 g = w[i];
263 h = rms ( f, g );
264 w[i] = h;
265 c = g / h;
266 s = - f / h;
267 for ( j = 0, aj = a; j < m; j++, aj += np ) {
268 y = aj[l1];
269 z = aj[i];
270 aj[l1] = y * c + z * s;
271 aj[i] = - y * s + z * c;
272 }
273 }
274 }
275
276 /* Converged? */
277 z = w[k];
278 if ( l == k ) {
279
280 /* Yes: ensure singular values non-negative */
281 if ( z < 0.0 ) {
282 w[k] = -z;
283 for ( j = 0, vj = v; j < n; j++, vj += np ) {
284 vj[k] = -vj[k];
285 }
286 }
287
288 /* Stop iterating */
289 break;
290
291 } else {
292
293 /* Not converged yet: set status if iteration limit reached */
294 if ( its >= ITMAX ) {
295 *jstat = k + 1;
296 }
297
298 /* Shift from bottom 2 x 2 minor */
299 x = w[l];
300 y = w[k1];
301 g = work[k1];
302 h = work[k];
303 f = ( ( y - z ) * ( y + z )
304 + ( g - h ) * ( g + h ) ) / ( 2.0 * h * y );
305 g = ( fabs ( f ) <= 1e15 ) ? rms ( f, 1.0 ) : fabs ( f );
306 f = ( ( x - z ) * ( x + z )
307 + h * ( y / ( f + dsign ( g, f ) ) - h ) ) / x;
308
309 /* Next QR transformation */
310 c = 1.0;
311 s = 1.0;
312 for ( i1 = l; i1 <= k1; i1++ ) {
313 i = i1 + 1;
314 g = work[i];
315 y = w[i];
316 h = s * g;
317 g = c * g;
318 z = rms ( f, h );
319 work[i1] = z;
320 if ( z != 0.0 ) {
321 c = f / z;
322 s = h / z;
323 } else {
324 c = 1.0;
325 s = 0.0;
326 }
327 f = x * c + g * s;
328 g = - x * s + g * c;
329 h = y * s;
330 y = y * c;
331 for ( j = 0, vj = v; j < n; j++, vj += np ) {
332 x = vj[i1];
333 z = vj[i];
334 vj[i1] = x * c + z * s;
335 vj[i] = - x * s + z * c;
336 }
337 z = rms ( f, h );
338 w[i1] = z;
339 if ( z != 0.0 ) {
340 c = f / z;
341 s = h / z;
342 }
343 f = c * g + s * y;
344 x = - s * g + c * y;
345 for ( j = 0, aj = a; j < m; j++, aj += np ) {
346 y = aj[i1];
347 z = aj[i];
348 aj[i1] = y * c + z * s;
349 aj[i] = - y * s + z * c;
350 }
351 }
352 work[l] = 0.0;
353 work[k] = f;
354 w[k] = x;
355 }
356 }
357 }
358 }
359}
360
361double rms ( double a, double b )
362
363/* sqrt(a*a+b*b) with protection against under/overflow */
364
365{
366 double wa, wb, w;
367
368 wa = fabs ( a );
369 wb = fabs ( b );
370
371 if ( wa > wb ) {
372 w = wa;
373 wa = wb;
374 wb = w;
375 }
376
377 if ( wb == 0.0 ) {
378 return 0.0;
379 } else {
380 w = wa / wb;
381 return wb * sqrt ( 1.0 + w * w );
382 }
383}
Note: See TracBrowser for help on using the repository browser.