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