1 /* svdecomp - SVD decomposition routines. */
2 /* Taken from Numerical Recipies. */
3 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
4 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
5 /* You may give out copies of this software; for conditions see the */
6 /* file COPYING included with this distribution. */
10 static double PYTHAG(a
, b
)
13 double at
= fabs(a
), bt
= fabs(b
), ct
, result
;
15 if (at
> bt
) { ct
= bt
/ at
; result
= at
* sqrt(1.0 + ct
* ct
); }
16 else if (bt
> 0.0) { ct
= at
/ bt
; result
= bt
* sqrt(1.0 + ct
* ct
); }
21 #define SWAPD(a, b) (temp = (a), (a) = (b), (b) = temp)
24 sort_sv(int m
, int n
, int k
,
25 double **a
, double *w
, double **v
)
30 for (i
= k
; (i
< n
- 1) && (w
[i
] < w
[i
+1]); i
++) {
32 for (j
= 0; j
< m
; j
++) SWAPD(a
[j
][i
], a
[j
][i
+1]);
33 for (j
= 0; j
< n
; j
++) SWAPD(v
[j
][i
], v
[j
][i
+1]);
37 static double maxarg1
, maxarg2
;
38 #define Max(a, b) (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
39 #define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
42 svdcmp(double **a
, int m
, int n
, double *w
, double **v
)
44 int flag
, i
, its
, j
, jj
, k
, l
, nm
;
45 double c
, f
, h
, s
, x
, y
, z
;
46 double anorm
= 0.0, g
= 0.0, scale
= 0.0;
49 if (m
< n
) return(FALSE
); /* flag an error if m < n */
53 /* Householder reduction to bidiagonal form */
54 for (i
= 0; i
< n
; i
++) {
56 /* left-hand reduction */
61 for (k
= i
; k
< m
; k
++) scale
+= fabs(a
[k
][i
]);
63 for (k
= i
; k
< m
; k
++) {
65 s
+= a
[k
][i
] * a
[k
][i
];
68 g
= -SIGN(sqrt(s
), f
);
72 for (j
= l
; j
< n
; j
++) {
73 for (s
= 0.0, k
= i
; k
< m
; k
++) s
+= a
[k
][i
] * a
[k
][j
];
75 for (k
= i
; k
< m
; k
++) a
[k
][j
] += f
* a
[k
][i
];
78 for (k
= i
; k
< m
; k
++) a
[k
][i
] *= scale
;
83 /* right-hand reduction */
85 if (i
< m
&& i
!= n
- 1) {
86 for (k
= l
; k
< n
; k
++) scale
+= fabs(a
[i
][k
]);
88 for (k
= l
; k
< n
; k
++) {
90 s
+= a
[i
][k
] * a
[i
][k
];
93 g
= -SIGN(sqrt(s
), f
);
96 for (k
= l
; k
< n
; k
++) rv1
[k
] = a
[i
][k
] / h
;
98 for (j
= l
; j
< m
; j
++) {
99 for (s
= 0.0, k
= l
; k
< n
; k
++) s
+= a
[j
][k
] * a
[i
][k
];
100 for (k
= l
; k
< n
; k
++) a
[j
][k
] += s
* rv1
[k
];
103 for (k
= l
; k
< n
; k
++) a
[i
][k
] *= scale
;
106 anorm
= Max(anorm
, (fabs(w
[i
]) + fabs(rv1
[i
])));
109 /* accumulate the right-hand transformation */
110 for (i
= n
- 1; i
>= 0; i
--) {
113 for (j
= l
; j
< n
; j
++)
114 v
[j
][i
] = (a
[i
][j
] / a
[i
][l
]) / g
;
115 for (j
= l
; j
< n
; j
++) {
116 for (s
= 0.0, k
= l
; k
< n
; k
++) s
+= a
[i
][k
] * v
[k
][j
];
117 for (k
= l
; k
< n
; k
++) v
[k
][j
] += s
* v
[k
][i
];
120 for (j
= l
; j
< n
; j
++) v
[i
][j
] = v
[j
][i
] = 0.0;
127 /* accumulate the left-hand transformation */
128 for (i
= n
- 1; i
>= 0; i
--) {
132 for (j
= l
; j
< n
; j
++) a
[i
][j
] = 0.0;
136 for (j
= l
; j
< n
; j
++) {
137 for (s
= 0.0, k
= l
; k
< m
; k
++) s
+= a
[k
][i
] * a
[k
][j
];
138 f
= (s
/ a
[i
][i
]) * g
;
139 for (k
= i
; k
< m
; k
++) a
[k
][j
] += f
* a
[k
][i
];
142 for (j
= i
; j
< m
; j
++) a
[j
][i
] *= g
;
145 for (j
= i
; j
< m
; j
++) a
[j
][i
] = 0.0;
150 /* diagonalize the bidiagonal form */
151 for (k
= n
- 1; k
>= 0; k
--) { /* loop over singular values */
152 for (its
= 0; its
< 30; its
++) { /* loop over allowed iterations */
154 for (l
= k
; l
>= 0; l
--) { /* test for splitting */
156 if (fabs(rv1
[l
]) + anorm
== anorm
) {
160 if (fabs(w
[nm
]) + anorm
== anorm
) break;
165 for (i
= l
; i
<= k
; i
++) {
167 if (fabs(f
) + anorm
!= anorm
) {
173 sprintf(s
, "h = %f, f = %f, g = %f\n", h
, f
, g
);
174 fprintf(stdout
,"%s",s
);
179 for (j
= 0; j
< m
; j
++) {
182 a
[j
][nm
] = y
* c
+ z
* s
;
183 a
[j
][i
] = z
* c
- y
* s
;
189 if (l
== k
) { /* convergence */
190 if (z
< 0.0) { /* make singular value nonnegative */
192 for (j
= 0; j
< n
; j
++) v
[j
][k
] = (-v
[j
][k
]);
194 sort_sv(m
, n
, k
, a
, w
, v
);
199 return(FALSE
); /* return an error flag */
202 /* shift from bottom 2 x 2 minor */
208 f
= ((y
- z
) * (y
+ z
) + (g
- h
) * (g
+ h
)) / (2.0 * h
* y
);
210 f
= ((x
- z
) * (x
+ z
) + h
* ((y
/ (f
+ SIGN(g
, f
))) - h
)) / x
;
212 /* next QR transformation */
214 for (j
= l
; j
<= nm
; j
++) {
228 for (jj
= 0; jj
< n
; jj
++) {
231 v
[jj
][j
] = x
* c
+ z
* s
;
232 v
[jj
][i
] = z
* c
- x
* s
;
241 f
= (c
* g
) + (s
* y
);
242 x
= (c
* y
) - (s
* g
);
243 for (jj
= 0; jj
< m
; jj
++) {
246 a
[jj
][j
] = y
* c
+ z
* s
;
247 a
[jj
][i
] = z
* c
- y
* s
;