LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cbdsqr.f
Go to the documentation of this file.
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
 All Files Functions