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