LAPACK: SRC/zbbcsd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
00002      $                   THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
00003      $                   V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
00004      $                   B22D, B22E, RWORK, LRWORK, INFO )
00005       IMPLICIT NONE
00006 *
00007 *  -- LAPACK routine (version 3.3.0) --
00008 *
00009 *  -- Contributed by Brian Sutton of the Randolph-Macon College --
00010 *  -- November 2010
00011 *
00012 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00013 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--     
00014 *
00015 *     .. Scalar Arguments ..
00016       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
00017       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
00018 *     ..
00019 *     .. Array Arguments ..
00020       DOUBLE PRECISION   B11D( * ), B11E( * ), B12D( * ), B12E( * ),
00021      $                   B21D( * ), B21E( * ), B22D( * ), B22E( * ),
00022      $                   PHI( * ), THETA( * ), RWORK( * )
00023       COMPLEX*16         U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00024      $                   V2T( LDV2T, * )
00025 *     ..
00026 *
00027 *  Purpose
00028 *  =======
00029 *
00030 *  ZBBCSD computes the CS decomposition of a unitary matrix in
00031 *  bidiagonal-block form,
00032 *
00033 *
00034 *      [ B11 | B12 0  0 ]
00035 *      [  0  |  0 -I  0 ]
00036 *  X = [----------------]
00037 *      [ B21 | B22 0  0 ]
00038 *      [  0  |  0  0  I ]
00039 *
00040 *                                [  C | -S  0  0 ]
00041 *                    [ U1 |    ] [  0 |  0 -I  0 ] [ V1 |    ]**H
00042 *                  = [---------] [---------------] [---------]   .
00043 *                    [    | U2 ] [  S |  C  0  0 ] [    | V2 ]
00044 *                                [  0 |  0  0  I ]
00045 *
00046 *  X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
00047 *  than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
00048 *  transposed and/or permuted. This can be done in constant time using
00049 *  the TRANS and SIGNS options. See ZUNCSD for details.)
00050 *
00051 *  The bidiagonal matrices B11, B12, B21, and B22 are represented
00052 *  implicitly by angles THETA(1:Q) and PHI(1:Q-1).
00053 *
00054 *  The unitary matrices U1, U2, V1T, and V2T are input/output.
00055 *  The input matrices are pre- or post-multiplied by the appropriate
00056 *  singular vector matrices.
00057 *
00058 *  Arguments
00059 *  =========
00060 *
00061 *  JOBU1   (input) CHARACTER
00062 *          = 'Y':      U1 is updated;
00063 *          otherwise:  U1 is not updated.
00064 *
00065 *  JOBU2   (input) CHARACTER
00066 *          = 'Y':      U2 is updated;
00067 *          otherwise:  U2 is not updated.
00068 *
00069 *  JOBV1T  (input) CHARACTER
00070 *          = 'Y':      V1T is updated;
00071 *          otherwise:  V1T is not updated.
00072 *
00073 *  JOBV2T  (input) CHARACTER
00074 *          = 'Y':      V2T is updated;
00075 *          otherwise:  V2T is not updated.
00076 *
00077 *  TRANS   (input) CHARACTER
00078 *          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
00079 *                      order;
00080 *          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
00081 *                      major order.
00082 *
00083 *  M       (input) INTEGER
00084 *          The number of rows and columns in X, the unitary matrix in
00085 *          bidiagonal-block form.
00086 *
00087 *  P       (input) INTEGER
00088 *          The number of rows in the top-left block of X. 0 <= P <= M.
00089 *
00090 *  Q       (input) INTEGER
00091 *          The number of columns in the top-left block of X.
00092 *          0 <= Q <= MIN(P,M-P,M-Q).
00093 *
00094 *  THETA   (input/output) DOUBLE PRECISION array, dimension (Q)
00095 *          On entry, the angles THETA(1),...,THETA(Q) that, along with
00096 *          PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
00097 *          form. On exit, the angles whose cosines and sines define the
00098 *          diagonal blocks in the CS decomposition.
00099 *
00100 *  PHI     (input/workspace) DOUBLE PRECISION array, dimension (Q-1)
00101 *          The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
00102 *          THETA(Q), define the matrix in bidiagonal-block form.
00103 *
00104 *  U1      (input/output) COMPLEX*16 array, dimension (LDU1,P)
00105 *          On entry, an LDU1-by-P matrix. On exit, U1 is postmultiplied
00106 *          by the left singular vector matrix common to [ B11 ; 0 ] and
00107 *          [ B12 0 0 ; 0 -I 0 0 ].
00108 *
00109 *  LDU1    (input) INTEGER
00110 *          The leading dimension of the array U1.
00111 *
00112 *  U2      (input/output) COMPLEX*16 array, dimension (LDU2,M-P)
00113 *          On entry, an LDU2-by-(M-P) matrix. On exit, U2 is
00114 *          postmultiplied by the left singular vector matrix common to
00115 *          [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
00116 *
00117 *  LDU2    (input) INTEGER
00118 *          The leading dimension of the array U2.
00119 *
00120 *  V1T     (input/output) COMPLEX*16 array, dimension (LDV1T,Q)
00121 *          On entry, a LDV1T-by-Q matrix. On exit, V1T is premultiplied
00122 *          by the conjugate transpose of the right singular vector
00123 *          matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
00124 *
00125 *  LDV1T   (input) INTEGER
00126 *          The leading dimension of the array V1T.
00127 *
00128 *  V2T     (input/output) COMPLEX*16 array, dimenison (LDV2T,M-Q)
00129 *          On entry, a LDV2T-by-(M-Q) matrix. On exit, V2T is
00130 *          premultiplied by the conjugate transpose of the right
00131 *          singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
00132 *          [ B22 0 0 ; 0 0 I ].
00133 *
00134 *  LDV2T   (input) INTEGER
00135 *          The leading dimension of the array V2T.
00136 *
00137 *  B11D    (output) DOUBLE PRECISION array, dimension (Q)
00138 *          When ZBBCSD converges, B11D contains the cosines of THETA(1),
00139 *          ..., THETA(Q). If ZBBCSD fails to converge, then B11D
00140 *          contains the diagonal of the partially reduced top-left
00141 *          block.
00142 *
00143 *  B11E    (output) DOUBLE PRECISION array, dimension (Q-1)
00144 *          When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails
00145 *          to converge, then B11E contains the superdiagonal of the
00146 *          partially reduced top-left block.
00147 *
00148 *  B12D    (output) DOUBLE PRECISION array, dimension (Q)
00149 *          When ZBBCSD converges, B12D contains the negative sines of
00150 *          THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
00151 *          B12D contains the diagonal of the partially reduced top-right
00152 *          block.
00153 *
00154 *  B12E    (output) DOUBLE PRECISION array, dimension (Q-1)
00155 *          When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails
00156 *          to converge, then B12E contains the subdiagonal of the
00157 *          partially reduced top-right block.
00158 *
00159 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00160 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00161 *
00162 *  LRWORK  (input) INTEGER
00163 *          The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).
00164 *
00165 *          If LRWORK = -1, then a workspace query is assumed; the
00166 *          routine only calculates the optimal size of the RWORK array,
00167 *          returns this value as the first entry of the work array, and
00168 *          no error message related to LRWORK is issued by XERBLA.
00169 *
00170 *  INFO    (output) INTEGER
00171 *          = 0:  successful exit.
00172 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00173 *          > 0:  if ZBBCSD did not converge, INFO specifies the number
00174 *                of nonzero entries in PHI, and B11D, B11E, etc.,
00175 *                contain the partially reduced matrix.
00176 *
00177 *  Reference
00178 *  =========
00179 *
00180 *  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
00181 *      Algorithms, 50(1):33-65, 2009.
00182 *
00183 *  Internal Parameters
00184 *  ===================
00185 *
00186 *  TOLMUL  DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))
00187 *          TOLMUL controls the convergence criterion of the QR loop.
00188 *          Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
00189 *          are within TOLMUL*EPS of either bound.
00190 *
00191 *  ===================================================================
00192 *
00193 *     .. Parameters ..
00194       INTEGER            MAXITR
00195       PARAMETER          ( MAXITR = 6 )
00196       DOUBLE PRECISION   HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO
00197       PARAMETER          ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0,
00198      $                     ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0,
00199      $                     TEN = 10.0D0, ZERO = 0.0D0 )
00200       COMPLEX*16         NEGONECOMPLEX
00201       PARAMETER          ( NEGONECOMPLEX = (-1.0D0,0.0D0) )
00202 *     ..
00203 *     .. Local Scalars ..
00204       LOGICAL            COLMAJOR, LQUERY, RESTART11, RESTART12,
00205      $                   RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
00206      $                   WANTV2T
00207       INTEGER            I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
00208      $                   IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
00209      $                   LRWORKMIN, LRWORKOPT, MAXIT, MINI
00210       DOUBLE PRECISION   B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
00211      $                   EPS, MU, NU, R, SIGMA11, SIGMA21,
00212      $                   TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
00213      $                   UNFL, X1, X2, Y1, Y2
00214 *
00215       EXTERNAL           DLARTGP, DLARTGS, DLAS2, XERBLA, ZLASR, ZSCAL,
00216      $                   ZSWAP
00217 *     ..
00218 *     .. External Functions ..
00219       DOUBLE PRECISION   DLAMCH
00220       LOGICAL            LSAME
00221       EXTERNAL           LSAME, DLAMCH
00222 *     ..
00223 *     .. Intrinsic Functions ..
00224       INTRINSIC          ABS, ATAN2, COS, MAX, MIN, SIN, SQRT
00225 *     ..
00226 *     .. Executable Statements ..
00227 *
00228 *     Test input arguments
00229 *
00230       INFO = 0
00231       LQUERY = LRWORK .EQ. -1
00232       WANTU1 = LSAME( JOBU1, 'Y' )
00233       WANTU2 = LSAME( JOBU2, 'Y' )
00234       WANTV1T = LSAME( JOBV1T, 'Y' )
00235       WANTV2T = LSAME( JOBV2T, 'Y' )
00236       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00237 *
00238       IF( M .LT. 0 ) THEN
00239          INFO = -6
00240       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00241          INFO = -7
00242       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
00243          INFO = -8
00244       ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN
00245          INFO = -8
00246       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
00247          INFO = -12
00248       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
00249          INFO = -14
00250       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
00251          INFO = -16
00252       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
00253          INFO = -18
00254       END IF
00255 *
00256 *     Quick return if Q = 0
00257 *
00258       IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN
00259          LRWORKMIN = 1
00260          RWORK(1) = LRWORKMIN
00261          RETURN
00262       END IF
00263 *
00264 *     Compute workspace
00265 *
00266       IF( INFO .EQ. 0 ) THEN
00267          IU1CS = 1
00268          IU1SN = IU1CS + Q
00269          IU2CS = IU1SN + Q
00270          IU2SN = IU2CS + Q
00271          IV1TCS = IU2SN + Q
00272          IV1TSN = IV1TCS + Q
00273          IV2TCS = IV1TSN + Q
00274          IV2TSN = IV2TCS + Q
00275          LRWORKOPT = IV2TSN + Q - 1
00276          LRWORKMIN = LRWORKOPT
00277          RWORK(1) = LRWORKOPT
00278          IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN
00279             INFO = -28
00280          END IF
00281       END IF
00282 *
00283       IF( INFO .NE. 0 ) THEN
00284          CALL XERBLA( 'ZBBCSD', -INFO )
00285          RETURN
00286       ELSE IF( LQUERY ) THEN
00287          RETURN
00288       END IF
00289 *
00290 *     Get machine constants
00291 *
00292       EPS = DLAMCH( 'Epsilon' )
00293       UNFL = DLAMCH( 'Safe minimum' )
00294       TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) )
00295       TOL = TOLMUL*EPS
00296       THRESH = MAX( TOL, MAXITR*Q*Q*UNFL )
00297 *
00298 *     Test for negligible sines or cosines
00299 *
00300       DO I = 1, Q
00301          IF( THETA(I) .LT. THRESH ) THEN
00302             THETA(I) = ZERO
00303          ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
00304             THETA(I) = PIOVER2
00305          END IF
00306       END DO
00307       DO I = 1, Q-1
00308          IF( PHI(I) .LT. THRESH ) THEN
00309             PHI(I) = ZERO
00310          ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
00311             PHI(I) = PIOVER2
00312          END IF
00313       END DO
00314 *
00315 *     Initial deflation
00316 *
00317       IMAX = Q
00318       DO WHILE( ( IMAX .GT. 1 ) .AND. ( PHI(IMAX-1) .EQ. ZERO ) )
00319          IMAX = IMAX - 1
00320       END DO
00321       IMIN = IMAX - 1
00322       IF  ( IMIN .GT. 1 ) THEN
00323          DO WHILE( PHI(IMIN-1) .NE. ZERO )
00324             IMIN = IMIN - 1
00325             IF  ( IMIN .LE. 1 ) EXIT
00326          END DO
00327       END IF
00328 *
00329 *     Initialize iteration counter
00330 *
00331       MAXIT = MAXITR*Q*Q
00332       ITER = 0
00333 *
00334 *     Begin main iteration loop
00335 *
00336       DO WHILE( IMAX .GT. 1 )
00337 *
00338 *        Compute the matrix entries
00339 *
00340          B11D(IMIN) = COS( THETA(IMIN) )
00341          B21D(IMIN) = -SIN( THETA(IMIN) )
00342          DO I = IMIN, IMAX - 1
00343             B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) )
00344             B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) )
00345             B12D(I) = SIN( THETA(I) ) * COS( PHI(I) )
00346             B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) )
00347             B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) )
00348             B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) )
00349             B22D(I) = COS( THETA(I) ) * COS( PHI(I) )
00350             B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) )
00351          END DO
00352          B12D(IMAX) = SIN( THETA(IMAX) )
00353          B22D(IMAX) = COS( THETA(IMAX) )
00354 *
00355 *        Abort if not converging; otherwise, increment ITER
00356 *
00357          IF( ITER .GT. MAXIT ) THEN
00358             INFO = 0
00359             DO I = 1, Q
00360                IF( PHI(I) .NE. ZERO )
00361      $            INFO = INFO + 1
00362             END DO
00363             RETURN
00364          END IF
00365 *
00366          ITER = ITER + IMAX - IMIN
00367 *
00368 *        Compute shifts
00369 *
00370          THETAMAX = THETA(IMIN)
00371          THETAMIN = THETA(IMIN)
00372          DO I = IMIN+1, IMAX
00373             IF( THETA(I) > THETAMAX )
00374      $         THETAMAX = THETA(I)
00375             IF( THETA(I) < THETAMIN )
00376      $         THETAMIN = THETA(I)
00377          END DO
00378 *
00379          IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN
00380 *
00381 *           Zero on diagonals of B11 and B22; induce deflation with a
00382 *           zero shift
00383 *
00384             MU = ZERO
00385             NU = ONE
00386 *
00387          ELSE IF( THETAMIN .LT. THRESH ) THEN
00388 *
00389 *           Zero on diagonals of B12 and B22; induce deflation with a
00390 *           zero shift
00391 *
00392             MU = ONE
00393             NU = ZERO
00394 *
00395          ELSE
00396 *
00397 *           Compute shifts for B11 and B21 and use the lesser
00398 *
00399             CALL DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11,
00400      $                  DUMMY )
00401             CALL DLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21,
00402      $                  DUMMY )
00403 *
00404             IF( SIGMA11 .LE. SIGMA21 ) THEN
00405                MU = SIGMA11
00406                NU = SQRT( ONE - MU**2 )
00407                IF( MU .LT. THRESH ) THEN
00408                   MU = ZERO
00409                   NU = ONE
00410                END IF
00411             ELSE
00412                NU = SIGMA21
00413                MU = SQRT( 1.0 - NU**2 )
00414                IF( NU .LT. THRESH ) THEN
00415                   MU = ONE
00416                   NU = ZERO
00417                END IF
00418             END IF
00419          END IF
00420 *
00421 *        Rotate to produce bulges in B11 and B21
00422 *
00423          IF( MU .LE. NU ) THEN
00424             CALL DLARTGS( B11D(IMIN), B11E(IMIN), MU,
00425      $                    RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
00426          ELSE
00427             CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU,
00428      $                    RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
00429          END IF
00430 *
00431          TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) +
00432      $          RWORK(IV1TSN+IMIN-1)*B11E(IMIN)
00433          B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) -
00434      $                RWORK(IV1TSN+IMIN-1)*B11D(IMIN)
00435          B11D(IMIN) = TEMP
00436          B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1)
00437          B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1)
00438          TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) +
00439      $          RWORK(IV1TSN+IMIN-1)*B21E(IMIN)
00440          B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) -
00441      $                RWORK(IV1TSN+IMIN-1)*B21D(IMIN)
00442          B21D(IMIN) = TEMP
00443          B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1)
00444          B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1)
00445 *
00446 *        Compute THETA(IMIN)
00447 *
00448          THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ),
00449      $                   SQRT( B11D(IMIN)**2+B11BULGE**2 ) )
00450 *
00451 *        Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
00452 *
00453          IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
00454             CALL DLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1),
00455      $                    RWORK(IU1CS+IMIN-1), R )
00456          ELSE IF( MU .LE. NU ) THEN
00457             CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU,
00458      $                    RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
00459          ELSE
00460             CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
00461      $                    RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
00462          END IF
00463          IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
00464             CALL DLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1),
00465      $                    RWORK(IU2CS+IMIN-1), R )
00466          ELSE IF( NU .LT. MU ) THEN
00467             CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU,
00468      $                    RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
00469          ELSE
00470             CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU,
00471      $                    RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
00472          END IF
00473          RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1)
00474          RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1)
00475 *
00476          TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) +
00477      $          RWORK(IU1SN+IMIN-1)*B11D(IMIN+1)
00478          B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) -
00479      $                  RWORK(IU1SN+IMIN-1)*B11E(IMIN)
00480          B11E(IMIN) = TEMP
00481          IF( IMAX .GT. IMIN+1 ) THEN
00482             B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1)
00483             B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1)
00484          END IF
00485          TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) +
00486      $          RWORK(IU1SN+IMIN-1)*B12E(IMIN)
00487          B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) -
00488      $                RWORK(IU1SN+IMIN-1)*B12D(IMIN)
00489          B12D(IMIN) = TEMP
00490          B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1)
00491          B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1)
00492          TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) +
00493      $          RWORK(IU2SN+IMIN-1)*B21D(IMIN+1)
00494          B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) -
00495      $                  RWORK(IU2SN+IMIN-1)*B21E(IMIN)
00496          B21E(IMIN) = TEMP
00497          IF( IMAX .GT. IMIN+1 ) THEN
00498             B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1)
00499             B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1)
00500          END IF
00501          TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) +
00502      $          RWORK(IU2SN+IMIN-1)*B22E(IMIN)
00503          B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) -
00504      $                RWORK(IU2SN+IMIN-1)*B22D(IMIN)
00505          B22D(IMIN) = TEMP
00506          B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1)
00507          B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1)
00508 *
00509 *        Inner loop: chase bulges from B11(IMIN,IMIN+2),
00510 *        B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
00511 *        bottom-right
00512 *
00513          DO I = IMIN+1, IMAX-1
00514 *
00515 *           Compute PHI(I-1)
00516 *
00517             X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1)
00518             X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE
00519             Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1)
00520             Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE
00521 *
00522             PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) )
00523 *
00524 *           Determine if there are bulges to chase or if a new direct
00525 *           summand has been reached
00526 *
00527             RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
00528             RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
00529             RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
00530             RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
00531 *
00532 *           If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
00533 *           B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
00534 *           chasing by applying the original shift again.
00535 *
00536             IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN
00537                CALL DLARTGP( X2, X1, RWORK(IV1TSN+I-1),
00538      $                       RWORK(IV1TCS+I-1), R )
00539             ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN
00540                CALL DLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1),
00541      $                       RWORK(IV1TCS+I-1), R )
00542             ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN
00543                CALL DLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1),
00544      $                       RWORK(IV1TCS+I-1), R )
00545             ELSE IF( MU .LE. NU ) THEN
00546                CALL DLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1),
00547      $                       RWORK(IV1TSN+I-1) )
00548             ELSE
00549                CALL DLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1),
00550      $                       RWORK(IV1TSN+I-1) )
00551             END IF
00552             RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1)
00553             RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1)
00554             IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
00555                CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1),
00556      $                       RWORK(IV2TCS+I-1-1), R )
00557             ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
00558                CALL DLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1),
00559      $                       RWORK(IV2TCS+I-1-1), R )
00560             ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
00561                CALL DLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1),
00562      $                       RWORK(IV2TCS+I-1-1), R )
00563             ELSE IF( NU .LT. MU ) THEN
00564                CALL DLARTGS( B12E(I-1), B12D(I), NU,
00565      $                       RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
00566             ELSE
00567                CALL DLARTGS( B22E(I-1), B22D(I), MU,
00568      $                       RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
00569             END IF
00570 *
00571             TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I)
00572             B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) -
00573      $                RWORK(IV1TSN+I-1)*B11D(I)
00574             B11D(I) = TEMP
00575             B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1)
00576             B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1)
00577             TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I)
00578             B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) -
00579      $                RWORK(IV1TSN+I-1)*B21D(I)
00580             B21D(I) = TEMP
00581             B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1)
00582             B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1)
00583             TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) +
00584      $             RWORK(IV2TSN+I-1-1)*B12D(I)
00585             B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) -
00586      $                RWORK(IV2TSN+I-1-1)*B12E(I-1)
00587             B12E(I-1) = TEMP
00588             B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I)
00589             B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I)
00590             TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) +
00591      $             RWORK(IV2TSN+I-1-1)*B22D(I)
00592             B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) -
00593      $                RWORK(IV2TSN+I-1-1)*B22E(I-1)
00594             B22E(I-1) = TEMP
00595             B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I)
00596             B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I)
00597 *
00598 *           Compute THETA(I)
00599 *
00600             X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1)
00601             X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE
00602             Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1)
00603             Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE
00604 *
00605             THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) )
00606 *
00607 *           Determine if there are bulges to chase or if a new direct
00608 *           summand has been reached
00609 *
00610             RESTART11 =   B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
00611             RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
00612             RESTART21 =   B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
00613             RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
00614 *
00615 *           If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
00616 *           B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
00617 *           chasing by applying the original shift again.
00618 *
00619             IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN
00620                CALL DLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1),
00621      $                       R )
00622             ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN
00623                CALL DLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1),
00624      $                       RWORK(IU1CS+I-1), R )
00625             ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN
00626                CALL DLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1),
00627      $                       RWORK(IU1CS+I-1), R )
00628             ELSE IF( MU .LE. NU ) THEN
00629                CALL DLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1),
00630      $                       RWORK(IU1SN+I-1) )
00631             ELSE
00632                CALL DLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1),
00633      $                       RWORK(IU1SN+I-1) )
00634             END IF
00635             IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN
00636                CALL DLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1),
00637      $                       R )
00638             ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN
00639                CALL DLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1),
00640      $                       RWORK(IU2CS+I-1), R )
00641             ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN
00642                CALL DLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1),
00643      $                       RWORK(IU2CS+I-1), R )
00644             ELSE IF( NU .LT. MU ) THEN
00645                CALL DLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1),
00646      $                       RWORK(IU2SN+I-1) )
00647             ELSE
00648                CALL DLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1),
00649      $                       RWORK(IU2SN+I-1) )
00650             END IF
00651             RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1)
00652             RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1)
00653 *
00654             TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1)
00655             B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) -
00656      $                  RWORK(IU1SN+I-1)*B11E(I)
00657             B11E(I) = TEMP
00658             IF( I .LT. IMAX - 1 ) THEN
00659                B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1)
00660                B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1)
00661             END IF
00662             TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1)
00663             B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) -
00664      $                  RWORK(IU2SN+I-1)*B21E(I)
00665             B21E(I) = TEMP
00666             IF( I .LT. IMAX - 1 ) THEN
00667                B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1)
00668                B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1)
00669             END IF
00670             TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I)
00671             B12E(I) = RWORK(IU1CS+I-1)*B12E(I) -
00672      $                RWORK(IU1SN+I-1)*B12D(I)
00673             B12D(I) = TEMP
00674             B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1)
00675             B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1)
00676             TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I)
00677             B22E(I) = RWORK(IU2CS+I-1)*B22E(I) -
00678      $                RWORK(IU2SN+I-1)*B22D(I)
00679             B22D(I) = TEMP
00680             B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1)
00681             B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1)
00682 *
00683          END DO
00684 *
00685 *        Compute PHI(IMAX-1)
00686 *
00687          X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) +
00688      $        COS(THETA(IMAX-1))*B21E(IMAX-1)
00689          Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) +
00690      $        COS(THETA(IMAX-1))*B22D(IMAX-1)
00691          Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE
00692 *
00693          PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) )
00694 *
00695 *        Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
00696 *
00697          RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
00698          RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
00699 *
00700          IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
00701             CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1),
00702      $                    RWORK(IV2TCS+IMAX-1-1), R )
00703          ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
00704             CALL DLARTGP( B12BULGE, B12D(IMAX-1),
00705      $                    RWORK(IV2TSN+IMAX-1-1),
00706      $                    RWORK(IV2TCS+IMAX-1-1), R )
00707          ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
00708             CALL DLARTGP( B22BULGE, B22D(IMAX-1),
00709      $                    RWORK(IV2TSN+IMAX-1-1),
00710      $                    RWORK(IV2TCS+IMAX-1-1), R )
00711          ELSE IF( NU .LT. MU ) THEN
00712             CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU,
00713      $                    RWORK(IV2TCS+IMAX-1-1),
00714      $                    RWORK(IV2TSN+IMAX-1-1) )
00715          ELSE
00716             CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU,
00717      $                    RWORK(IV2TCS+IMAX-1-1),
00718      $                    RWORK(IV2TSN+IMAX-1-1) )
00719          END IF
00720 *
00721          TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) +
00722      $          RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX)
00723          B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) -
00724      $                RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1)
00725          B12E(IMAX-1) = TEMP
00726          TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) +
00727      $          RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX)
00728          B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) -
00729      $                RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1)
00730          B22E(IMAX-1) = TEMP
00731 *
00732 *        Update singular vectors
00733 *
00734          IF( WANTU1 ) THEN
00735             IF( COLMAJOR ) THEN
00736                CALL ZLASR( 'R', 'V', 'F', P, IMAX-IMIN+1,
00737      $                     RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
00738      $                     U1(1,IMIN), LDU1 )
00739             ELSE
00740                CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, P,
00741      $                     RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
00742      $                     U1(IMIN,1), LDU1 )
00743             END IF
00744          END IF
00745          IF( WANTU2 ) THEN
00746             IF( COLMAJOR ) THEN
00747                CALL ZLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1,
00748      $                     RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
00749      $                     U2(1,IMIN), LDU2 )
00750             ELSE
00751                CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P,
00752      $                     RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
00753      $                     U2(IMIN,1), LDU2 )
00754             END IF
00755          END IF
00756          IF( WANTV1T ) THEN
00757             IF( COLMAJOR ) THEN
00758                CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q,
00759      $                     RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
00760      $                     V1T(IMIN,1), LDV1T )
00761             ELSE
00762                CALL ZLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1,
00763      $                     RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
00764      $                     V1T(1,IMIN), LDV1T )
00765             END IF
00766          END IF
00767          IF( WANTV2T ) THEN
00768             IF( COLMAJOR ) THEN
00769                CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q,
00770      $                     RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
00771      $                     V2T(IMIN,1), LDV2T )
00772             ELSE
00773                CALL ZLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1,
00774      $                     RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
00775      $                     V2T(1,IMIN), LDV2T )
00776             END IF
00777          END IF
00778 *
00779 *        Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
00780 *
00781          IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN
00782             B11D(IMAX) = -B11D(IMAX)
00783             B21D(IMAX) = -B21D(IMAX)
00784             IF( WANTV1T ) THEN
00785                IF( COLMAJOR ) THEN
00786                   CALL ZSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
00787                ELSE
00788                   CALL ZSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
00789                END IF
00790             END IF
00791          END IF
00792 *
00793 *        Compute THETA(IMAX)
00794 *
00795          X1 = COS(PHI(IMAX-1))*B11D(IMAX) +
00796      $        SIN(PHI(IMAX-1))*B12E(IMAX-1)
00797          Y1 = COS(PHI(IMAX-1))*B21D(IMAX) +
00798      $        SIN(PHI(IMAX-1))*B22E(IMAX-1)
00799 *
00800          THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) )
00801 *
00802 *        Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
00803 *        and B22(IMAX,IMAX-1)
00804 *
00805          IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN
00806             B12D(IMAX) = -B12D(IMAX)
00807             IF( WANTU1 ) THEN
00808                IF( COLMAJOR ) THEN
00809                   CALL ZSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
00810                ELSE
00811                   CALL ZSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
00812                END IF
00813             END IF
00814          END IF
00815          IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN
00816             B22D(IMAX) = -B22D(IMAX)
00817             IF( WANTU2 ) THEN
00818                IF( COLMAJOR ) THEN
00819                   CALL ZSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
00820                ELSE
00821                   CALL ZSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
00822                END IF
00823             END IF
00824          END IF
00825 *
00826 *        Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
00827 *
00828          IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
00829             IF( WANTV2T ) THEN
00830                IF( COLMAJOR ) THEN
00831                   CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
00832                ELSE
00833                   CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
00834                END IF
00835             END IF
00836          END IF
00837 *
00838 *        Test for negligible sines or cosines
00839 *
00840          DO I = IMIN, IMAX
00841             IF( THETA(I) .LT. THRESH ) THEN
00842                THETA(I) = ZERO
00843             ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
00844                THETA(I) = PIOVER2
00845             END IF
00846          END DO
00847          DO I = IMIN, IMAX-1
00848             IF( PHI(I) .LT. THRESH ) THEN
00849                PHI(I) = ZERO
00850             ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
00851                PHI(I) = PIOVER2
00852             END IF
00853          END DO
00854 *
00855 *        Deflate
00856 *
00857          IF (IMAX .GT. 1) THEN
00858             DO WHILE( PHI(IMAX-1) .EQ. ZERO )
00859                IMAX = IMAX - 1
00860                IF (IMAX .LE. 1) EXIT
00861             END DO
00862          END IF
00863          IF( IMIN .GT. IMAX - 1 )
00864      $      IMIN = IMAX - 1
00865          IF (IMIN .GT. 1) THEN
00866             DO WHILE (PHI(IMIN-1) .NE. ZERO)
00867                 IMIN = IMIN - 1
00868                 IF (IMIN .LE. 1) EXIT
00869             END DO
00870          END IF
00871 *
00872 *        Repeat main iteration loop
00873 *
00874       END DO
00875 *
00876 *     Postprocessing: order THETA from least to greatest
00877 *
00878       DO I = 1, Q
00879 *
00880          MINI = I
00881          THETAMIN = THETA(I)
00882          DO J = I+1, Q
00883             IF( THETA(J) .LT. THETAMIN ) THEN
00884                MINI = J
00885                THETAMIN = THETA(J)
00886             END IF
00887          END DO
00888 *
00889          IF( MINI .NE. I ) THEN
00890             THETA(MINI) = THETA(I)
00891             THETA(I) = THETAMIN
00892             IF( COLMAJOR ) THEN
00893                IF( WANTU1 )
00894      $            CALL ZSWAP( P, U1(1,I), 1, U1(1,MINI), 1 )
00895                IF( WANTU2 )
00896      $            CALL ZSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 )
00897                IF( WANTV1T )
00898      $            CALL ZSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T )
00899                IF( WANTV2T )
00900      $            CALL ZSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1),
00901      $               LDV2T )
00902             ELSE
00903                IF( WANTU1 )
00904      $            CALL ZSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 )
00905                IF( WANTU2 )
00906      $            CALL ZSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 )
00907                IF( WANTV1T )
00908      $            CALL ZSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 )
00909                IF( WANTV2T )
00910      $            CALL ZSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 )
00911             END IF
00912          END IF
00913 *
00914       END DO
00915 *
00916       RETURN
00917 *
00918 *     End of ZBBCSD
00919 *
00920       END
00921