 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, 00002 $ LDU, C, LDC, WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * January 2007 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER UPLO 00011 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), 00015 $ VT( LDVT, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DBDSQR computes the singular values and, optionally, the right and/or 00022 * left singular vectors from the singular value decomposition (SVD) of 00023 * a real N-by-N (upper or lower) bidiagonal matrix B using the implicit 00024 * zero-shift QR algorithm. The SVD of B has the form 00025 * 00026 * B = Q * S * P**T 00027 * 00028 * where S is the diagonal matrix of singular values, Q is an orthogonal 00029 * matrix of left singular vectors, and P is an orthogonal matrix of 00030 * right singular vectors. If left singular vectors are requested, this 00031 * subroutine actually returns U*Q instead of Q, and, if right singular 00032 * vectors are requested, this subroutine returns P**T*VT instead of 00033 * P**T, for given real input matrices U and VT. When U and VT are the 00034 * orthogonal matrices that reduce a general matrix A to bidiagonal 00035 * form: A = U*B*VT, as computed by DGEBRD, then 00036 * 00037 * A = (U*Q) * S * (P**T*VT) 00038 * 00039 * is the SVD of A. Optionally, the subroutine may also compute Q**T*C 00040 * for a given real input matrix C. 00041 * 00042 * See "Computing Small Singular Values of Bidiagonal Matrices With 00043 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, 00044 * LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, 00045 * no. 5, pp. 873-912, Sept 1990) and 00046 * "Accurate singular values and differential qd algorithms," by 00047 * B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics 00048 * Department, University of California at Berkeley, July 1992 00049 * for a detailed description of the algorithm. 00050 * 00051 * Arguments 00052 * ========= 00053 * 00054 * UPLO (input) CHARACTER*1 00055 * = 'U': B is upper bidiagonal; 00056 * = 'L': B is lower bidiagonal. 00057 * 00058 * N (input) INTEGER 00059 * The order of the matrix B. N >= 0. 00060 * 00061 * NCVT (input) INTEGER 00062 * The number of columns of the matrix VT. NCVT >= 0. 00063 * 00064 * NRU (input) INTEGER 00065 * The number of rows of the matrix U. NRU >= 0. 00066 * 00067 * NCC (input) INTEGER 00068 * The number of columns of the matrix C. NCC >= 0. 00069 * 00070 * D (input/output) DOUBLE PRECISION array, dimension (N) 00071 * On entry, the n diagonal elements of the bidiagonal matrix B. 00072 * On exit, if INFO=0, the singular values of B in decreasing 00073 * order. 00074 * 00075 * E (input/output) DOUBLE PRECISION array, dimension (N-1) 00076 * On entry, the N-1 offdiagonal elements of the bidiagonal 00077 * matrix B. 00078 * On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E 00079 * will contain the diagonal and superdiagonal elements of a 00080 * bidiagonal matrix orthogonally equivalent to the one given 00081 * as input. 00082 * 00083 * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) 00084 * On entry, an N-by-NCVT matrix VT. 00085 * On exit, VT is overwritten by P**T * VT. 00086 * Not referenced if NCVT = 0. 00087 * 00088 * LDVT (input) INTEGER 00089 * The leading dimension of the array VT. 00090 * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. 00091 * 00092 * U (input/output) DOUBLE PRECISION array, dimension (LDU, N) 00093 * On entry, an NRU-by-N matrix U. 00094 * On exit, U is overwritten by U * Q. 00095 * Not referenced if NRU = 0. 00096 * 00097 * LDU (input) INTEGER 00098 * The leading dimension of the array U. LDU >= max(1,NRU). 00099 * 00100 * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) 00101 * On entry, an N-by-NCC matrix C. 00102 * On exit, C is overwritten by Q**T * C. 00103 * Not referenced if NCC = 0. 00104 * 00105 * LDC (input) INTEGER 00106 * The leading dimension of the array C. 00107 * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. 00108 * 00109 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 00110 * 00111 * INFO (output) INTEGER 00112 * = 0: successful exit 00113 * < 0: If INFO = -i, the i-th argument had an illegal value 00114 * > 0: 00115 * if NCVT = NRU = NCC = 0, 00116 * = 1, a split was marked by a positive value in E 00117 * = 2, current block of Z not diagonalized after 30*N 00118 * iterations (in inner while loop) 00119 * = 3, termination criterion of outer while loop not met 00120 * (program created more than N unreduced blocks) 00121 * else NCVT = NRU = NCC = 0, 00122 * the algorithm did not converge; D and E contain the 00123 * elements of a bidiagonal matrix which is orthogonally 00124 * similar to the input matrix B; if INFO = i, i 00125 * elements of E have not converged to zero. 00126 * 00127 * Internal Parameters 00128 * =================== 00129 * 00130 * TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8))) 00131 * TOLMUL controls the convergence criterion of the QR loop. 00132 * If it is positive, TOLMUL*EPS is the desired relative 00133 * precision in the computed singular values. 00134 * If it is negative, abs(TOLMUL*EPS*sigma_max) is the 00135 * desired absolute accuracy in the computed singular 00136 * values (corresponds to relative accuracy 00137 * abs(TOLMUL*EPS) in the largest singular value. 00138 * abs(TOLMUL) should be between 1 and 1/EPS, and preferably 00139 * between 10 (for fast convergence) and .1/EPS 00140 * (for there to be some accuracy in the results). 00141 * Default is to lose at either one eighth or 2 of the 00142 * available decimal digits in each computed singular value 00143 * (whichever is smaller). 00144 * 00145 * MAXITR INTEGER, default = 6 00146 * MAXITR controls the maximum number of passes of the 00147 * algorithm through its inner loop. The algorithms stops 00148 * (and so fails to converge) if the number of passes 00149 * through the inner loop exceeds MAXITR*N**2. 00150 * 00151 * ===================================================================== 00152 * 00153 * .. Parameters .. 00154 DOUBLE PRECISION ZERO 00155 PARAMETER ( ZERO = 0.0D0 ) 00156 DOUBLE PRECISION ONE 00157 PARAMETER ( ONE = 1.0D0 ) 00158 DOUBLE PRECISION NEGONE 00159 PARAMETER ( NEGONE = -1.0D0 ) 00160 DOUBLE PRECISION HNDRTH 00161 PARAMETER ( HNDRTH = 0.01D0 ) 00162 DOUBLE PRECISION TEN 00163 PARAMETER ( TEN = 10.0D0 ) 00164 DOUBLE PRECISION HNDRD 00165 PARAMETER ( HNDRD = 100.0D0 ) 00166 DOUBLE PRECISION MEIGTH 00167 PARAMETER ( MEIGTH = -0.125D0 ) 00168 INTEGER MAXITR 00169 PARAMETER ( MAXITR = 6 ) 00170 * .. 00171 * .. Local Scalars .. 00172 LOGICAL LOWER, ROTATE 00173 INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1, 00174 $ NM12, NM13, OLDLL, OLDM 00175 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, 00176 $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, 00177 $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA, 00178 $ SN, THRESH, TOL, TOLMUL, UNFL 00179 * .. 00180 * .. External Functions .. 00181 LOGICAL LSAME 00182 DOUBLE PRECISION DLAMCH 00183 EXTERNAL LSAME, DLAMCH 00184 * .. 00185 * .. External Subroutines .. 00186 EXTERNAL DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT, 00187 $ DSCAL, DSWAP, XERBLA 00188 * .. 00189 * .. Intrinsic Functions .. 00190 INTRINSIC ABS, DBLE, MAX, MIN, SIGN, SQRT 00191 * .. 00192 * .. Executable Statements .. 00193 * 00194 * Test the input parameters. 00195 * 00196 INFO = 0 00197 LOWER = LSAME( UPLO, 'L' ) 00198 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN 00199 INFO = -1 00200 ELSE IF( N.LT.0 ) THEN 00201 INFO = -2 00202 ELSE IF( NCVT.LT.0 ) THEN 00203 INFO = -3 00204 ELSE IF( NRU.LT.0 ) THEN 00205 INFO = -4 00206 ELSE IF( NCC.LT.0 ) THEN 00207 INFO = -5 00208 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. 00209 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN 00210 INFO = -9 00211 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN 00212 INFO = -11 00213 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. 00214 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN 00215 INFO = -13 00216 END IF 00217 IF( INFO.NE.0 ) THEN 00218 CALL XERBLA( 'DBDSQR', -INFO ) 00219 RETURN 00220 END IF 00221 IF( N.EQ.0 ) 00222 $ RETURN 00223 IF( N.EQ.1 ) 00224 $ GO TO 160 00225 * 00226 * ROTATE is true if any singular vectors desired, false otherwise 00227 * 00228 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) 00229 * 00230 * If no singular vectors desired, use qd algorithm 00231 * 00232 IF( .NOT.ROTATE ) THEN 00233 CALL DLASQ1( N, D, E, WORK, INFO ) 00234 RETURN 00235 END IF 00236 * 00237 NM1 = N - 1 00238 NM12 = NM1 + NM1 00239 NM13 = NM12 + NM1 00240 IDIR = 0 00241 * 00242 * Get machine constants 00243 * 00244 EPS = DLAMCH( 'Epsilon' ) 00245 UNFL = DLAMCH( 'Safe minimum' ) 00246 * 00247 * If matrix lower bidiagonal, rotate to be upper bidiagonal 00248 * by applying Givens rotations on the left 00249 * 00250 IF( LOWER ) THEN 00251 DO 10 I = 1, N - 1 00252 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 00253 D( I ) = R 00254 E( I ) = SN*D( I+1 ) 00255 D( I+1 ) = CS*D( I+1 ) 00256 WORK( I ) = CS 00257 WORK( NM1+I ) = SN 00258 10 CONTINUE 00259 * 00260 * Update singular vectors if desired 00261 * 00262 IF( NRU.GT.0 ) 00263 $ CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U, 00264 $ LDU ) 00265 IF( NCC.GT.0 ) 00266 $ CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C, 00267 $ LDC ) 00268 END IF 00269 * 00270 * Compute singular values to relative accuracy TOL 00271 * (By setting TOL to be negative, algorithm will compute 00272 * singular values to absolute accuracy ABS(TOL)*norm(input matrix)) 00273 * 00274 TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) ) 00275 TOL = TOLMUL*EPS 00276 * 00277 * Compute approximate maximum, minimum singular values 00278 * 00279 SMAX = ZERO 00280 DO 20 I = 1, N 00281 SMAX = MAX( SMAX, ABS( D( I ) ) ) 00282 20 CONTINUE 00283 DO 30 I = 1, N - 1 00284 SMAX = MAX( SMAX, ABS( E( I ) ) ) 00285 30 CONTINUE 00286 SMINL = ZERO 00287 IF( TOL.GE.ZERO ) THEN 00288 * 00289 * Relative accuracy desired 00290 * 00291 SMINOA = ABS( D( 1 ) ) 00292 IF( SMINOA.EQ.ZERO ) 00293 $ GO TO 50 00294 MU = SMINOA 00295 DO 40 I = 2, N 00296 MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) ) 00297 SMINOA = MIN( SMINOA, MU ) 00298 IF( SMINOA.EQ.ZERO ) 00299 $ GO TO 50 00300 40 CONTINUE 00301 50 CONTINUE 00302 SMINOA = SMINOA / SQRT( DBLE( N ) ) 00303 THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL ) 00304 ELSE 00305 * 00306 * Absolute accuracy desired 00307 * 00308 THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL ) 00309 END IF 00310 * 00311 * Prepare for main iteration loop for the singular values 00312 * (MAXIT is the maximum number of passes through the inner 00313 * loop permitted before nonconvergence signalled.) 00314 * 00315 MAXIT = MAXITR*N*N 00316 ITER = 0 00317 OLDLL = -1 00318 OLDM = -1 00319 * 00320 * M points to last element of unconverged part of matrix 00321 * 00322 M = N 00323 * 00324 * Begin main iteration loop 00325 * 00326 60 CONTINUE 00327 * 00328 * Check for convergence or exceeding iteration count 00329 * 00330 IF( M.LE.1 ) 00331 $ GO TO 160 00332 IF( ITER.GT.MAXIT ) 00333 $ GO TO 200 00334 * 00335 * Find diagonal block of matrix to work on 00336 * 00337 IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH ) 00338 $ D( M ) = ZERO 00339 SMAX = ABS( D( M ) ) 00340 SMIN = SMAX 00341 DO 70 LLL = 1, M - 1 00342 LL = M - LLL 00343 ABSS = ABS( D( LL ) ) 00344 ABSE = ABS( E( LL ) ) 00345 IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH ) 00346 $ D( LL ) = ZERO 00347 IF( ABSE.LE.THRESH ) 00348 $ GO TO 80 00349 SMIN = MIN( SMIN, ABSS ) 00350 SMAX = MAX( SMAX, ABSS, ABSE ) 00351 70 CONTINUE 00352 LL = 0 00353 GO TO 90 00354 80 CONTINUE 00355 E( LL ) = ZERO 00356 * 00357 * Matrix splits since E(LL) = 0 00358 * 00359 IF( LL.EQ.M-1 ) THEN 00360 * 00361 * Convergence of bottom singular value, return to top of loop 00362 * 00363 M = M - 1 00364 GO TO 60 00365 END IF 00366 90 CONTINUE 00367 LL = LL + 1 00368 * 00369 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero 00370 * 00371 IF( LL.EQ.M-1 ) THEN 00372 * 00373 * 2 by 2 block, handle separately 00374 * 00375 CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR, 00376 $ COSR, SINL, COSL ) 00377 D( M-1 ) = SIGMX 00378 E( M-1 ) = ZERO 00379 D( M ) = SIGMN 00380 * 00381 * Compute singular vectors, if desired 00382 * 00383 IF( NCVT.GT.0 ) 00384 $ CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR, 00385 $ SINR ) 00386 IF( NRU.GT.0 ) 00387 $ CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL ) 00388 IF( NCC.GT.0 ) 00389 $ CALL DROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL, 00390 $ SINL ) 00391 M = M - 2 00392 GO TO 60 00393 END IF 00394 * 00395 * If working on new submatrix, choose shift direction 00396 * (from larger end diagonal element towards smaller) 00397 * 00398 IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN 00399 IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN 00400 * 00401 * Chase bulge from top (big end) to bottom (small end) 00402 * 00403 IDIR = 1 00404 ELSE 00405 * 00406 * Chase bulge from bottom (big end) to top (small end) 00407 * 00408 IDIR = 2 00409 END IF 00410 END IF 00411 * 00412 * Apply convergence tests 00413 * 00414 IF( IDIR.EQ.1 ) THEN 00415 * 00416 * Run convergence test in forward direction 00417 * First apply standard test to bottom of matrix 00418 * 00419 IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR. 00420 $ ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN 00421 E( M-1 ) = ZERO 00422 GO TO 60 00423 END IF 00424 * 00425 IF( TOL.GE.ZERO ) THEN 00426 * 00427 * If relative accuracy desired, 00428 * apply convergence criterion forward 00429 * 00430 MU = ABS( D( LL ) ) 00431 SMINL = MU 00432 DO 100 LLL = LL, M - 1 00433 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN 00434 E( LLL ) = ZERO 00435 GO TO 60 00436 END IF 00437 MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) 00438 SMINL = MIN( SMINL, MU ) 00439 100 CONTINUE 00440 END IF 00441 * 00442 ELSE 00443 * 00444 * Run convergence test in backward direction 00445 * First apply standard test to top of matrix 00446 * 00447 IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR. 00448 $ ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN 00449 E( LL ) = ZERO 00450 GO TO 60 00451 END IF 00452 * 00453 IF( TOL.GE.ZERO ) THEN 00454 * 00455 * If relative accuracy desired, 00456 * apply convergence criterion backward 00457 * 00458 MU = ABS( D( M ) ) 00459 SMINL = MU 00460 DO 110 LLL = M - 1, LL, -1 00461 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN 00462 E( LLL ) = ZERO 00463 GO TO 60 00464 END IF 00465 MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) 00466 SMINL = MIN( SMINL, MU ) 00467 110 CONTINUE 00468 END IF 00469 END IF 00470 OLDLL = LL 00471 OLDM = M 00472 * 00473 * Compute shift. First, test if shifting would ruin relative 00474 * accuracy, and if so set the shift to zero. 00475 * 00476 IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE. 00477 $ MAX( EPS, HNDRTH*TOL ) ) THEN 00478 * 00479 * Use a zero shift to avoid loss of relative accuracy 00480 * 00481 SHIFT = ZERO 00482 ELSE 00483 * 00484 * Compute the shift from 2-by-2 block at end of matrix 00485 * 00486 IF( IDIR.EQ.1 ) THEN 00487 SLL = ABS( D( LL ) ) 00488 CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R ) 00489 ELSE 00490 SLL = ABS( D( M ) ) 00491 CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R ) 00492 END IF 00493 * 00494 * Test if shift negligible, and if so set to zero 00495 * 00496 IF( SLL.GT.ZERO ) THEN 00497 IF( ( SHIFT / SLL )**2.LT.EPS ) 00498 $ SHIFT = ZERO 00499 END IF 00500 END IF 00501 * 00502 * Increment iteration count 00503 * 00504 ITER = ITER + M - LL 00505 * 00506 * If SHIFT = 0, do simplified QR iteration 00507 * 00508 IF( SHIFT.EQ.ZERO ) THEN 00509 IF( IDIR.EQ.1 ) THEN 00510 * 00511 * Chase bulge from top to bottom 00512 * Save cosines and sines for later singular vector updates 00513 * 00514 CS = ONE 00515 OLDCS = ONE 00516 DO 120 I = LL, M - 1 00517 CALL DLARTG( D( I )*CS, E( I ), CS, SN, R ) 00518 IF( I.GT.LL ) 00519 $ E( I-1 ) = OLDSN*R 00520 CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) ) 00521 WORK( I-LL+1 ) = CS 00522 WORK( I-LL+1+NM1 ) = SN 00523 WORK( I-LL+1+NM12 ) = OLDCS 00524 WORK( I-LL+1+NM13 ) = OLDSN 00525 120 CONTINUE 00526 H = D( M )*CS 00527 D( M ) = H*OLDCS 00528 E( M-1 ) = H*OLDSN 00529 * 00530 * Update singular vectors 00531 * 00532 IF( NCVT.GT.0 ) 00533 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), 00534 $ WORK( N ), VT( LL, 1 ), LDVT ) 00535 IF( NRU.GT.0 ) 00536 $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), 00537 $ WORK( NM13+1 ), U( 1, LL ), LDU ) 00538 IF( NCC.GT.0 ) 00539 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), 00540 $ WORK( NM13+1 ), C( LL, 1 ), LDC ) 00541 * 00542 * Test convergence 00543 * 00544 IF( ABS( E( M-1 ) ).LE.THRESH ) 00545 $ E( M-1 ) = ZERO 00546 * 00547 ELSE 00548 * 00549 * Chase bulge from bottom to top 00550 * Save cosines and sines for later singular vector updates 00551 * 00552 CS = ONE 00553 OLDCS = ONE 00554 DO 130 I = M, LL + 1, -1 00555 CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R ) 00556 IF( I.LT.M ) 00557 $ E( I ) = OLDSN*R 00558 CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) ) 00559 WORK( I-LL ) = CS 00560 WORK( I-LL+NM1 ) = -SN 00561 WORK( I-LL+NM12 ) = OLDCS 00562 WORK( I-LL+NM13 ) = -OLDSN 00563 130 CONTINUE 00564 H = D( LL )*CS 00565 D( LL ) = H*OLDCS 00566 E( LL ) = H*OLDSN 00567 * 00568 * Update singular vectors 00569 * 00570 IF( NCVT.GT.0 ) 00571 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), 00572 $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) 00573 IF( NRU.GT.0 ) 00574 $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), 00575 $ WORK( N ), U( 1, LL ), LDU ) 00576 IF( NCC.GT.0 ) 00577 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ), 00578 $ WORK( N ), C( LL, 1 ), LDC ) 00579 * 00580 * Test convergence 00581 * 00582 IF( ABS( E( LL ) ).LE.THRESH ) 00583 $ E( LL ) = ZERO 00584 END IF 00585 ELSE 00586 * 00587 * Use nonzero shift 00588 * 00589 IF( IDIR.EQ.1 ) THEN 00590 * 00591 * Chase bulge from top to bottom 00592 * Save cosines and sines for later singular vector updates 00593 * 00594 F = ( ABS( D( LL ) )-SHIFT )* 00595 $ ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) ) 00596 G = E( LL ) 00597 DO 140 I = LL, M - 1 00598 CALL DLARTG( F, G, COSR, SINR, R ) 00599 IF( I.GT.LL ) 00600 $ E( I-1 ) = R 00601 F = COSR*D( I ) + SINR*E( I ) 00602 E( I ) = COSR*E( I ) - SINR*D( I ) 00603 G = SINR*D( I+1 ) 00604 D( I+1 ) = COSR*D( I+1 ) 00605 CALL DLARTG( F, G, COSL, SINL, R ) 00606 D( I ) = R 00607 F = COSL*E( I ) + SINL*D( I+1 ) 00608 D( I+1 ) = COSL*D( I+1 ) - SINL*E( I ) 00609 IF( I.LT.M-1 ) THEN 00610 G = SINL*E( I+1 ) 00611 E( I+1 ) = COSL*E( I+1 ) 00612 END IF 00613 WORK( I-LL+1 ) = COSR 00614 WORK( I-LL+1+NM1 ) = SINR 00615 WORK( I-LL+1+NM12 ) = COSL 00616 WORK( I-LL+1+NM13 ) = SINL 00617 140 CONTINUE 00618 E( M-1 ) = F 00619 * 00620 * Update singular vectors 00621 * 00622 IF( NCVT.GT.0 ) 00623 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), 00624 $ WORK( N ), VT( LL, 1 ), LDVT ) 00625 IF( NRU.GT.0 ) 00626 $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), 00627 $ WORK( NM13+1 ), U( 1, LL ), LDU ) 00628 IF( NCC.GT.0 ) 00629 $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), 00630 $ WORK( NM13+1 ), C( LL, 1 ), LDC ) 00631 * 00632 * Test convergence 00633 * 00634 IF( ABS( E( M-1 ) ).LE.THRESH ) 00635 $ E( M-1 ) = ZERO 00636 * 00637 ELSE 00638 * 00639 * Chase bulge from bottom to top 00640 * Save cosines and sines for later singular vector updates 00641 * 00642 F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT / 00643 $ D( M ) ) 00644 G = E( M-1 ) 00645 DO 150 I = M, LL + 1, -1 00646 CALL DLARTG( F, G, COSR, SINR, R ) 00647 IF( I.LT.M ) 00648 $ E( I ) = R 00649 F = COSR*D( I ) + SINR*E( I-1 ) 00650 E( I-1 ) = COSR*E( I-1 ) - SINR*D( I ) 00651 G = SINR*D( I-1 ) 00652 D( I-1 ) = COSR*D( I-1 ) 00653 CALL DLARTG( F, G, COSL, SINL, R ) 00654 D( I ) = R 00655 F = COSL*E( I-1 ) + SINL*D( I-1 ) 00656 D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 ) 00657 IF( I.GT.LL+1 ) THEN 00658 G = SINL*E( I-2 ) 00659 E( I-2 ) = COSL*E( I-2 ) 00660 END IF 00661 WORK( I-LL ) = COSR 00662 WORK( I-LL+NM1 ) = -SINR 00663 WORK( I-LL+NM12 ) = COSL 00664 WORK( I-LL+NM13 ) = -SINL 00665 150 CONTINUE 00666 E( LL ) = F 00667 * 00668 * Test convergence 00669 * 00670 IF( ABS( E( LL ) ).LE.THRESH ) 00671 $ E( LL ) = ZERO 00672 * 00673 * Update singular vectors if desired 00674 * 00675 IF( NCVT.GT.0 ) 00676 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), 00677 $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) 00678 IF( NRU.GT.0 ) 00679 $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), 00680 $ WORK( N ), U( 1, LL ), LDU ) 00681 IF( NCC.GT.0 ) 00682 $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ), 00683 $ WORK( N ), C( LL, 1 ), LDC ) 00684 END IF 00685 END IF 00686 * 00687 * QR iteration finished, go back and check convergence 00688 * 00689 GO TO 60 00690 * 00691 * All singular values converged, so make them positive 00692 * 00693 160 CONTINUE 00694 DO 170 I = 1, N 00695 IF( D( I ).LT.ZERO ) THEN 00696 D( I ) = -D( I ) 00697 * 00698 * Change sign of singular vectors, if desired 00699 * 00700 IF( NCVT.GT.0 ) 00701 $ CALL DSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT ) 00702 END IF 00703 170 CONTINUE 00704 * 00705 * Sort the singular values into decreasing order (insertion sort on 00706 * singular values, but only one transposition per singular vector) 00707 * 00708 DO 190 I = 1, N - 1 00709 * 00710 * Scan for smallest D(I) 00711 * 00712 ISUB = 1 00713 SMIN = D( 1 ) 00714 DO 180 J = 2, N + 1 - I 00715 IF( D( J ).LE.SMIN ) THEN 00716 ISUB = J 00717 SMIN = D( J ) 00718 END IF 00719 180 CONTINUE 00720 IF( ISUB.NE.N+1-I ) THEN 00721 * 00722 * Swap singular values and vectors 00723 * 00724 D( ISUB ) = D( N+1-I ) 00725 D( N+1-I ) = SMIN 00726 IF( NCVT.GT.0 ) 00727 $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ), 00728 $ LDVT ) 00729 IF( NRU.GT.0 ) 00730 $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 ) 00731 IF( NCC.GT.0 ) 00732 $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC ) 00733 END IF 00734 190 CONTINUE 00735 GO TO 220 00736 * 00737 * Maximum number of iterations exceeded, failure to converge 00738 * 00739 200 CONTINUE 00740 INFO = 0 00741 DO 210 I = 1, N - 1 00742 IF( E( I ).NE.ZERO ) 00743 $ INFO = INFO + 1 00744 210 CONTINUE 00745 220 CONTINUE 00746 RETURN 00747 * 00748 * End of DBDSQR 00749 * 00750 END
1.7.3 
