3 * Copyright (C) 2007 Brian Gough
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
21 chop_small_elements (gsl_vector * d, gsl_vector * f)
23 const size_t N = d->size;
24 double d_i = gsl_vector_get (d, 0);
28 for (i = 0; i < N - 1; i++)
30 double f_i = gsl_vector_get (f, i);
31 double d_ip1 = gsl_vector_get (d, i + 1);
33 if (fabs (f_i) < GSL_DBL_EPSILON * (fabs (d_i) + fabs (d_ip1)))
35 gsl_vector_set (f, i, 0.0);
43 trailing_eigenvalue (const gsl_vector * d, const gsl_vector * f)
45 const size_t n = d->size;
47 double da = gsl_vector_get (d, n - 2);
48 double db = gsl_vector_get (d, n - 1);
49 double fa = (n > 2) ? gsl_vector_get (f, n - 3) : 0.0;
50 double fb = gsl_vector_get (f, n - 2);
52 double ta = da * da + fa * fa;
53 double tb = db * db + fb * fb;
56 double dt = (ta - tb) / 2.0;
62 mu = tb - (tab * tab) / (dt + hypot (dt, tab));
66 mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab));
73 create_schur (double d0, double f0, double d1, double * c, double * s)
75 double apq = 2.0 * d0 * f0;
77 if (d0 == 0 || f0 == 0)
84 /* Check if we need to rescale to avoid underflow/overflow */
85 if (fabs(d0) < GSL_SQRT_DBL_MIN || fabs(d0) > GSL_SQRT_DBL_MAX
86 || fabs(f0) < GSL_SQRT_DBL_MIN || fabs(f0) > GSL_SQRT_DBL_MAX
87 || fabs(d1) < GSL_SQRT_DBL_MIN || fabs(d1) > GSL_SQRT_DBL_MAX)
93 /* Bring |d0*f0| into the range GSL_DBL_MIN to GSL_DBL_MAX */
94 scale = ldexp(1.0, -(d0_exp + f0_exp)/4);
104 double tau = (f0*f0 + (d1 + d0)*(d1 - d0)) / apq;
108 t = 1.0/(tau + hypot(1.0, tau));
112 t = -1.0/(-tau + hypot(1.0, tau));
115 *c = 1.0 / hypot(1.0, t);
126 svd2 (gsl_vector * d, gsl_vector * f, gsl_matrix * U, gsl_matrix * V)
129 double c, s, a11, a12, a21, a22;
131 const size_t M = U->size1;
132 const size_t N = V->size1;
134 double d0 = gsl_vector_get (d, 0);
135 double f0 = gsl_vector_get (f, 0);
137 double d1 = gsl_vector_get (d, 1);
141 /* Eliminate off-diagonal element in [0,f0;0,d1] to make [d,0;0,0] */
143 create_givens (f0, d1, &c, &s);
145 /* compute B <= G^T B X, where X = [0,1;1,0] */
147 gsl_vector_set (d, 0, c * f0 - s * d1);
148 gsl_vector_set (f, 0, s * f0 + c * d1);
149 gsl_vector_set (d, 1, 0.0);
151 /* Compute U <= U G */
153 for (i = 0; i < M; i++)
155 double Uip = gsl_matrix_get (U, i, 0);
156 double Uiq = gsl_matrix_get (U, i, 1);
157 gsl_matrix_set (U, i, 0, c * Uip - s * Uiq);
158 gsl_matrix_set (U, i, 1, s * Uip + c * Uiq);
161 /* Compute V <= V X */
163 gsl_matrix_swap_columns (V, 0, 1);
169 /* Eliminate off-diagonal element in [d0,f0;0,0] */
171 create_givens (d0, f0, &c, &s);
173 /* compute B <= B G */
175 gsl_vector_set (d, 0, d0 * c - f0 * s);
176 gsl_vector_set (f, 0, 0.0);
178 /* Compute V <= V G */
180 for (i = 0; i < N; i++)
182 double Vip = gsl_matrix_get (V, i, 0);
183 double Viq = gsl_matrix_get (V, i, 1);
184 gsl_matrix_set (V, i, 0, c * Vip - s * Viq);
185 gsl_matrix_set (V, i, 1, s * Vip + c * Viq);
192 /* Make columns orthogonal, A = [d0, f0; 0, d1] * G */
194 create_schur (d0, f0, d1, &c, &s);
196 /* compute B <= B G */
198 a11 = c * d0 - s * f0;
201 a12 = s * d0 + c * f0;
204 /* Compute V <= V G */
206 for (i = 0; i < N; i++)
208 double Vip = gsl_matrix_get (V, i, 0);
209 double Viq = gsl_matrix_get (V, i, 1);
210 gsl_matrix_set (V, i, 0, c * Vip - s * Viq);
211 gsl_matrix_set (V, i, 1, s * Vip + c * Viq);
214 /* Eliminate off-diagonal elements, bring column with largest
215 norm to first column */
217 if (hypot(a11, a21) < hypot(a12,a22))
223 t1 = a11; a11 = a12; a12 = t1;
224 t2 = a21; a21 = a22; a22 = t2;
228 gsl_matrix_swap_columns(V, 0, 1);
231 create_givens (a11, a21, &c, &s);
233 /* compute B <= G^T B */
235 gsl_vector_set (d, 0, c * a11 - s * a21);
236 gsl_vector_set (f, 0, c * a12 - s * a22);
237 gsl_vector_set (d, 1, s * a12 + c * a22);
239 /* Compute U <= U G */
241 for (i = 0; i < M; i++)
243 double Uip = gsl_matrix_get (U, i, 0);
244 double Uiq = gsl_matrix_get (U, i, 1);
245 gsl_matrix_set (U, i, 0, c * Uip - s * Uiq);
246 gsl_matrix_set (U, i, 1, s * Uip + c * Uiq);
255 chase_out_intermediate_zero (gsl_vector * d, gsl_vector * f, gsl_matrix * U, size_t k0)
258 const size_t M = U->size1;
260 const size_t n = d->size;
265 x = gsl_vector_get (f, k0);
266 y = gsl_vector_get (d, k0+1);
268 for (k = k0; k < n - 1; k++)
270 create_givens (y, -x, &c, &s);
272 /* Compute U <= U G */
276 gsl_vector_view Uk0 = gsl_matrix_column(U,k0);
277 gsl_vector_view Ukp1 = gsl_matrix_column(U,k+1);
278 gsl_blas_drot(&Uk0.vector, &Ukp1.vector, c, -s);
284 for (i = 0; i < M; i++)
286 double Uip = gsl_matrix_get (U, i, k0);
287 double Uiq = gsl_matrix_get (U, i, k + 1);
288 gsl_matrix_set (U, i, k0, c * Uip - s * Uiq);
289 gsl_matrix_set (U, i, k + 1, s * Uip + c * Uiq);
294 /* compute B <= G^T B */
296 gsl_vector_set (d, k + 1, s * x + c * y);
299 gsl_vector_set (f, k, c * x - s * y );
303 double z = gsl_vector_get (f, k + 1);
304 gsl_vector_set (f, k + 1, c * z);
307 y = gsl_vector_get (d, k + 2);
313 chase_out_trailing_zero (gsl_vector * d, gsl_vector * f, gsl_matrix * V)
316 const size_t N = V->size1;
318 const size_t n = d->size;
323 x = gsl_vector_get (d, n - 2);
324 y = gsl_vector_get (f, n - 2);
326 for (k = n - 1; k > 0 && k--;)
328 create_givens (x, y, &c, &s);
330 /* Compute V <= V G where G = [c, s ; -s, c] */
334 gsl_vector_view Vp = gsl_matrix_column(V,k);
335 gsl_vector_view Vq = gsl_matrix_column(V,n-1);
336 gsl_blas_drot(&Vp.vector, &Vq.vector, c, -s);
342 for (i = 0; i < N; i++)
344 double Vip = gsl_matrix_get (V, i, k);
345 double Viq = gsl_matrix_get (V, i, n - 1);
346 gsl_matrix_set (V, i, k, c * Vip - s * Viq);
347 gsl_matrix_set (V, i, n - 1, s * Vip + c * Viq);
352 /* compute B <= B G */
354 gsl_vector_set (d, k, c * x - s * y);
357 gsl_vector_set (f, k, s * x + c * y );
361 double z = gsl_vector_get (f, k - 1);
362 gsl_vector_set (f, k - 1, c * z);
364 x = gsl_vector_get (d, k - 1);
371 qrstep (gsl_vector * d, gsl_vector * f, gsl_matrix * U, gsl_matrix * V)
374 const size_t M = U->size1;
375 const size_t N = V->size1;
377 const size_t n = d->size;
379 double ak, bk, zk, ap, bp, aq, bq;
383 return; /* shouldn't happen */
385 /* Compute 2x2 svd directly */
393 /* Chase out any zeroes on the diagonal */
395 for (i = 0; i < n - 1; i++)
397 double d_i = gsl_vector_get (d, i);
401 chase_out_intermediate_zero (d, f, U, i);
406 /* Chase out any zero at the end of the diagonal */
409 double d_nm1 = gsl_vector_get (d, n - 1);
413 chase_out_trailing_zero (d, f, V);
419 /* Apply QR reduction steps to the diagonal and offdiagonal */
422 double d0 = gsl_vector_get (d, 0);
423 double f0 = gsl_vector_get (f, 0);
425 double d1 = gsl_vector_get (d, 1);
426 double f1 = gsl_vector_get (f, 1);
429 double mu = trailing_eigenvalue (d, f);
435 /* Set up the recurrence for Givens rotations on a bidiagonal matrix */
447 for (k = 0; k < n - 1; k++)
450 create_givens (y, z, &c, &s);
452 /* Compute V <= V G */
456 gsl_vector_view Vk = gsl_matrix_column(V,k);
457 gsl_vector_view Vkp1 = gsl_matrix_column(V,k+1);
458 gsl_blas_drot(&Vk.vector, &Vkp1.vector, c, -s);
461 for (i = 0; i < N; i++)
463 double Vip = gsl_matrix_get (V, i, k);
464 double Viq = gsl_matrix_get (V, i, k + 1);
465 gsl_matrix_set (V, i, k, c * Vip - s * Viq);
466 gsl_matrix_set (V, i, k + 1, s * Vip + c * Viq);
470 /* compute B <= B G */
473 double bk1 = c * bk - s * z;
475 double ap1 = c * ap - s * bp;
476 double bp1 = s * ap + c * bp;
477 double zp1 = -s * aq;
483 gsl_vector_set (f, k - 1, bk1);
494 bp = gsl_vector_get (f, k + 1);
505 create_givens (y, z, &c, &s);
507 /* Compute U <= U G */
511 gsl_vector_view Uk = gsl_matrix_column(U,k);
512 gsl_vector_view Ukp1 = gsl_matrix_column(U,k+1);
513 gsl_blas_drot(&Uk.vector, &Ukp1.vector, c, -s);
516 for (i = 0; i < M; i++)
518 double Uip = gsl_matrix_get (U, i, k);
519 double Uiq = gsl_matrix_get (U, i, k + 1);
520 gsl_matrix_set (U, i, k, c * Uip - s * Uiq);
521 gsl_matrix_set (U, i, k + 1, s * Uip + c * Uiq);
525 /* compute B <= G^T B */
528 double ak1 = c * ak - s * zk;
529 double bk1 = c * bk - s * ap;
530 double zk1 = -s * bp;
532 double ap1 = s * bk + c * ap;
535 gsl_vector_set (d, k, ak1);
546 aq = gsl_vector_get (d, k + 2);
558 gsl_vector_set (f, n - 2, bk);
559 gsl_vector_set (d, n - 1, ap);