"Fossies" - the Fresh Open Source Software Archive

Member "lapack-3.9.1/SRC/cbdsqr.f" (25 Mar 2021, 25815 Bytes) of package /linux/misc/lapack-3.9.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Fortran 77 source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the latest Fossies "Diffs" side-by-side code changes report for "cbdsqr.f": 3.9.0_vs_3.9.1.

    1 *> \brief \b CBDSQR
    2 *
    3 *  =========== DOCUMENTATION ===========
    4 *
    5 * Online html documentation available at
    6 *            http://www.netlib.org/lapack/explore-html/
    7 *
    8 *> \htmlonly
    9 *> Download CBDSQR + dependencies
   10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cbdsqr.f">
   11 *> [TGZ]</a>
   12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cbdsqr.f">
   13 *> [ZIP]</a>
   14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cbdsqr.f">
   15 *> [TXT]</a>
   16 *> \endhtmlonly
   17 *
   18 *  Definition:
   19 *  ===========
   20 *
   21 *       SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
   22 *                          LDU, C, LDC, RWORK, INFO )
   23 *
   24 *       .. Scalar Arguments ..
   25 *       CHARACTER          UPLO
   26 *       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
   27 *       ..
   28 *       .. Array Arguments ..
   29 *       REAL               D( * ), E( * ), RWORK( * )
   30 *       COMPLEX            C( LDC, * ), U( LDU, * ), VT( LDVT, * )
   31 *       ..
   32 *
   33 *
   34 *> \par Purpose:
   35 *  =============
   36 *>
   37 *> \verbatim
   38 *>
   39 *> CBDSQR computes the singular values and, optionally, the right and/or
   40 *> left singular vectors from the singular value decomposition (SVD) of
   41 *> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
   42 *> zero-shift QR algorithm.  The SVD of B has the form
   43 *>
   44 *>    B = Q * S * P**H
   45 *>
   46 *> where S is the diagonal matrix of singular values, Q is an orthogonal
   47 *> matrix of left singular vectors, and P is an orthogonal matrix of
   48 *> right singular vectors.  If left singular vectors are requested, this
   49 *> subroutine actually returns U*Q instead of Q, and, if right singular
   50 *> vectors are requested, this subroutine returns P**H*VT instead of
   51 *> P**H, for given complex input matrices U and VT.  When U and VT are
   52 *> the unitary matrices that reduce a general matrix A to bidiagonal
   53 *> form: A = U*B*VT, as computed by CGEBRD, then
   54 *>
   55 *>    A = (U*Q) * S * (P**H*VT)
   56 *>
   57 *> is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
   58 *> for a given complex input matrix C.
   59 *>
   60 *> See "Computing  Small Singular Values of Bidiagonal Matrices With
   61 *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
   62 *> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
   63 *> no. 5, pp. 873-912, Sept 1990) and
   64 *> "Accurate singular values and differential qd algorithms," by
   65 *> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
   66 *> Department, University of California at Berkeley, July 1992
   67 *> for a detailed description of the algorithm.
   68 *> \endverbatim
   69 *
   70 *  Arguments:
   71 *  ==========
   72 *
   73 *> \param[in] UPLO
   74 *> \verbatim
   75 *>          UPLO is CHARACTER*1
   76 *>          = 'U':  B is upper bidiagonal;
   77 *>          = 'L':  B is lower bidiagonal.
   78 *> \endverbatim
   79 *>
   80 *> \param[in] N
   81 *> \verbatim
   82 *>          N is INTEGER
   83 *>          The order of the matrix B.  N >= 0.
   84 *> \endverbatim
   85 *>
   86 *> \param[in] NCVT
   87 *> \verbatim
   88 *>          NCVT is INTEGER
   89 *>          The number of columns of the matrix VT. NCVT >= 0.
   90 *> \endverbatim
   91 *>
   92 *> \param[in] NRU
   93 *> \verbatim
   94 *>          NRU is INTEGER
   95 *>          The number of rows of the matrix U. NRU >= 0.
   96 *> \endverbatim
   97 *>
   98 *> \param[in] NCC
   99 *> \verbatim
  100 *>          NCC is INTEGER
  101 *>          The number of columns of the matrix C. NCC >= 0.
  102 *> \endverbatim
  103 *>
  104 *> \param[in,out] D
  105 *> \verbatim
  106 *>          D is REAL array, dimension (N)
  107 *>          On entry, the n diagonal elements of the bidiagonal matrix B.
  108 *>          On exit, if INFO=0, the singular values of B in decreasing
  109 *>          order.
  110 *> \endverbatim
  111 *>
  112 *> \param[in,out] E
  113 *> \verbatim
  114 *>          E is REAL array, dimension (N-1)
  115 *>          On entry, the N-1 offdiagonal elements of the bidiagonal
  116 *>          matrix B.
  117 *>          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
  118 *>          will contain the diagonal and superdiagonal elements of a
  119 *>          bidiagonal matrix orthogonally equivalent to the one given
  120 *>          as input.
  121 *> \endverbatim
  122 *>
  123 *> \param[in,out] VT
  124 *> \verbatim
  125 *>          VT is COMPLEX array, dimension (LDVT, NCVT)
  126 *>          On entry, an N-by-NCVT matrix VT.
  127 *>          On exit, VT is overwritten by P**H * VT.
  128 *>          Not referenced if NCVT = 0.
  129 *> \endverbatim
  130 *>
  131 *> \param[in] LDVT
  132 *> \verbatim
  133 *>          LDVT is INTEGER
  134 *>          The leading dimension of the array VT.
  135 *>          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
  136 *> \endverbatim
  137 *>
  138 *> \param[in,out] U
  139 *> \verbatim
  140 *>          U is COMPLEX array, dimension (LDU, N)
  141 *>          On entry, an NRU-by-N matrix U.
  142 *>          On exit, U is overwritten by U * Q.
  143 *>          Not referenced if NRU = 0.
  144 *> \endverbatim
  145 *>
  146 *> \param[in] LDU
  147 *> \verbatim
  148 *>          LDU is INTEGER
  149 *>          The leading dimension of the array U.  LDU >= max(1,NRU).
  150 *> \endverbatim
  151 *>
  152 *> \param[in,out] C
  153 *> \verbatim
  154 *>          C is COMPLEX array, dimension (LDC, NCC)
  155 *>          On entry, an N-by-NCC matrix C.
  156 *>          On exit, C is overwritten by Q**H * C.
  157 *>          Not referenced if NCC = 0.
  158 *> \endverbatim
  159 *>
  160 *> \param[in] LDC
  161 *> \verbatim
  162 *>          LDC is INTEGER
  163 *>          The leading dimension of the array C.
  164 *>          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
  165 *> \endverbatim
  166 *>
  167 *> \param[out] RWORK
  168 *> \verbatim
  169 *>          RWORK is REAL array, dimension (4*N)
  170 *> \endverbatim
  171 *>
  172 *> \param[out] INFO
  173 *> \verbatim
  174 *>          INFO is INTEGER
  175 *>          = 0:  successful exit
  176 *>          < 0:  If INFO = -i, the i-th argument had an illegal value
  177 *>          > 0:  the algorithm did not converge; D and E contain the
  178 *>                elements of a bidiagonal matrix which is orthogonally
  179 *>                similar to the input matrix B;  if INFO = i, i
  180 *>                elements of E have not converged to zero.
  181 *> \endverbatim
  182 *
  183 *> \par Internal Parameters:
  184 *  =========================
  185 *>
  186 *> \verbatim
  187 *>  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
  188 *>          TOLMUL controls the convergence criterion of the QR loop.
  189 *>          If it is positive, TOLMUL*EPS is the desired relative
  190 *>             precision in the computed singular values.
  191 *>          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
  192 *>             desired absolute accuracy in the computed singular
  193 *>             values (corresponds to relative accuracy
  194 *>             abs(TOLMUL*EPS) in the largest singular value.
  195 *>          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
  196 *>             between 10 (for fast convergence) and .1/EPS
  197 *>             (for there to be some accuracy in the results).
  198 *>          Default is to lose at either one eighth or 2 of the
  199 *>             available decimal digits in each computed singular value
  200 *>             (whichever is smaller).
  201 *>
  202 *>  MAXITR  INTEGER, default = 6
  203 *>          MAXITR controls the maximum number of passes of the
  204 *>          algorithm through its inner loop. The algorithms stops
  205 *>          (and so fails to converge) if the number of passes
  206 *>          through the inner loop exceeds MAXITR*N**2.
  207 *> \endverbatim
  208 *
  209 *  Authors:
  210 *  ========
  211 *
  212 *> \author Univ. of Tennessee
  213 *> \author Univ. of California Berkeley
  214 *> \author Univ. of Colorado Denver
  215 *> \author NAG Ltd.
  216 *
  217 *> \ingroup complexOTHERcomputational
  218 *
  219 *  =====================================================================
  220       SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
  221      $                   LDU, C, LDC, RWORK, INFO )
  222 *
  223 *  -- LAPACK computational routine --
  224 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  225 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  226 *
  227 *     .. Scalar Arguments ..
  228       CHARACTER          UPLO
  229       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
  230 *     ..
  231 *     .. Array Arguments ..
  232       REAL               D( * ), E( * ), RWORK( * )
  233       COMPLEX            C( LDC, * ), U( LDU, * ), VT( LDVT, * )
  234 *     ..
  235 *
  236 *  =====================================================================
  237 *
  238 *     .. Parameters ..
  239       REAL               ZERO
  240       PARAMETER          ( ZERO = 0.0E0 )
  241       REAL               ONE
  242       PARAMETER          ( ONE = 1.0E0 )
  243       REAL               NEGONE
  244       PARAMETER          ( NEGONE = -1.0E0 )
  245       REAL               HNDRTH
  246       PARAMETER          ( HNDRTH = 0.01E0 )
  247       REAL               TEN
  248       PARAMETER          ( TEN = 10.0E0 )
  249       REAL               HNDRD
  250       PARAMETER          ( HNDRD = 100.0E0 )
  251       REAL               MEIGTH
  252       PARAMETER          ( MEIGTH = -0.125E0 )
  253       INTEGER            MAXITR
  254       PARAMETER          ( MAXITR = 6 )
  255 *     ..
  256 *     .. Local Scalars ..
  257       LOGICAL            LOWER, ROTATE
  258       INTEGER            I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
  259      $                   NM12, NM13, OLDLL, OLDM
  260       REAL               ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
  261      $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
  262      $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
  263      $                   SN, THRESH, TOL, TOLMUL, UNFL
  264 *     ..
  265 *     .. External Functions ..
  266       LOGICAL            LSAME
  267       REAL               SLAMCH
  268       EXTERNAL           LSAME, SLAMCH
  269 *     ..
  270 *     .. External Subroutines ..
  271       EXTERNAL           CLASR, CSROT, CSSCAL, CSWAP, SLARTG, SLAS2,
  272      $                   SLASQ1, SLASV2, XERBLA
  273 *     ..
  274 *     .. Intrinsic Functions ..
  275       INTRINSIC          ABS, MAX, MIN, REAL, SIGN, SQRT
  276 *     ..
  277 *     .. Executable Statements ..
  278 *
  279 *     Test the input parameters.
  280 *
  281       INFO = 0
  282       LOWER = LSAME( UPLO, 'L' )
  283       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
  284          INFO = -1
  285       ELSE IF( N.LT.0 ) THEN
  286          INFO = -2
  287       ELSE IF( NCVT.LT.0 ) THEN
  288          INFO = -3
  289       ELSE IF( NRU.LT.0 ) THEN
  290          INFO = -4
  291       ELSE IF( NCC.LT.0 ) THEN
  292          INFO = -5
  293       ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
  294      $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
  295          INFO = -9
  296       ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
  297          INFO = -11
  298       ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
  299      $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
  300          INFO = -13
  301       END IF
  302       IF( INFO.NE.0 ) THEN
  303          CALL XERBLA( 'CBDSQR', -INFO )
  304          RETURN
  305       END IF
  306       IF( N.EQ.0 )
  307      $   RETURN
  308       IF( N.EQ.1 )
  309      $   GO TO 160
  310 *
  311 *     ROTATE is true if any singular vectors desired, false otherwise
  312 *
  313       ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
  314 *
  315 *     If no singular vectors desired, use qd algorithm
  316 *
  317       IF( .NOT.ROTATE ) THEN
  318          CALL SLASQ1( N, D, E, RWORK, INFO )
  319 *
  320 *     If INFO equals 2, dqds didn't finish, try to finish
  321 *
  322          IF( INFO .NE. 2 ) RETURN
  323          INFO = 0
  324       END IF
  325 *
  326       NM1 = N - 1
  327       NM12 = NM1 + NM1
  328       NM13 = NM12 + NM1
  329       IDIR = 0
  330 *
  331 *     Get machine constants
  332 *
  333       EPS = SLAMCH( 'Epsilon' )
  334       UNFL = SLAMCH( 'Safe minimum' )
  335 *
  336 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
  337 *     by applying Givens rotations on the left
  338 *
  339       IF( LOWER ) THEN
  340          DO 10 I = 1, N - 1
  341             CALL SLARTG( D( I ), E( I ), CS, SN, R )
  342             D( I ) = R
  343             E( I ) = SN*D( I+1 )
  344             D( I+1 ) = CS*D( I+1 )
  345             RWORK( I ) = CS
  346             RWORK( NM1+I ) = SN
  347    10    CONTINUE
  348 *
  349 *        Update singular vectors if desired
  350 *
  351          IF( NRU.GT.0 )
  352      $      CALL CLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ),
  353      $                  U, LDU )
  354          IF( NCC.GT.0 )
  355      $      CALL CLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), RWORK( N ),
  356      $                  C, LDC )
  357       END IF
  358 *
  359 *     Compute singular values to relative accuracy TOL
  360 *     (By setting TOL to be negative, algorithm will compute
  361 *     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
  362 *
  363       TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
  364       TOL = TOLMUL*EPS
  365 *
  366 *     Compute approximate maximum, minimum singular values
  367 *
  368       SMAX = ZERO
  369       DO 20 I = 1, N
  370          SMAX = MAX( SMAX, ABS( D( I ) ) )
  371    20 CONTINUE
  372       DO 30 I = 1, N - 1
  373          SMAX = MAX( SMAX, ABS( E( I ) ) )
  374    30 CONTINUE
  375       SMINL = ZERO
  376       IF( TOL.GE.ZERO ) THEN
  377 *
  378 *        Relative accuracy desired
  379 *
  380          SMINOA = ABS( D( 1 ) )
  381          IF( SMINOA.EQ.ZERO )
  382      $      GO TO 50
  383          MU = SMINOA
  384          DO 40 I = 2, N
  385             MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
  386             SMINOA = MIN( SMINOA, MU )
  387             IF( SMINOA.EQ.ZERO )
  388      $         GO TO 50
  389    40    CONTINUE
  390    50    CONTINUE
  391          SMINOA = SMINOA / SQRT( REAL( N ) )
  392          THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
  393       ELSE
  394 *
  395 *        Absolute accuracy desired
  396 *
  397          THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
  398       END IF
  399 *
  400 *     Prepare for main iteration loop for the singular values
  401 *     (MAXIT is the maximum number of passes through the inner
  402 *     loop permitted before nonconvergence signalled.)
  403 *
  404       MAXIT = MAXITR*N*N
  405       ITER = 0
  406       OLDLL = -1
  407       OLDM = -1
  408 *
  409 *     M points to last element of unconverged part of matrix
  410 *
  411       M = N
  412 *
  413 *     Begin main iteration loop
  414 *
  415    60 CONTINUE
  416 *
  417 *     Check for convergence or exceeding iteration count
  418 *
  419       IF( M.LE.1 )
  420      $   GO TO 160
  421       IF( ITER.GT.MAXIT )
  422      $   GO TO 200
  423 *
  424 *     Find diagonal block of matrix to work on
  425 *
  426       IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
  427      $   D( M ) = ZERO
  428       SMAX = ABS( D( M ) )
  429       SMIN = SMAX
  430       DO 70 LLL = 1, M - 1
  431          LL = M - LLL
  432          ABSS = ABS( D( LL ) )
  433          ABSE = ABS( E( LL ) )
  434          IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
  435      $      D( LL ) = ZERO
  436          IF( ABSE.LE.THRESH )
  437      $      GO TO 80
  438          SMIN = MIN( SMIN, ABSS )
  439          SMAX = MAX( SMAX, ABSS, ABSE )
  440    70 CONTINUE
  441       LL = 0
  442       GO TO 90
  443    80 CONTINUE
  444       E( LL ) = ZERO
  445 *
  446 *     Matrix splits since E(LL) = 0
  447 *
  448       IF( LL.EQ.M-1 ) THEN
  449 *
  450 *        Convergence of bottom singular value, return to top of loop
  451 *
  452          M = M - 1
  453          GO TO 60
  454       END IF
  455    90 CONTINUE
  456       LL = LL + 1
  457 *
  458 *     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
  459 *
  460       IF( LL.EQ.M-1 ) THEN
  461 *
  462 *        2 by 2 block, handle separately
  463 *
  464          CALL SLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
  465      $                COSR, SINL, COSL )
  466          D( M-1 ) = SIGMX
  467          E( M-1 ) = ZERO
  468          D( M ) = SIGMN
  469 *
  470 *        Compute singular vectors, if desired
  471 *
  472          IF( NCVT.GT.0 )
  473      $      CALL CSROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT,
  474      $                  COSR, SINR )
  475          IF( NRU.GT.0 )
  476      $      CALL CSROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
  477          IF( NCC.GT.0 )
  478      $      CALL CSROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
  479      $                  SINL )
  480          M = M - 2
  481          GO TO 60
  482       END IF
  483 *
  484 *     If working on new submatrix, choose shift direction
  485 *     (from larger end diagonal element towards smaller)
  486 *
  487       IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
  488          IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
  489 *
  490 *           Chase bulge from top (big end) to bottom (small end)
  491 *
  492             IDIR = 1
  493          ELSE
  494 *
  495 *           Chase bulge from bottom (big end) to top (small end)
  496 *
  497             IDIR = 2
  498          END IF
  499       END IF
  500 *
  501 *     Apply convergence tests
  502 *
  503       IF( IDIR.EQ.1 ) THEN
  504 *
  505 *        Run convergence test in forward direction
  506 *        First apply standard test to bottom of matrix
  507 *
  508          IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
  509      $       ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
  510             E( M-1 ) = ZERO
  511             GO TO 60
  512          END IF
  513 *
  514          IF( TOL.GE.ZERO ) THEN
  515 *
  516 *           If relative accuracy desired,
  517 *           apply convergence criterion forward
  518 *
  519             MU = ABS( D( LL ) )
  520             SMINL = MU
  521             DO 100 LLL = LL, M - 1
  522                IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
  523                   E( LLL ) = ZERO
  524                   GO TO 60
  525                END IF
  526                MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
  527                SMINL = MIN( SMINL, MU )
  528   100       CONTINUE
  529          END IF
  530 *
  531       ELSE
  532 *
  533 *        Run convergence test in backward direction
  534 *        First apply standard test to top of matrix
  535 *
  536          IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
  537      $       ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
  538             E( LL ) = ZERO
  539             GO TO 60
  540          END IF
  541 *
  542          IF( TOL.GE.ZERO ) THEN
  543 *
  544 *           If relative accuracy desired,
  545 *           apply convergence criterion backward
  546 *
  547             MU = ABS( D( M ) )
  548             SMINL = MU
  549             DO 110 LLL = M - 1, LL, -1
  550                IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
  551                   E( LLL ) = ZERO
  552                   GO TO 60
  553                END IF
  554                MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
  555                SMINL = MIN( SMINL, MU )
  556   110       CONTINUE
  557          END IF
  558       END IF
  559       OLDLL = LL
  560       OLDM = M
  561 *
  562 *     Compute shift.  First, test if shifting would ruin relative
  563 *     accuracy, and if so set the shift to zero.
  564 *
  565       IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
  566      $    MAX( EPS, HNDRTH*TOL ) ) THEN
  567 *
  568 *        Use a zero shift to avoid loss of relative accuracy
  569 *
  570          SHIFT = ZERO
  571       ELSE
  572 *
  573 *        Compute the shift from 2-by-2 block at end of matrix
  574 *
  575          IF( IDIR.EQ.1 ) THEN
  576             SLL = ABS( D( LL ) )
  577             CALL SLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
  578          ELSE
  579             SLL = ABS( D( M ) )
  580             CALL SLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
  581          END IF
  582 *
  583 *        Test if shift negligible, and if so set to zero
  584 *
  585          IF( SLL.GT.ZERO ) THEN
  586             IF( ( SHIFT / SLL )**2.LT.EPS )
  587      $         SHIFT = ZERO
  588          END IF
  589       END IF
  590 *
  591 *     Increment iteration count
  592 *
  593       ITER = ITER + M - LL
  594 *
  595 *     If SHIFT = 0, do simplified QR iteration
  596 *
  597       IF( SHIFT.EQ.ZERO ) THEN
  598          IF( IDIR.EQ.1 ) THEN
  599 *
  600 *           Chase bulge from top to bottom
  601 *           Save cosines and sines for later singular vector updates
  602 *
  603             CS = ONE
  604             OLDCS = ONE
  605             DO 120 I = LL, M - 1
  606                CALL SLARTG( D( I )*CS, E( I ), CS, SN, R )
  607                IF( I.GT.LL )
  608      $            E( I-1 ) = OLDSN*R
  609                CALL SLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
  610                RWORK( I-LL+1 ) = CS
  611                RWORK( I-LL+1+NM1 ) = SN
  612                RWORK( I-LL+1+NM12 ) = OLDCS
  613                RWORK( I-LL+1+NM13 ) = OLDSN
  614   120       CONTINUE
  615             H = D( M )*CS
  616             D( M ) = H*OLDCS
  617             E( M-1 ) = H*OLDSN
  618 *
  619 *           Update singular vectors
  620 *
  621             IF( NCVT.GT.0 )
  622      $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
  623      $                     RWORK( N ), VT( LL, 1 ), LDVT )
  624             IF( NRU.GT.0 )
  625      $         CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ),
  626      $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
  627             IF( NCC.GT.0 )
  628      $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ),
  629      $                     RWORK( NM13+1 ), C( LL, 1 ), LDC )
  630 *
  631 *           Test convergence
  632 *
  633             IF( ABS( E( M-1 ) ).LE.THRESH )
  634      $         E( M-1 ) = ZERO
  635 *
  636          ELSE
  637 *
  638 *           Chase bulge from bottom to top
  639 *           Save cosines and sines for later singular vector updates
  640 *
  641             CS = ONE
  642             OLDCS = ONE
  643             DO 130 I = M, LL + 1, -1
  644                CALL SLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
  645                IF( I.LT.M )
  646      $            E( I ) = OLDSN*R
  647                CALL SLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
  648                RWORK( I-LL ) = CS
  649                RWORK( I-LL+NM1 ) = -SN
  650                RWORK( I-LL+NM12 ) = OLDCS
  651                RWORK( I-LL+NM13 ) = -OLDSN
  652   130       CONTINUE
  653             H = D( LL )*CS
  654             D( LL ) = H*OLDCS
  655             E( LL ) = H*OLDSN
  656 *
  657 *           Update singular vectors
  658 *
  659             IF( NCVT.GT.0 )
  660      $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ),
  661      $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
  662             IF( NRU.GT.0 )
  663      $         CALL CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
  664      $                     RWORK( N ), U( 1, LL ), LDU )
  665             IF( NCC.GT.0 )
  666      $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ),
  667      $                     RWORK( N ), C( LL, 1 ), LDC )
  668 *
  669 *           Test convergence
  670 *
  671             IF( ABS( E( LL ) ).LE.THRESH )
  672      $         E( LL ) = ZERO
  673          END IF
  674       ELSE
  675 *
  676 *        Use nonzero shift
  677 *
  678          IF( IDIR.EQ.1 ) THEN
  679 *
  680 *           Chase bulge from top to bottom
  681 *           Save cosines and sines for later singular vector updates
  682 *
  683             F = ( ABS( D( LL ) )-SHIFT )*
  684      $          ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
  685             G = E( LL )
  686             DO 140 I = LL, M - 1
  687                CALL SLARTG( F, G, COSR, SINR, R )
  688                IF( I.GT.LL )
  689      $            E( I-1 ) = R
  690                F = COSR*D( I ) + SINR*E( I )
  691                E( I ) = COSR*E( I ) - SINR*D( I )
  692                G = SINR*D( I+1 )
  693                D( I+1 ) = COSR*D( I+1 )
  694                CALL SLARTG( F, G, COSL, SINL, R )
  695                D( I ) = R
  696                F = COSL*E( I ) + SINL*D( I+1 )
  697                D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
  698                IF( I.LT.M-1 ) THEN
  699                   G = SINL*E( I+1 )
  700                   E( I+1 ) = COSL*E( I+1 )
  701                END IF
  702                RWORK( I-LL+1 ) = COSR
  703                RWORK( I-LL+1+NM1 ) = SINR
  704                RWORK( I-LL+1+NM12 ) = COSL
  705                RWORK( I-LL+1+NM13 ) = SINL
  706   140       CONTINUE
  707             E( M-1 ) = F
  708 *
  709 *           Update singular vectors
  710 *
  711             IF( NCVT.GT.0 )
  712      $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
  713      $                     RWORK( N ), VT( LL, 1 ), LDVT )
  714             IF( NRU.GT.0 )
  715      $         CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ),
  716      $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
  717             IF( NCC.GT.0 )
  718      $         CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ),
  719      $                     RWORK( NM13+1 ), C( LL, 1 ), LDC )
  720 *
  721 *           Test convergence
  722 *
  723             IF( ABS( E( M-1 ) ).LE.THRESH )
  724      $         E( M-1 ) = ZERO
  725 *
  726          ELSE
  727 *
  728 *           Chase bulge from bottom to top
  729 *           Save cosines and sines for later singular vector updates
  730 *
  731             F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
  732      $          D( M ) )
  733             G = E( M-1 )
  734             DO 150 I = M, LL + 1, -1
  735                CALL SLARTG( F, G, COSR, SINR, R )
  736                IF( I.LT.M )
  737      $            E( I ) = R
  738                F = COSR*D( I ) + SINR*E( I-1 )
  739                E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
  740                G = SINR*D( I-1 )
  741                D( I-1 ) = COSR*D( I-1 )
  742                CALL SLARTG( F, G, COSL, SINL, R )
  743                D( I ) = R
  744                F = COSL*E( I-1 ) + SINL*D( I-1 )
  745                D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
  746                IF( I.GT.LL+1 ) THEN
  747                   G = SINL*E( I-2 )
  748                   E( I-2 ) = COSL*E( I-2 )
  749                END IF
  750                RWORK( I-LL ) = COSR
  751                RWORK( I-LL+NM1 ) = -SINR
  752                RWORK( I-LL+NM12 ) = COSL
  753                RWORK( I-LL+NM13 ) = -SINL
  754   150       CONTINUE
  755             E( LL ) = F
  756 *
  757 *           Test convergence
  758 *
  759             IF( ABS( E( LL ) ).LE.THRESH )
  760      $         E( LL ) = ZERO
  761 *
  762 *           Update singular vectors if desired
  763 *
  764             IF( NCVT.GT.0 )
  765      $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ),
  766      $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
  767             IF( NRU.GT.0 )
  768      $         CALL CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
  769      $                     RWORK( N ), U( 1, LL ), LDU )
  770             IF( NCC.GT.0 )
  771      $         CALL CLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ),
  772      $                     RWORK( N ), C( LL, 1 ), LDC )
  773          END IF
  774       END IF
  775 *
  776 *     QR iteration finished, go back and check convergence
  777 *
  778       GO TO 60
  779 *
  780 *     All singular values converged, so make them positive
  781 *
  782   160 CONTINUE
  783       DO 170 I = 1, N
  784          IF( D( I ).LT.ZERO ) THEN
  785             D( I ) = -D( I )
  786 *
  787 *           Change sign of singular vectors, if desired
  788 *
  789             IF( NCVT.GT.0 )
  790      $         CALL CSSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
  791          END IF
  792   170 CONTINUE
  793 *
  794 *     Sort the singular values into decreasing order (insertion sort on
  795 *     singular values, but only one transposition per singular vector)
  796 *
  797       DO 190 I = 1, N - 1
  798 *
  799 *        Scan for smallest D(I)
  800 *
  801          ISUB = 1
  802          SMIN = D( 1 )
  803          DO 180 J = 2, N + 1 - I
  804             IF( D( J ).LE.SMIN ) THEN
  805                ISUB = J
  806                SMIN = D( J )
  807             END IF
  808   180    CONTINUE
  809          IF( ISUB.NE.N+1-I ) THEN
  810 *
  811 *           Swap singular values and vectors
  812 *
  813             D( ISUB ) = D( N+1-I )
  814             D( N+1-I ) = SMIN
  815             IF( NCVT.GT.0 )
  816      $         CALL CSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
  817      $                     LDVT )
  818             IF( NRU.GT.0 )
  819      $         CALL CSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
  820             IF( NCC.GT.0 )
  821      $         CALL CSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
  822          END IF
  823   190 CONTINUE
  824       GO TO 220
  825 *
  826 *     Maximum number of iterations exceeded, failure to converge
  827 *
  828   200 CONTINUE
  829       INFO = 0
  830       DO 210 I = 1, N - 1
  831          IF( E( I ).NE.ZERO )
  832      $      INFO = INFO + 1
  833   210 CONTINUE
  834   220 CONTINUE
  835       RETURN
  836 *
  837 *     End of CBDSQR
  838 *
  839       END