LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dbdsqr.f
Go to the documentation of this file.
00001 *> \brief \b DBDSQR
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DBDSQR + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DBDSQR( 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 *       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
00030 *      $                   VT( LDVT, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DBDSQR 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 DGEBRD, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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  DOUBLE PRECISION, 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 DBDSQR( 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       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
00244      $                   VT( LDVT, * ), WORK( * )
00245 *     ..
00246 *
00247 *  =====================================================================
00248 *
00249 *     .. Parameters ..
00250       DOUBLE PRECISION   ZERO
00251       PARAMETER          ( ZERO = 0.0D0 )
00252       DOUBLE PRECISION   ONE
00253       PARAMETER          ( ONE = 1.0D0 )
00254       DOUBLE PRECISION   NEGONE
00255       PARAMETER          ( NEGONE = -1.0D0 )
00256       DOUBLE PRECISION   HNDRTH
00257       PARAMETER          ( HNDRTH = 0.01D0 )
00258       DOUBLE PRECISION   TEN
00259       PARAMETER          ( TEN = 10.0D0 )
00260       DOUBLE PRECISION   HNDRD
00261       PARAMETER          ( HNDRD = 100.0D0 )
00262       DOUBLE PRECISION   MEIGTH
00263       PARAMETER          ( MEIGTH = -0.125D0 )
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       DOUBLE PRECISION   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       DOUBLE PRECISION   DLAMCH
00279       EXTERNAL           LSAME, DLAMCH
00280 *     ..
00281 *     .. External Subroutines ..
00282       EXTERNAL           DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT,
00283      $                   DSCAL, DSWAP, XERBLA
00284 *     ..
00285 *     .. Intrinsic Functions ..
00286       INTRINSIC          ABS, DBLE, MAX, MIN, 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( 'DBDSQR', -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 DLASQ1( 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 = DLAMCH( 'Epsilon' )
00345       UNFL = DLAMCH( '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 DLARTG( 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 DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U,
00364      $                  LDU )
00365          IF( NCC.GT.0 )
00366      $      CALL DLASR( '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( DBLE( 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 DLASV2( 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 DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,
00485      $                 SINR )
00486          IF( NRU.GT.0 )
00487      $      CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
00488          IF( NCC.GT.0 )
00489      $      CALL DROT( 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 DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
00589          ELSE
00590             SLL = ABS( D( M ) )
00591             CALL DLAS2( 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 DLARTG( D( I )*CS, E( I ), CS, SN, R )
00618                IF( I.GT.LL )
00619      $            E( I-1 ) = OLDSN*R
00620                CALL DLARTG( 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 DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
00634      $                     WORK( N ), VT( LL, 1 ), LDVT )
00635             IF( NRU.GT.0 )
00636      $         CALL DLASR( '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 DLASR( '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 DLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
00656                IF( I.LT.M )
00657      $            E( I ) = OLDSN*R
00658                CALL DLARTG( 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 DLASR( '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 DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
00675      $                     WORK( N ), U( 1, LL ), LDU )
00676             IF( NCC.GT.0 )
00677      $         CALL DLASR( '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 DLARTG( 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 DLARTG( 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 DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
00724      $                     WORK( N ), VT( LL, 1 ), LDVT )
00725             IF( NRU.GT.0 )
00726      $         CALL DLASR( '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 DLASR( '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 DLARTG( 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 DLARTG( 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 DLASR( '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 DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
00780      $                     WORK( N ), U( 1, LL ), LDU )
00781             IF( NCC.GT.0 )
00782      $         CALL DLASR( '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 DSCAL( 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 DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
00828      $                     LDVT )
00829             IF( NRU.GT.0 )
00830      $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
00831             IF( NCC.GT.0 )
00832      $         CALL DSWAP( 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 DBDSQR
00849 *
00850       END
 All Files Functions