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