Mercurial > octave-nkf

comparison libcruft/lapack/dbdsqr.f @ 7034:68db500cb558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-16 18:54:19 by jwe]
author jwe
date Tue, 16 Oct 2007 18:54:23 +0000
parents edcaebe1b81b
children
comparison
equal deleted inserted replaced
7033:f0142f2afdc6 7034:68db500cb558
1 SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, 1 SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
2 $ LDU, C, LDC, WORK, INFO ) 2 $ LDU, C, LDC, WORK, INFO )
3 * 3 *
4 * -- LAPACK routine (version 3.0) -- 4 * -- LAPACK routine (version 3.1.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * Courant Institute, Argonne National Lab, and Rice University 6 * January 2007
7 * October 31, 1999
8 * 7 *
9 * .. Scalar Arguments .. 8 * .. Scalar Arguments ..
10 CHARACTER UPLO 9 CHARACTER UPLO
11 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU 10 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
12 * .. 11 * ..
16 * .. 15 * ..
17 * 16 *
18 * Purpose 17 * Purpose
19 * ======= 18 * =======
20 * 19 *
21 * DBDSQR computes the singular value decomposition (SVD) of a real 20 * DBDSQR computes the singular values and, optionally, the right and/or
22 * N-by-N (upper or lower) bidiagonal matrix B: B = Q * S * P' (P' 21 * left singular vectors from the singular value decomposition (SVD) of
23 * denotes the transpose of P), where S is a diagonal matrix with 22 * a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
24 * non-negative diagonal elements (the singular values of B), and Q 23 * zero-shift QR algorithm. The SVD of B has the form
25 * and P are orthogonal matrices. 24 *
26 * 25 * B = Q * S * P**T
27 * The routine computes S, and optionally computes U * Q, P' * VT, 26 *
28 * or Q' * C, for given real input matrices U, VT, and C. 27 * where S is the diagonal matrix of singular values, Q is an orthogonal
28 * matrix of left singular vectors, and P is an orthogonal matrix of
29 * right singular vectors. If left singular vectors are requested, this
30 * subroutine actually returns U*Q instead of Q, and, if right singular
31 * vectors are requested, this subroutine returns P**T*VT instead of
32 * P**T, for given real input matrices U and VT. When U and VT are the
33 * orthogonal matrices that reduce a general matrix A to bidiagonal
34 * form: A = U*B*VT, as computed by DGEBRD, then
35 *
36 * A = (U*Q) * S * (P**T*VT)
37 *
38 * is the SVD of A. Optionally, the subroutine may also compute Q**T*C
39 * for a given real input matrix C.
29 * 40 *
30 * See "Computing Small Singular Values of Bidiagonal Matrices With 41 * See "Computing Small Singular Values of Bidiagonal Matrices With
31 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, 42 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
32 * LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, 43 * LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
33 * no. 5, pp. 873-912, Sept 1990) and 44 * no. 5, pp. 873-912, Sept 1990) and
58 * D (input/output) DOUBLE PRECISION array, dimension (N) 69 * D (input/output) DOUBLE PRECISION array, dimension (N)
59 * On entry, the n diagonal elements of the bidiagonal matrix B. 70 * On entry, the n diagonal elements of the bidiagonal matrix B.
60 * On exit, if INFO=0, the singular values of B in decreasing 71 * On exit, if INFO=0, the singular values of B in decreasing
61 * order. 72 * order.
62 * 73 *
63 * E (input/output) DOUBLE PRECISION array, dimension (N) 74 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
64 * On entry, the elements of E contain the 75 * On entry, the N-1 offdiagonal elements of the bidiagonal
65 * offdiagonal elements of the bidiagonal matrix whose SVD 76 * matrix B.
66 * is desired. On normal exit (INFO = 0), E is destroyed. 77 * On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
67 * If the algorithm does not converge (INFO > 0), D and E
68 * will contain the diagonal and superdiagonal elements of a 78 * will contain the diagonal and superdiagonal elements of a
69 * bidiagonal matrix orthogonally equivalent to the one given 79 * bidiagonal matrix orthogonally equivalent to the one given
70 * as input. E(N) is used for workspace. 80 * as input.
71 * 81 *
72 * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) 82 * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
73 * On entry, an N-by-NCVT matrix VT. 83 * On entry, an N-by-NCVT matrix VT.
74 * On exit, VT is overwritten by P' * VT. 84 * On exit, VT is overwritten by P**T * VT.
75 * VT is not referenced if NCVT = 0. 85 * Not referenced if NCVT = 0.
76 * 86 *
77 * LDVT (input) INTEGER 87 * LDVT (input) INTEGER
78 * The leading dimension of the array VT. 88 * The leading dimension of the array VT.
79 * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. 89 * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
80 * 90 *
81 * U (input/output) DOUBLE PRECISION array, dimension (LDU, N) 91 * U (input/output) DOUBLE PRECISION array, dimension (LDU, N)
82 * On entry, an NRU-by-N matrix U. 92 * On entry, an NRU-by-N matrix U.
83 * On exit, U is overwritten by U * Q. 93 * On exit, U is overwritten by U * Q.
84 * U is not referenced if NRU = 0. 94 * Not referenced if NRU = 0.
85 * 95 *
86 * LDU (input) INTEGER 96 * LDU (input) INTEGER
87 * The leading dimension of the array U. LDU >= max(1,NRU). 97 * The leading dimension of the array U. LDU >= max(1,NRU).
88 * 98 *
89 * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) 99 * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
90 * On entry, an N-by-NCC matrix C. 100 * On entry, an N-by-NCC matrix C.
91 * On exit, C is overwritten by Q' * C. 101 * On exit, C is overwritten by Q**T * C.
92 * C is not referenced if NCC = 0. 102 * Not referenced if NCC = 0.
93 * 103 *
94 * LDC (input) INTEGER 104 * LDC (input) INTEGER
95 * The leading dimension of the array C. 105 * The leading dimension of the array C.
96 * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. 106 * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
97 * 107 *
98 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 108 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
109 * if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise
99 * 110 *
100 * INFO (output) INTEGER 111 * INFO (output) INTEGER
101 * = 0: successful exit 112 * = 0: successful exit
102 * < 0: If INFO = -i, the i-th argument had an illegal value 113 * < 0: If INFO = -i, the i-th argument had an illegal value
103 * > 0: the algorithm did not converge; D and E contain the 114 * > 0: the algorithm did not converge; D and E contain the
153 LOGICAL LOWER, ROTATE 164 LOGICAL LOWER, ROTATE
154 INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1, 165 INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
155 $ NM12, NM13, OLDLL, OLDM 166 $ NM12, NM13, OLDLL, OLDM
156 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, 167 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
157 $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, 168 $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
158 $ SINR, SLL, SMAX, SMIN, SMINL, SMINLO, SMINOA, 169 $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
159 $ SN, THRESH, TOL, TOLMUL, UNFL 170 $ SN, THRESH, TOL, TOLMUL, UNFL
160 * .. 171 * ..
161 * .. External Functions .. 172 * .. External Functions ..
162 LOGICAL LSAME 173 LOGICAL LSAME
163 DOUBLE PRECISION DLAMCH 174 DOUBLE PRECISION DLAMCH
413 DO 100 LLL = LL, M - 1 424 DO 100 LLL = LL, M - 1
414 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN 425 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
415 E( LLL ) = ZERO 426 E( LLL ) = ZERO
416 GO TO 60 427 GO TO 60
417 END IF 428 END IF
418 SMINLO = SMINL
419 MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) 429 MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
420 SMINL = MIN( SMINL, MU ) 430 SMINL = MIN( SMINL, MU )
421 100 CONTINUE 431 100 CONTINUE
422 END IF 432 END IF
423 * 433 *
442 DO 110 LLL = M - 1, LL, -1 452 DO 110 LLL = M - 1, LL, -1
443 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN 453 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
444 E( LLL ) = ZERO 454 E( LLL ) = ZERO
445 GO TO 60 455 GO TO 60
446 END IF 456 END IF
447 SMINLO = SMINL
448 MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) 457 MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
449 SMINL = MIN( SMINL, MU ) 458 SMINL = MIN( SMINL, MU )
450 110 CONTINUE 459 110 CONTINUE
451 END IF 460 END IF
452 END IF 461 END IF