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. */
12 static double PYTHAG(double a
, double b
)
14 double at
= fabs(a
), bt
= fabs(b
), ct
, result
;
16 if (at
> bt
) { ct
= bt
/ at
; result
= at
* sqrt(1.0 + ct
* ct
); }
17 else if (bt
> 0.0) { ct
= at
/ bt
; result
= bt
* sqrt(1.0 + ct
* ct
); }
22 #define SWAPD(a, b) (temp = (a), (a) = (b), (b) = temp)
25 sort_sv(size_t m
, size_t n
, size_t k
,
26 double **a
, double *w
, double **v
)
31 for (i
= k
; (i
< n
- 1) && (w
[i
] < w
[i
+1]); i
++) {
33 for (j
= 0; j
< m
; j
++) SWAPD(a
[j
][i
], a
[j
][i
+1]);
34 for (j
= 0; j
< n
; j
++) SWAPD(v
[j
][i
], v
[j
][i
+1]);
38 static double maxarg1
, maxarg2
;
39 #define Max(a, b) (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
40 #define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
43 svdcmp(double **a
, size_t m
, size_t n
, double *w
, double **v
)
46 long i
, k
, l
; /* These must be signed */
48 double c
, f
, h
, s
, x
, y
, z
;
49 double anorm
= 0.0, g
= 0.0, scale
= 0.0;
52 if (m
< n
) return(FALSE
); /* flag an error if m < n */
56 /* Householder reduction to bidiagonal form */
57 for (i
= 0; i
< n
; i
++) {
59 /* left-hand reduction */
64 for (k
= i
; k
< m
; k
++) scale
+= fabs(a
[k
][i
]);
66 for (k
= i
; k
< m
; k
++) {
68 s
+= a
[k
][i
] * a
[k
][i
];
71 g
= -SIGN(sqrt(s
), f
);
75 for (j
= l
; j
< n
; j
++) {
76 for (s
= 0.0, k
= i
; k
< m
; k
++) s
+= a
[k
][i
] * a
[k
][j
];
78 for (k
= i
; k
< m
; k
++) a
[k
][j
] += f
* a
[k
][i
];
81 for (k
= i
; k
< m
; k
++) a
[k
][i
] *= scale
;
86 /* right-hand reduction */
88 if (i
< m
&& i
!= n
- 1) {
89 for (k
= l
; k
< n
; k
++) scale
+= fabs(a
[i
][k
]);
91 for (k
= l
; k
< n
; k
++) {
93 s
+= a
[i
][k
] * a
[i
][k
];
96 g
= -SIGN(sqrt(s
), f
);
99 for (k
= l
; k
< n
; k
++) rv1
[k
] = a
[i
][k
] / h
;
101 for (j
= l
; j
< m
; j
++) {
102 for (s
= 0.0, k
= l
; k
< n
; k
++) s
+= a
[j
][k
] * a
[i
][k
];
103 for (k
= l
; k
< n
; k
++) a
[j
][k
] += s
* rv1
[k
];
106 for (k
= l
; k
< n
; k
++) a
[i
][k
] *= scale
;
109 anorm
= Max(anorm
, (fabs(w
[i
]) + fabs(rv1
[i
])));
112 /* accumulate the right-hand transformation */
113 for (i
= n
- 1; i
>= 0; i
--) {
116 for (j
= l
; j
< n
; j
++)
117 v
[j
][i
] = (a
[i
][j
] / a
[i
][l
]) / g
;
119 for (j
= l
; j
< n
; j
++) {
120 for (s
= 0.0, k
= l
; k
< n
; k
++) s
+= a
[i
][k
] * v
[k
][j
];
121 for (k
= l
; k
< n
; k
++) v
[k
][j
] += s
* v
[k
][i
];
124 for (j
= l
; j
< n
; j
++) v
[i
][j
] = v
[j
][i
] = 0.0;
131 /* accumulate the left-hand transformation */
132 for (i
= n
- 1; i
>= 0; i
--) {
136 for (j
= l
; j
< n
; j
++) a
[i
][j
] = 0.0;
140 for (j
= l
; j
< n
; j
++) {
141 for (s
= 0.0, k
= l
; k
< m
; k
++) s
+= a
[k
][i
] * a
[k
][j
];
142 f
= (s
/ a
[i
][i
]) * g
;
143 for (k
= i
; k
< m
; k
++) a
[k
][j
] += f
* a
[k
][i
];
146 for (j
= i
; j
< m
; j
++) a
[j
][i
] *= g
;
149 for (j
= i
; j
< m
; j
++) a
[j
][i
] = 0.0;
154 /* diagonalize the bidiagonal form */
155 for (k
= n
- 1; k
>= 0; k
--) { /* loop over singular values */
156 for (its
= 0; its
< 30; its
++) { /* loop over allowed iterations */
158 for (l
= k
; l
>= 0; l
--) { /* test for splitting */
160 if (fabs(rv1
[l
]) + anorm
== anorm
) {
164 if (fabs(w
[nm
]) + anorm
== anorm
) break;
169 for (i
= l
; i
<= k
; i
++) {
171 if (fabs(f
) + anorm
!= anorm
) {
177 sprintf(s
, "h = %f, f = %f, g = %f\n", h
, f
, g
);
178 fprintf(stdout
,"%s",s
);
183 for (j
= 0; j
< m
; j
++) {
186 a
[j
][nm
] = y
* c
+ z
* s
;
187 a
[j
][i
] = z
* c
- y
* s
;
193 if (l
== k
) { /* convergence */
194 if (z
< 0.0) { /* make singular value nonnegative */
196 for (j
= 0; j
< n
; j
++) v
[j
][k
] = (-v
[j
][k
]);
198 sort_sv(m
, n
, k
, a
, w
, v
);
203 return(FALSE
); /* return an error flag */
206 /* shift from bottom 2 x 2 minor */
212 f
= ((y
- z
) * (y
+ z
) + (g
- h
) * (g
+ h
)) / (2.0 * h
* y
);
214 f
= ((x
- z
) * (x
+ z
) + h
* ((y
/ (f
+ SIGN(g
, f
))) - h
)) / x
;
216 /* next QR transformation */
218 for (j
= l
; j
<= nm
; j
++) {
232 for (jj
= 0; jj
< n
; jj
++) {
235 v
[jj
][j
] = x
* c
+ z
* s
;
236 v
[jj
][i
] = z
* c
- x
* s
;
245 f
= (c
* g
) + (s
* y
);
246 x
= (c
* y
) - (s
* g
);
247 for (jj
= 0; jj
< m
; jj
++) {
250 a
[jj
][j
] = y
* c
+ z
* s
;
251 a
[jj
][i
] = z
* c
- y
* s
;