"Fossies" - the Fresh Open Source Software Archive

Member "lapack-3.9.1/SRC/cbbcsd.f" (25 Mar 2021, 38992 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 "cbbcsd.f": 3.9.0_vs_3.9.1.

    1 *> \brief \b CBBCSD
    2 *
    3 *  =========== DOCUMENTATION ===========
    4 *
    5 * Online html documentation available at
    6 *            http://www.netlib.org/lapack/explore-html/
    7 *
    8 *> \htmlonly
    9 *> Download CBBCSD + dependencies
   10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cbbcsd.f">
   11 *> [TGZ]</a>
   12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cbbcsd.f">
   13 *> [ZIP]</a>
   14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cbbcsd.f">
   15 *> [TXT]</a>
   16 *> \endhtmlonly
   17 *
   18 *  Definition:
   19 *  ===========
   20 *
   21 *       SUBROUTINE CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
   22 *                          THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
   23 *                          V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
   24 *                          B22D, B22E, RWORK, LRWORK, INFO )
   25 *
   26 *       .. Scalar Arguments ..
   27 *       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
   28 *       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
   29 *       ..
   30 *       .. Array Arguments ..
   31 *       REAL               B11D( * ), B11E( * ), B12D( * ), B12E( * ),
   32 *      $                   B21D( * ), B21E( * ), B22D( * ), B22E( * ),
   33 *      $                   PHI( * ), THETA( * ), RWORK( * )
   34 *       COMPLEX            U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
   35 *      $                   V2T( LDV2T, * )
   36 *       ..
   37 *
   38 *
   39 *> \par Purpose:
   40 *  =============
   41 *>
   42 *> \verbatim
   43 *>
   44 *> CBBCSD computes the CS decomposition of a unitary matrix in
   45 *> bidiagonal-block form,
   46 *>
   47 *>
   48 *>     [ B11 | B12 0  0 ]
   49 *>     [  0  |  0 -I  0 ]
   50 *> X = [----------------]
   51 *>     [ B21 | B22 0  0 ]
   52 *>     [  0  |  0  0  I ]
   53 *>
   54 *>                               [  C | -S  0  0 ]
   55 *>                   [ U1 |    ] [  0 |  0 -I  0 ] [ V1 |    ]**H
   56 *>                 = [---------] [---------------] [---------]   .
   57 *>                   [    | U2 ] [  S |  C  0  0 ] [    | V2 ]
   58 *>                               [  0 |  0  0  I ]
   59 *>
   60 *> X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
   61 *> than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
   62 *> transposed and/or permuted. This can be done in constant time using
   63 *> the TRANS and SIGNS options. See CUNCSD for details.)
   64 *>
   65 *> The bidiagonal matrices B11, B12, B21, and B22 are represented
   66 *> implicitly by angles THETA(1:Q) and PHI(1:Q-1).
   67 *>
   68 *> The unitary matrices U1, U2, V1T, and V2T are input/output.
   69 *> The input matrices are pre- or post-multiplied by the appropriate
   70 *> singular vector matrices.
   71 *> \endverbatim
   72 *
   73 *  Arguments:
   74 *  ==========
   75 *
   76 *> \param[in] JOBU1
   77 *> \verbatim
   78 *>          JOBU1 is CHARACTER
   79 *>          = 'Y':      U1 is updated;
   80 *>          otherwise:  U1 is not updated.
   81 *> \endverbatim
   82 *>
   83 *> \param[in] JOBU2
   84 *> \verbatim
   85 *>          JOBU2 is CHARACTER
   86 *>          = 'Y':      U2 is updated;
   87 *>          otherwise:  U2 is not updated.
   88 *> \endverbatim
   89 *>
   90 *> \param[in] JOBV1T
   91 *> \verbatim
   92 *>          JOBV1T is CHARACTER
   93 *>          = 'Y':      V1T is updated;
   94 *>          otherwise:  V1T is not updated.
   95 *> \endverbatim
   96 *>
   97 *> \param[in] JOBV2T
   98 *> \verbatim
   99 *>          JOBV2T is CHARACTER
  100 *>          = 'Y':      V2T is updated;
  101 *>          otherwise:  V2T is not updated.
  102 *> \endverbatim
  103 *>
  104 *> \param[in] TRANS
  105 *> \verbatim
  106 *>          TRANS is CHARACTER
  107 *>          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
  108 *>                      order;
  109 *>          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
  110 *>                      major order.
  111 *> \endverbatim
  112 *>
  113 *> \param[in] M
  114 *> \verbatim
  115 *>          M is INTEGER
  116 *>          The number of rows and columns in X, the unitary matrix in
  117 *>          bidiagonal-block form.
  118 *> \endverbatim
  119 *>
  120 *> \param[in] P
  121 *> \verbatim
  122 *>          P is INTEGER
  123 *>          The number of rows in the top-left block of X. 0 <= P <= M.
  124 *> \endverbatim
  125 *>
  126 *> \param[in] Q
  127 *> \verbatim
  128 *>          Q is INTEGER
  129 *>          The number of columns in the top-left block of X.
  130 *>          0 <= Q <= MIN(P,M-P,M-Q).
  131 *> \endverbatim
  132 *>
  133 *> \param[in,out] THETA
  134 *> \verbatim
  135 *>          THETA is REAL array, dimension (Q)
  136 *>          On entry, the angles THETA(1),...,THETA(Q) that, along with
  137 *>          PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
  138 *>          form. On exit, the angles whose cosines and sines define the
  139 *>          diagonal blocks in the CS decomposition.
  140 *> \endverbatim
  141 *>
  142 *> \param[in,out] PHI
  143 *> \verbatim
  144 *>          PHI is REAL array, dimension (Q-1)
  145 *>          The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
  146 *>          THETA(Q), define the matrix in bidiagonal-block form.
  147 *> \endverbatim
  148 *>
  149 *> \param[in,out] U1
  150 *> \verbatim
  151 *>          U1 is COMPLEX array, dimension (LDU1,P)
  152 *>          On entry, a P-by-P matrix. On exit, U1 is postmultiplied
  153 *>          by the left singular vector matrix common to [ B11 ; 0 ] and
  154 *>          [ B12 0 0 ; 0 -I 0 0 ].
  155 *> \endverbatim
  156 *>
  157 *> \param[in] LDU1
  158 *> \verbatim
  159 *>          LDU1 is INTEGER
  160 *>          The leading dimension of the array U1, LDU1 >= MAX(1,P).
  161 *> \endverbatim
  162 *>
  163 *> \param[in,out] U2
  164 *> \verbatim
  165 *>          U2 is COMPLEX array, dimension (LDU2,M-P)
  166 *>          On entry, an (M-P)-by-(M-P) matrix. On exit, U2 is
  167 *>          postmultiplied by the left singular vector matrix common to
  168 *>          [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
  169 *> \endverbatim
  170 *>
  171 *> \param[in] LDU2
  172 *> \verbatim
  173 *>          LDU2 is INTEGER
  174 *>          The leading dimension of the array U2, LDU2 >= MAX(1,M-P).
  175 *> \endverbatim
  176 *>
  177 *> \param[in,out] V1T
  178 *> \verbatim
  179 *>          V1T is COMPLEX array, dimension (LDV1T,Q)
  180 *>          On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
  181 *>          by the conjugate transpose of the right singular vector
  182 *>          matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
  183 *> \endverbatim
  184 *>
  185 *> \param[in] LDV1T
  186 *> \verbatim
  187 *>          LDV1T is INTEGER
  188 *>          The leading dimension of the array V1T, LDV1T >= MAX(1,Q).
  189 *> \endverbatim
  190 *>
  191 *> \param[in,out] V2T
  192 *> \verbatim
  193 *>          V2T is COMPLEX array, dimension (LDV2T,M-Q)
  194 *>          On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
  195 *>          premultiplied by the conjugate transpose of the right
  196 *>          singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
  197 *>          [ B22 0 0 ; 0 0 I ].
  198 *> \endverbatim
  199 *>
  200 *> \param[in] LDV2T
  201 *> \verbatim
  202 *>          LDV2T is INTEGER
  203 *>          The leading dimension of the array V2T, LDV2T >= MAX(1,M-Q).
  204 *> \endverbatim
  205 *>
  206 *> \param[out] B11D
  207 *> \verbatim
  208 *>          B11D is REAL array, dimension (Q)
  209 *>          When CBBCSD converges, B11D contains the cosines of THETA(1),
  210 *>          ..., THETA(Q). If CBBCSD fails to converge, then B11D
  211 *>          contains the diagonal of the partially reduced top-left
  212 *>          block.
  213 *> \endverbatim
  214 *>
  215 *> \param[out] B11E
  216 *> \verbatim
  217 *>          B11E is REAL array, dimension (Q-1)
  218 *>          When CBBCSD converges, B11E contains zeros. If CBBCSD fails
  219 *>          to converge, then B11E contains the superdiagonal of the
  220 *>          partially reduced top-left block.
  221 *> \endverbatim
  222 *>
  223 *> \param[out] B12D
  224 *> \verbatim
  225 *>          B12D is REAL array, dimension (Q)
  226 *>          When CBBCSD converges, B12D contains the negative sines of
  227 *>          THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then
  228 *>          B12D contains the diagonal of the partially reduced top-right
  229 *>          block.
  230 *> \endverbatim
  231 *>
  232 *> \param[out] B12E
  233 *> \verbatim
  234 *>          B12E is REAL array, dimension (Q-1)
  235 *>          When CBBCSD converges, B12E contains zeros. If CBBCSD fails
  236 *>          to converge, then B12E contains the subdiagonal of the
  237 *>          partially reduced top-right block.
  238 *> \endverbatim
  239 *>
  240 *> \param[out] B21D
  241 *> \verbatim
  242 *>          B21D is REAL array, dimension (Q)
  243 *>          When CBBCSD converges, B21D contains the negative sines of
  244 *>          THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then
  245 *>          B21D contains the diagonal of the partially reduced bottom-left
  246 *>          block.
  247 *> \endverbatim
  248 *>
  249 *> \param[out] B21E
  250 *> \verbatim
  251 *>          B21E is REAL array, dimension (Q-1)
  252 *>          When CBBCSD converges, B21E contains zeros. If CBBCSD fails
  253 *>          to converge, then B21E contains the subdiagonal of the
  254 *>          partially reduced bottom-left block.
  255 *> \endverbatim
  256 *>
  257 *> \param[out] B22D
  258 *> \verbatim
  259 *>          B22D is REAL array, dimension (Q)
  260 *>          When CBBCSD converges, B22D contains the negative sines of
  261 *>          THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then
  262 *>          B22D contains the diagonal of the partially reduced bottom-right
  263 *>          block.
  264 *> \endverbatim
  265 *>
  266 *> \param[out] B22E
  267 *> \verbatim
  268 *>          B22E is REAL array, dimension (Q-1)
  269 *>          When CBBCSD converges, B22E contains zeros. If CBBCSD fails
  270 *>          to converge, then B22E contains the subdiagonal of the
  271 *>          partially reduced bottom-right block.
  272 *> \endverbatim
  273 *>
  274 *> \param[out] RWORK
  275 *> \verbatim
  276 *>          RWORK is REAL array, dimension (MAX(1,LRWORK))
  277 *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
  278 *> \endverbatim
  279 *>
  280 *> \param[in] LRWORK
  281 *> \verbatim
  282 *>          LRWORK is INTEGER
  283 *>          The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).
  284 *>
  285 *>          If LRWORK = -1, then a workspace query is assumed; the
  286 *>          routine only calculates the optimal size of the RWORK array,
  287 *>          returns this value as the first entry of the work array, and
  288 *>          no error message related to LRWORK is issued by XERBLA.
  289 *> \endverbatim
  290 *>
  291 *> \param[out] INFO
  292 *> \verbatim
  293 *>          INFO is INTEGER
  294 *>          = 0:  successful exit.
  295 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  296 *>          > 0:  if CBBCSD did not converge, INFO specifies the number
  297 *>                of nonzero entries in PHI, and B11D, B11E, etc.,
  298 *>                contain the partially reduced matrix.
  299 *> \endverbatim
  300 *
  301 *> \par Internal Parameters:
  302 *  =========================
  303 *>
  304 *> \verbatim
  305 *>  TOLMUL  REAL, default = MAX(10,MIN(100,EPS**(-1/8)))
  306 *>          TOLMUL controls the convergence criterion of the QR loop.
  307 *>          Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
  308 *>          are within TOLMUL*EPS of either bound.
  309 *> \endverbatim
  310 *
  311 *> \par References:
  312 *  ================
  313 *>
  314 *>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
  315 *>      Algorithms, 50(1):33-65, 2009.
  316 *
  317 *  Authors:
  318 *  ========
  319 *
  320 *> \author Univ. of Tennessee
  321 *> \author Univ. of California Berkeley
  322 *> \author Univ. of Colorado Denver
  323 *> \author NAG Ltd.
  324 *
  325 *> \ingroup complexOTHERcomputational
  326 *
  327 *  =====================================================================
  328       SUBROUTINE CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
  329      $                   THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
  330      $                   V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
  331      $                   B22D, B22E, RWORK, LRWORK, INFO )
  332 *
  333 *  -- LAPACK computational routine --
  334 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  335 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  336 *
  337 *     .. Scalar Arguments ..
  338       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
  339       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
  340 *     ..
  341 *     .. Array Arguments ..
  342       REAL               B11D( * ), B11E( * ), B12D( * ), B12E( * ),
  343      $                   B21D( * ), B21E( * ), B22D( * ), B22E( * ),
  344      $                   PHI( * ), THETA( * ), RWORK( * )
  345       COMPLEX            U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
  346      $                   V2T( LDV2T, * )
  347 *     ..
  348 *
  349 *  ===================================================================
  350 *
  351 *     .. Parameters ..
  352       INTEGER            MAXITR
  353       PARAMETER          ( MAXITR = 6 )
  354       REAL               HUNDRED, MEIGHTH, ONE, TEN, ZERO
  355       PARAMETER          ( HUNDRED = 100.0E0, MEIGHTH = -0.125E0,
  356      $                     ONE = 1.0E0, TEN = 10.0E0, ZERO = 0.0E0 )
  357       COMPLEX            NEGONECOMPLEX
  358       PARAMETER          ( NEGONECOMPLEX = (-1.0E0,0.0E0) )
  359       REAL               PIOVER2
  360       PARAMETER ( PIOVER2 = 1.57079632679489661923132169163975144210E0 )
  361 *     ..
  362 *     .. Local Scalars ..
  363       LOGICAL            COLMAJOR, LQUERY, RESTART11, RESTART12,
  364      $                   RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
  365      $                   WANTV2T
  366       INTEGER            I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
  367      $                   IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
  368      $                   LRWORKMIN, LRWORKOPT, MAXIT, MINI
  369       REAL               B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
  370      $                   EPS, MU, NU, R, SIGMA11, SIGMA21,
  371      $                   TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
  372      $                   UNFL, X1, X2, Y1, Y2
  373 *
  374 *     .. External Subroutines ..
  375       EXTERNAL           CLASR, CSCAL, CSWAP, SLARTGP, SLARTGS, SLAS2,
  376      $                   XERBLA
  377 *     ..
  378 *     .. External Functions ..
  379       REAL               SLAMCH
  380       LOGICAL            LSAME
  381       EXTERNAL           LSAME, SLAMCH
  382 *     ..
  383 *     .. Intrinsic Functions ..
  384       INTRINSIC          ABS, ATAN2, COS, MAX, MIN, SIN, SQRT
  385 *     ..
  386 *     .. Executable Statements ..
  387 *
  388 *     Test input arguments
  389 *
  390       INFO = 0
  391       LQUERY = LRWORK .EQ. -1
  392       WANTU1 = LSAME( JOBU1, 'Y' )
  393       WANTU2 = LSAME( JOBU2, 'Y' )
  394       WANTV1T = LSAME( JOBV1T, 'Y' )
  395       WANTV2T = LSAME( JOBV2T, 'Y' )
  396       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
  397 *
  398       IF( M .LT. 0 ) THEN
  399          INFO = -6
  400       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
  401          INFO = -7
  402       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
  403          INFO = -8
  404       ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN
  405          INFO = -8
  406       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
  407          INFO = -12
  408       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
  409          INFO = -14
  410       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
  411          INFO = -16
  412       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
  413          INFO = -18
  414       END IF
  415 *
  416 *     Quick return if Q = 0
  417 *
  418       IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN
  419          LRWORKMIN = 1
  420          RWORK(1) = LRWORKMIN
  421          RETURN
  422       END IF
  423 *
  424 *     Compute workspace
  425 *
  426       IF( INFO .EQ. 0 ) THEN
  427          IU1CS = 1
  428          IU1SN = IU1CS + Q
  429          IU2CS = IU1SN + Q
  430          IU2SN = IU2CS + Q
  431          IV1TCS = IU2SN + Q
  432          IV1TSN = IV1TCS + Q
  433          IV2TCS = IV1TSN + Q
  434          IV2TSN = IV2TCS + Q
  435          LRWORKOPT = IV2TSN + Q - 1
  436          LRWORKMIN = LRWORKOPT
  437          RWORK(1) = LRWORKOPT
  438          IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN
  439             INFO = -28
  440          END IF
  441       END IF
  442 *
  443       IF( INFO .NE. 0 ) THEN
  444          CALL XERBLA( 'CBBCSD', -INFO )
  445          RETURN
  446       ELSE IF( LQUERY ) THEN
  447          RETURN
  448       END IF
  449 *
  450 *     Get machine constants
  451 *
  452       EPS = SLAMCH( 'Epsilon' )
  453       UNFL = SLAMCH( 'Safe minimum' )
  454       TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) )
  455       TOL = TOLMUL*EPS
  456       THRESH = MAX( TOL, MAXITR*Q*Q*UNFL )
  457 *
  458 *     Test for negligible sines or cosines
  459 *
  460       DO I = 1, Q
  461          IF( THETA(I) .LT. THRESH ) THEN
  462             THETA(I) = ZERO
  463          ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
  464             THETA(I) = PIOVER2
  465          END IF
  466       END DO
  467       DO I = 1, Q-1
  468          IF( PHI(I) .LT. THRESH ) THEN
  469             PHI(I) = ZERO
  470          ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
  471             PHI(I) = PIOVER2
  472          END IF
  473       END DO
  474 *
  475 *     Initial deflation
  476 *
  477       IMAX = Q
  478       DO WHILE( IMAX .GT. 1 )
  479          IF( PHI(IMAX-1) .NE. ZERO ) THEN
  480             EXIT
  481          END IF
  482          IMAX = IMAX - 1
  483       END DO
  484       IMIN = IMAX - 1
  485       IF  ( IMIN .GT. 1 ) THEN
  486          DO WHILE( PHI(IMIN-1) .NE. ZERO )
  487             IMIN = IMIN - 1
  488             IF  ( IMIN .LE. 1 ) EXIT
  489          END DO
  490       END IF
  491 *
  492 *     Initialize iteration counter
  493 *
  494       MAXIT = MAXITR*Q*Q
  495       ITER = 0
  496 *
  497 *     Begin main iteration loop
  498 *
  499       DO WHILE( IMAX .GT. 1 )
  500 *
  501 *        Compute the matrix entries
  502 *
  503          B11D(IMIN) = COS( THETA(IMIN) )
  504          B21D(IMIN) = -SIN( THETA(IMIN) )
  505          DO I = IMIN, IMAX - 1
  506             B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) )
  507             B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) )
  508             B12D(I) = SIN( THETA(I) ) * COS( PHI(I) )
  509             B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) )
  510             B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) )
  511             B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) )
  512             B22D(I) = COS( THETA(I) ) * COS( PHI(I) )
  513             B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) )
  514          END DO
  515          B12D(IMAX) = SIN( THETA(IMAX) )
  516          B22D(IMAX) = COS( THETA(IMAX) )
  517 *
  518 *        Abort if not converging; otherwise, increment ITER
  519 *
  520          IF( ITER .GT. MAXIT ) THEN
  521             INFO = 0
  522             DO I = 1, Q
  523                IF( PHI(I) .NE. ZERO )
  524      $            INFO = INFO + 1
  525             END DO
  526             RETURN
  527          END IF
  528 *
  529          ITER = ITER + IMAX - IMIN
  530 *
  531 *        Compute shifts
  532 *
  533          THETAMAX = THETA(IMIN)
  534          THETAMIN = THETA(IMIN)
  535          DO I = IMIN+1, IMAX
  536             IF( THETA(I) > THETAMAX )
  537      $         THETAMAX = THETA(I)
  538             IF( THETA(I) < THETAMIN )
  539      $         THETAMIN = THETA(I)
  540          END DO
  541 *
  542          IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN
  543 *
  544 *           Zero on diagonals of B11 and B22; induce deflation with a
  545 *           zero shift
  546 *
  547             MU = ZERO
  548             NU = ONE
  549 *
  550          ELSE IF( THETAMIN .LT. THRESH ) THEN
  551 *
  552 *           Zero on diagonals of B12 and B22; induce deflation with a
  553 *           zero shift
  554 *
  555             MU = ONE
  556             NU = ZERO
  557 *
  558          ELSE
  559 *
  560 *           Compute shifts for B11 and B21 and use the lesser
  561 *
  562             CALL SLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11,
  563      $                  DUMMY )
  564             CALL SLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21,
  565      $                  DUMMY )
  566 *
  567             IF( SIGMA11 .LE. SIGMA21 ) THEN
  568                MU = SIGMA11
  569                NU = SQRT( ONE - MU**2 )
  570                IF( MU .LT. THRESH ) THEN
  571                   MU = ZERO
  572                   NU = ONE
  573                END IF
  574             ELSE
  575                NU = SIGMA21
  576                MU = SQRT( 1.0 - NU**2 )
  577                IF( NU .LT. THRESH ) THEN
  578                   MU = ONE
  579                   NU = ZERO
  580                END IF
  581             END IF
  582          END IF
  583 *
  584 *        Rotate to produce bulges in B11 and B21
  585 *
  586          IF( MU .LE. NU ) THEN
  587             CALL SLARTGS( B11D(IMIN), B11E(IMIN), MU,
  588      $                    RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
  589          ELSE
  590             CALL SLARTGS( B21D(IMIN), B21E(IMIN), NU,
  591      $                    RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
  592          END IF
  593 *
  594          TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) +
  595      $          RWORK(IV1TSN+IMIN-1)*B11E(IMIN)
  596          B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) -
  597      $                RWORK(IV1TSN+IMIN-1)*B11D(IMIN)
  598          B11D(IMIN) = TEMP
  599          B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1)
  600          B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1)
  601          TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) +
  602      $          RWORK(IV1TSN+IMIN-1)*B21E(IMIN)
  603          B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) -
  604      $                RWORK(IV1TSN+IMIN-1)*B21D(IMIN)
  605          B21D(IMIN) = TEMP
  606          B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1)
  607          B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1)
  608 *
  609 *        Compute THETA(IMIN)
  610 *
  611          THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ),
  612      $                   SQRT( B11D(IMIN)**2+B11BULGE**2 ) )
  613 *
  614 *        Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
  615 *
  616          IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
  617             CALL SLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1),
  618      $                    RWORK(IU1CS+IMIN-1), R )
  619          ELSE IF( MU .LE. NU ) THEN
  620             CALL SLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU,
  621      $                    RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
  622          ELSE
  623             CALL SLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
  624      $                    RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
  625          END IF
  626          IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
  627             CALL SLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1),
  628      $                    RWORK(IU2CS+IMIN-1), R )
  629          ELSE IF( NU .LT. MU ) THEN
  630             CALL SLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU,
  631      $                    RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
  632          ELSE
  633             CALL SLARTGS( B22D(IMIN), B22E(IMIN), MU,
  634      $                    RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
  635          END IF
  636          RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1)
  637          RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1)
  638 *
  639          TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) +
  640      $          RWORK(IU1SN+IMIN-1)*B11D(IMIN+1)
  641          B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) -
  642      $                  RWORK(IU1SN+IMIN-1)*B11E(IMIN)
  643          B11E(IMIN) = TEMP
  644          IF( IMAX .GT. IMIN+1 ) THEN
  645             B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1)
  646             B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1)
  647          END IF
  648          TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) +
  649      $          RWORK(IU1SN+IMIN-1)*B12E(IMIN)
  650          B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) -
  651      $                RWORK(IU1SN+IMIN-1)*B12D(IMIN)
  652          B12D(IMIN) = TEMP
  653          B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1)
  654          B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1)
  655          TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) +
  656      $          RWORK(IU2SN+IMIN-1)*B21D(IMIN+1)
  657          B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) -
  658      $                  RWORK(IU2SN+IMIN-1)*B21E(IMIN)
  659          B21E(IMIN) = TEMP
  660          IF( IMAX .GT. IMIN+1 ) THEN
  661             B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1)
  662             B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1)
  663          END IF
  664          TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) +
  665      $          RWORK(IU2SN+IMIN-1)*B22E(IMIN)
  666          B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) -
  667      $                RWORK(IU2SN+IMIN-1)*B22D(IMIN)
  668          B22D(IMIN) = TEMP
  669          B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1)
  670          B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1)
  671 *
  672 *        Inner loop: chase bulges from B11(IMIN,IMIN+2),
  673 *        B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
  674 *        bottom-right
  675 *
  676          DO I = IMIN+1, IMAX-1
  677 *
  678 *           Compute PHI(I-1)
  679 *
  680             X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1)
  681             X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE
  682             Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1)
  683             Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE
  684 *
  685             PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) )
  686 *
  687 *           Determine if there are bulges to chase or if a new direct
  688 *           summand has been reached
  689 *
  690             RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
  691             RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
  692             RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
  693             RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
  694 *
  695 *           If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
  696 *           B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
  697 *           chasing by applying the original shift again.
  698 *
  699             IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN
  700                CALL SLARTGP( X2, X1, RWORK(IV1TSN+I-1),
  701      $                       RWORK(IV1TCS+I-1), R )
  702             ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN
  703                CALL SLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1),
  704      $                       RWORK(IV1TCS+I-1), R )
  705             ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN
  706                CALL SLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1),
  707      $                       RWORK(IV1TCS+I-1), R )
  708             ELSE IF( MU .LE. NU ) THEN
  709                CALL SLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1),
  710      $                       RWORK(IV1TSN+I-1) )
  711             ELSE
  712                CALL SLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1),
  713      $                       RWORK(IV1TSN+I-1) )
  714             END IF
  715             RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1)
  716             RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1)
  717             IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
  718                CALL SLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1),
  719      $                       RWORK(IV2TCS+I-1-1), R )
  720             ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
  721                CALL SLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1),
  722      $                       RWORK(IV2TCS+I-1-1), R )
  723             ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
  724                CALL SLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1),
  725      $                       RWORK(IV2TCS+I-1-1), R )
  726             ELSE IF( NU .LT. MU ) THEN
  727                CALL SLARTGS( B12E(I-1), B12D(I), NU,
  728      $                       RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
  729             ELSE
  730                CALL SLARTGS( B22E(I-1), B22D(I), MU,
  731      $                       RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
  732             END IF
  733 *
  734             TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I)
  735             B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) -
  736      $                RWORK(IV1TSN+I-1)*B11D(I)
  737             B11D(I) = TEMP
  738             B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1)
  739             B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1)
  740             TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I)
  741             B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) -
  742      $                RWORK(IV1TSN+I-1)*B21D(I)
  743             B21D(I) = TEMP
  744             B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1)
  745             B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1)
  746             TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) +
  747      $             RWORK(IV2TSN+I-1-1)*B12D(I)
  748             B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) -
  749      $                RWORK(IV2TSN+I-1-1)*B12E(I-1)
  750             B12E(I-1) = TEMP
  751             B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I)
  752             B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I)
  753             TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) +
  754      $             RWORK(IV2TSN+I-1-1)*B22D(I)
  755             B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) -
  756      $                RWORK(IV2TSN+I-1-1)*B22E(I-1)
  757             B22E(I-1) = TEMP
  758             B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I)
  759             B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I)
  760 *
  761 *           Compute THETA(I)
  762 *
  763             X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1)
  764             X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE
  765             Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1)
  766             Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE
  767 *
  768             THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) )
  769 *
  770 *           Determine if there are bulges to chase or if a new direct
  771 *           summand has been reached
  772 *
  773             RESTART11 =   B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
  774             RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
  775             RESTART21 =   B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
  776             RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
  777 *
  778 *           If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
  779 *           B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
  780 *           chasing by applying the original shift again.
  781 *
  782             IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN
  783                CALL SLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1),
  784      $                       R )
  785             ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN
  786                CALL SLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1),
  787      $                       RWORK(IU1CS+I-1), R )
  788             ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN
  789                CALL SLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1),
  790      $                       RWORK(IU1CS+I-1), R )
  791             ELSE IF( MU .LE. NU ) THEN
  792                CALL SLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1),
  793      $                       RWORK(IU1SN+I-1) )
  794             ELSE
  795                CALL SLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1),
  796      $                       RWORK(IU1SN+I-1) )
  797             END IF
  798             IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN
  799                CALL SLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1),
  800      $                       R )
  801             ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN
  802                CALL SLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1),
  803      $                       RWORK(IU2CS+I-1), R )
  804             ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN
  805                CALL SLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1),
  806      $                       RWORK(IU2CS+I-1), R )
  807             ELSE IF( NU .LT. MU ) THEN
  808                CALL SLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1),
  809      $                       RWORK(IU2SN+I-1) )
  810             ELSE
  811                CALL SLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1),
  812      $                       RWORK(IU2SN+I-1) )
  813             END IF
  814             RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1)
  815             RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1)
  816 *
  817             TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1)
  818             B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) -
  819      $                  RWORK(IU1SN+I-1)*B11E(I)
  820             B11E(I) = TEMP
  821             IF( I .LT. IMAX - 1 ) THEN
  822                B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1)
  823                B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1)
  824             END IF
  825             TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1)
  826             B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) -
  827      $                  RWORK(IU2SN+I-1)*B21E(I)
  828             B21E(I) = TEMP
  829             IF( I .LT. IMAX - 1 ) THEN
  830                B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1)
  831                B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1)
  832             END IF
  833             TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I)
  834             B12E(I) = RWORK(IU1CS+I-1)*B12E(I) -
  835      $                RWORK(IU1SN+I-1)*B12D(I)
  836             B12D(I) = TEMP
  837             B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1)
  838             B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1)
  839             TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I)
  840             B22E(I) = RWORK(IU2CS+I-1)*B22E(I) -
  841      $                RWORK(IU2SN+I-1)*B22D(I)
  842             B22D(I) = TEMP
  843             B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1)
  844             B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1)
  845 *
  846          END DO
  847 *
  848 *        Compute PHI(IMAX-1)
  849 *
  850          X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) +
  851      $        COS(THETA(IMAX-1))*B21E(IMAX-1)
  852          Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) +
  853      $        COS(THETA(IMAX-1))*B22D(IMAX-1)
  854          Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE
  855 *
  856          PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) )
  857 *
  858 *        Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
  859 *
  860          RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
  861          RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
  862 *
  863          IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
  864             CALL SLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1),
  865      $                    RWORK(IV2TCS+IMAX-1-1), R )
  866          ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
  867             CALL SLARTGP( B12BULGE, B12D(IMAX-1),
  868      $                    RWORK(IV2TSN+IMAX-1-1),
  869      $                    RWORK(IV2TCS+IMAX-1-1), R )
  870          ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
  871             CALL SLARTGP( B22BULGE, B22D(IMAX-1),
  872      $                    RWORK(IV2TSN+IMAX-1-1),
  873      $                    RWORK(IV2TCS+IMAX-1-1), R )
  874          ELSE IF( NU .LT. MU ) THEN
  875             CALL SLARTGS( B12E(IMAX-1), B12D(IMAX), NU,
  876      $                    RWORK(IV2TCS+IMAX-1-1),
  877      $                    RWORK(IV2TSN+IMAX-1-1) )
  878          ELSE
  879             CALL SLARTGS( B22E(IMAX-1), B22D(IMAX), MU,
  880      $                    RWORK(IV2TCS+IMAX-1-1),
  881      $                    RWORK(IV2TSN+IMAX-1-1) )
  882          END IF
  883 *
  884          TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) +
  885      $          RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX)
  886          B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) -
  887      $                RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1)
  888          B12E(IMAX-1) = TEMP
  889          TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) +
  890      $          RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX)
  891          B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) -
  892      $                RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1)
  893          B22E(IMAX-1) = TEMP
  894 *
  895 *        Update singular vectors
  896 *
  897          IF( WANTU1 ) THEN
  898             IF( COLMAJOR ) THEN
  899                CALL CLASR( 'R', 'V', 'F', P, IMAX-IMIN+1,
  900      $                     RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
  901      $                     U1(1,IMIN), LDU1 )
  902             ELSE
  903                CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, P,
  904      $                     RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
  905      $                     U1(IMIN,1), LDU1 )
  906             END IF
  907          END IF
  908          IF( WANTU2 ) THEN
  909             IF( COLMAJOR ) THEN
  910                CALL CLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1,
  911      $                     RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
  912      $                     U2(1,IMIN), LDU2 )
  913             ELSE
  914                CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P,
  915      $                     RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
  916      $                     U2(IMIN,1), LDU2 )
  917             END IF
  918          END IF
  919          IF( WANTV1T ) THEN
  920             IF( COLMAJOR ) THEN
  921                CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q,
  922      $                     RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
  923      $                     V1T(IMIN,1), LDV1T )
  924             ELSE
  925                CALL CLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1,
  926      $                     RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
  927      $                     V1T(1,IMIN), LDV1T )
  928             END IF
  929          END IF
  930          IF( WANTV2T ) THEN
  931             IF( COLMAJOR ) THEN
  932                CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q,
  933      $                     RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
  934      $                     V2T(IMIN,1), LDV2T )
  935             ELSE
  936                CALL CLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1,
  937      $                     RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
  938      $                     V2T(1,IMIN), LDV2T )
  939             END IF
  940          END IF
  941 *
  942 *        Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
  943 *
  944          IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN
  945             B11D(IMAX) = -B11D(IMAX)
  946             B21D(IMAX) = -B21D(IMAX)
  947             IF( WANTV1T ) THEN
  948                IF( COLMAJOR ) THEN
  949                   CALL CSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
  950                ELSE
  951                   CALL CSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
  952                END IF
  953             END IF
  954          END IF
  955 *
  956 *        Compute THETA(IMAX)
  957 *
  958          X1 = COS(PHI(IMAX-1))*B11D(IMAX) +
  959      $        SIN(PHI(IMAX-1))*B12E(IMAX-1)
  960          Y1 = COS(PHI(IMAX-1))*B21D(IMAX) +
  961      $        SIN(PHI(IMAX-1))*B22E(IMAX-1)
  962 *
  963          THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) )
  964 *
  965 *        Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
  966 *        and B22(IMAX,IMAX-1)
  967 *
  968          IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN
  969             B12D(IMAX) = -B12D(IMAX)
  970             IF( WANTU1 ) THEN
  971                IF( COLMAJOR ) THEN
  972                   CALL CSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
  973                ELSE
  974                   CALL CSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
  975                END IF
  976             END IF
  977          END IF
  978          IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN
  979             B22D(IMAX) = -B22D(IMAX)
  980             IF( WANTU2 ) THEN
  981                IF( COLMAJOR ) THEN
  982                   CALL CSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
  983                ELSE
  984                   CALL CSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
  985                END IF
  986             END IF
  987          END IF
  988 *
  989 *        Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
  990 *
  991          IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
  992             IF( WANTV2T ) THEN
  993                IF( COLMAJOR ) THEN
  994                   CALL CSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
  995                ELSE
  996                   CALL CSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
  997                END IF
  998             END IF
  999          END IF
 1000 *
 1001 *        Test for negligible sines or cosines
 1002 *
 1003          DO I = IMIN, IMAX
 1004             IF( THETA(I) .LT. THRESH ) THEN
 1005                THETA(I) = ZERO
 1006             ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
 1007                THETA(I) = PIOVER2
 1008             END IF
 1009          END DO
 1010          DO I = IMIN, IMAX-1
 1011             IF( PHI(I) .LT. THRESH ) THEN
 1012                PHI(I) = ZERO
 1013             ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
 1014                PHI(I) = PIOVER2
 1015             END IF
 1016          END DO
 1017 *
 1018 *        Deflate
 1019 *
 1020          IF (IMAX .GT. 1) THEN
 1021             DO WHILE( PHI(IMAX-1) .EQ. ZERO )
 1022                IMAX = IMAX - 1
 1023                IF (IMAX .LE. 1) EXIT
 1024             END DO
 1025          END IF
 1026          IF( IMIN .GT. IMAX - 1 )
 1027      $      IMIN = IMAX - 1
 1028          IF (IMIN .GT. 1) THEN
 1029             DO WHILE (PHI(IMIN-1) .NE. ZERO)
 1030                 IMIN = IMIN - 1
 1031                 IF (IMIN .LE. 1) EXIT
 1032             END DO
 1033          END IF
 1034 *
 1035 *        Repeat main iteration loop
 1036 *
 1037       END DO
 1038 *
 1039 *     Postprocessing: order THETA from least to greatest
 1040 *
 1041       DO I = 1, Q
 1042 *
 1043          MINI = I
 1044          THETAMIN = THETA(I)
 1045          DO J = I+1, Q
 1046             IF( THETA(J) .LT. THETAMIN ) THEN
 1047                MINI = J
 1048                THETAMIN = THETA(J)
 1049             END IF
 1050          END DO
 1051 *
 1052          IF( MINI .NE. I ) THEN
 1053             THETA(MINI) = THETA(I)
 1054             THETA(I) = THETAMIN
 1055             IF( COLMAJOR ) THEN
 1056                IF( WANTU1 )
 1057      $            CALL CSWAP( P, U1(1,I), 1, U1(1,MINI), 1 )
 1058                IF( WANTU2 )
 1059      $            CALL CSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 )
 1060                IF( WANTV1T )
 1061      $            CALL CSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T )
 1062                IF( WANTV2T )
 1063      $            CALL CSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1),
 1064      $               LDV2T )
 1065             ELSE
 1066                IF( WANTU1 )
 1067      $            CALL CSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 )
 1068                IF( WANTU2 )
 1069      $            CALL CSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 )
 1070                IF( WANTV1T )
 1071      $            CALL CSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 )
 1072                IF( WANTV2T )
 1073      $            CALL CSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 )
 1074             END IF
 1075          END IF
 1076 *
 1077       END DO
 1078 *
 1079       RETURN
 1080 *
 1081 *     End of CBBCSD
 1082 *
 1083       END
 1084