SRC/dbbcsd.f

   1 *> \brief \b DBBCSD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DBBCSD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbbcsd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbbcsd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbbcsd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DBBCSD( 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, WORK, LWORK, INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
  28 *       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       DOUBLE PRECISION   B11D( * ), B11E( * ), B12D( * ), B12E( * ),
  32 *      $                   B21D( * ), B21E( * ), B22D( * ), B22E( * ),
  33 *      $                   PHI( * ), THETA( * ), WORK( * )
  34 *       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
  35 *      $                   V2T( LDV2T, * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> DBBCSD computes the CS decomposition of an orthogonal 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 |    ]**T
  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 DORCSD 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 orthogonal 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 orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDV1T,Q)
 180 *>          On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
 181 *>          by the 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 DOUBLE PRECISION array, dimenison (LDV2T,M-Q)
 194 *>          On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
 195 *>          premultiplied by the 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 DOUBLE PRECISION array, dimension (Q)
 209 *>          When DBBCSD converges, B11D contains the cosines of THETA(1),
 210 *>          ..., THETA(Q). If DBBCSD 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 DOUBLE PRECISION array, dimension (Q-1)
 218 *>          When DBBCSD converges, B11E contains zeros. If DBBCSD 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 DOUBLE PRECISION array, dimension (Q)
 226 *>          When DBBCSD converges, B12D contains the negative sines of
 227 *>          THETA(1), ..., THETA(Q). If DBBCSD 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 DOUBLE PRECISION array, dimension (Q-1)
 235 *>          When DBBCSD converges, B12E contains zeros. If DBBCSD 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 DOUBLE PRECISION  array, dimension (Q)
 243 *>          When DBBCSD converges, B21D contains the negative sines of
 244 *>          THETA(1), ..., THETA(Q). If DBBCSD 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 DOUBLE PRECISION  array, dimension (Q-1)
 252 *>          When DBBCSD converges, B21E contains zeros. If DBBCSD 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 DOUBLE PRECISION  array, dimension (Q)
 260 *>          When DBBCSD converges, B22D contains the negative sines of
 261 *>          THETA(1), ..., THETA(Q). If DBBCSD 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 DOUBLE PRECISION  array, dimension (Q-1)
 269 *>          When DBBCSD converges, B22E contains zeros. If DBBCSD fails
 270 *>          to converge, then B22E contains the subdiagonal of the
 271 *>          partially reduced bottom-right block.
 272 *> \endverbatim
 273 *>
 274 *> \param[out] WORK
 275 *> \verbatim
 276 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 277 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 278 *> \endverbatim
 279 *>
 280 *> \param[in] LWORK
 281 *> \verbatim
 282 *>          LWORK is INTEGER
 283 *>          The dimension of the array WORK. LWORK >= MAX(1,8*Q).
 284 *>
 285 *>          If LWORK = -1, then a workspace query is assumed; the
 286 *>          routine only calculates the optimal size of the WORK array,
 287 *>          returns this value as the first entry of the work array, and
 288 *>          no error message related to LWORK 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 DBBCSD 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  DOUBLE PRECISION, 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 *> \date June 2016
 326 *
 327 *> \ingroup doubleOTHERcomputational
 328 *
 329 *  =====================================================================
 330       SUBROUTINE DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
 331      $                   THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
 332      $                   V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
 333      $                   B22D, B22E, WORK, LWORK, INFO )
 334 *
 335 *  -- LAPACK computational routine (version 3.6.1) --
 336 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 337 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 338 *     June 2016
 339 *
 340 *     .. Scalar Arguments ..
 341       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
 342       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
 343 *     ..
 344 *     .. Array Arguments ..
 345       DOUBLE PRECISION   B11D( * ), B11E( * ), B12D( * ), B12E( * ),
 346      $                   B21D( * ), B21E( * ), B22D( * ), B22E( * ),
 347      $                   PHI( * ), THETA( * ), WORK( * )
 348       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
 349      $                   V2T( LDV2T, * )
 350 *     ..
 351 *
 352 *  ===================================================================
 353 *
 354 *     .. Parameters ..
 355       INTEGER            MAXITR
 356       PARAMETER          ( MAXITR = 6 )
 357       DOUBLE PRECISION   HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO
 358       PARAMETER          ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0,
 359      $                     ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0,
 360      $                     TEN = 10.0D0, ZERO = 0.0D0 )
 361       DOUBLE PRECISION   NEGONE
 362       PARAMETER          ( NEGONE = -1.0D0 )
 363 *     ..
 364 *     .. Local Scalars ..
 365       LOGICAL            COLMAJOR, LQUERY, RESTART11, RESTART12,
 366      $                   RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
 367      $                   WANTV2T
 368       INTEGER            I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
 369      $                   IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
 370      $                   LWORKMIN, LWORKOPT, MAXIT, MINI
 371       DOUBLE PRECISION   B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
 372      $                   EPS, MU, NU, R, SIGMA11, SIGMA21,
 373      $                   TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
 374      $                   UNFL, X1, X2, Y1, Y2
 375 *
 376 *     .. External Subroutines ..
 377       EXTERNAL           DLASR, DSCAL, DSWAP, DLARTGP, DLARTGS, DLAS2,
 378      $                   XERBLA
 379 *     ..
 380 *     .. External Functions ..
 381       DOUBLE PRECISION   DLAMCH
 382       LOGICAL            LSAME
 383       EXTERNAL           LSAME, DLAMCH
 384 *     ..
 385 *     .. Intrinsic Functions ..
 386       INTRINSIC          ABS, ATAN2, COS, MAX, MIN, SIN, SQRT
 387 *     ..
 388 *     .. Executable Statements ..
 389 *
 390 *     Test input arguments
 391 *
 392       INFO = 0
 393       LQUERY = LWORK .EQ. -1
 394       WANTU1 = LSAME( JOBU1, 'Y' )
 395       WANTU2 = LSAME( JOBU2, 'Y' )
 396       WANTV1T = LSAME( JOBV1T, 'Y' )
 397       WANTV2T = LSAME( JOBV2T, 'Y' )
 398       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
 399 *
 400       IF( M .LT. 0 ) THEN
 401          INFO = -6
 402       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
 403          INFO = -7
 404       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
 405          INFO = -8
 406       ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN
 407          INFO = -8
 408       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
 409          INFO = -12
 410       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
 411          INFO = -14
 412       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
 413          INFO = -16
 414       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
 415          INFO = -18
 416       END IF
 417 *
 418 *     Quick return if Q = 0
 419 *
 420       IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN
 421          LWORKMIN = 1
 422          WORK(1) = LWORKMIN
 423          RETURN
 424       END IF
 425 *
 426 *     Compute workspace
 427 *
 428       IF( INFO .EQ. 0 ) THEN
 429          IU1CS = 1
 430          IU1SN = IU1CS + Q
 431          IU2CS = IU1SN + Q
 432          IU2SN = IU2CS + Q
 433          IV1TCS = IU2SN + Q
 434          IV1TSN = IV1TCS + Q
 435          IV2TCS = IV1TSN + Q
 436          IV2TSN = IV2TCS + Q
 437          LWORKOPT = IV2TSN + Q - 1
 438          LWORKMIN = LWORKOPT
 439          WORK(1) = LWORKOPT
 440          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
 441             INFO = -28
 442          END IF
 443       END IF
 444 *
 445       IF( INFO .NE. 0 ) THEN
 446          CALL XERBLA( 'DBBCSD', -INFO )
 447          RETURN
 448       ELSE IF( LQUERY ) THEN
 449          RETURN
 450       END IF
 451 *
 452 *     Get machine constants
 453 *
 454       EPS = DLAMCH( 'Epsilon' )
 455       UNFL = DLAMCH( 'Safe minimum' )
 456       TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) )
 457       TOL = TOLMUL*EPS
 458       THRESH = MAX( TOL, MAXITR*Q*Q*UNFL )
 459 *
 460 *     Test for negligible sines or cosines
 461 *
 462       DO I = 1, Q
 463          IF( THETA(I) .LT. THRESH ) THEN
 464             THETA(I) = ZERO
 465          ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
 466             THETA(I) = PIOVER2
 467          END IF
 468       END DO
 469       DO I = 1, Q-1
 470          IF( PHI(I) .LT. THRESH ) THEN
 471             PHI(I) = ZERO
 472          ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
 473             PHI(I) = PIOVER2
 474          END IF
 475       END DO
 476 *
 477 *     Initial deflation
 478 *
 479       IMAX = Q
 480       DO WHILE( IMAX .GT. 1 )
 481          IF( PHI(IMAX-1) .NE. ZERO ) THEN
 482             EXIT
 483          END IF
 484          IMAX = IMAX - 1
 485       END DO
 486       IMIN = IMAX - 1
 487       IF  ( IMIN .GT. 1 ) THEN
 488          DO WHILE( PHI(IMIN-1) .NE. ZERO )
 489             IMIN = IMIN - 1
 490             IF  ( IMIN .LE. 1 ) EXIT
 491          END DO
 492       END IF
 493 *
 494 *     Initialize iteration counter
 495 *
 496       MAXIT = MAXITR*Q*Q
 497       ITER = 0
 498 *
 499 *     Begin main iteration loop
 500 *
 501       DO WHILE( IMAX .GT. 1 )
 502 *
 503 *        Compute the matrix entries
 504 *
 505          B11D(IMIN) = COS( THETA(IMIN) )
 506          B21D(IMIN) = -SIN( THETA(IMIN) )
 507          DO I = IMIN, IMAX - 1
 508             B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) )
 509             B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) )
 510             B12D(I) = SIN( THETA(I) ) * COS( PHI(I) )
 511             B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) )
 512             B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) )
 513             B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) )
 514             B22D(I) = COS( THETA(I) ) * COS( PHI(I) )
 515             B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) )
 516          END DO
 517          B12D(IMAX) = SIN( THETA(IMAX) )
 518          B22D(IMAX) = COS( THETA(IMAX) )
 519 *
 520 *        Abort if not converging; otherwise, increment ITER
 521 *
 522          IF( ITER .GT. MAXIT ) THEN
 523             INFO = 0
 524             DO I = 1, Q
 525                IF( PHI(I) .NE. ZERO )
 526      $            INFO = INFO + 1
 527             END DO
 528             RETURN
 529          END IF
 530 *
 531          ITER = ITER + IMAX - IMIN
 532 *
 533 *        Compute shifts
 534 *
 535          THETAMAX = THETA(IMIN)
 536          THETAMIN = THETA(IMIN)
 537          DO I = IMIN+1, IMAX
 538             IF( THETA(I) > THETAMAX )
 539      $         THETAMAX = THETA(I)
 540             IF( THETA(I) < THETAMIN )
 541      $         THETAMIN = THETA(I)
 542          END DO
 543 *
 544          IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN
 545 *
 546 *           Zero on diagonals of B11 and B22; induce deflation with a
 547 *           zero shift
 548 *
 549             MU = ZERO
 550             NU = ONE
 551 *
 552          ELSE IF( THETAMIN .LT. THRESH ) THEN
 553 *
 554 *           Zero on diagonals of B12 and B22; induce deflation with a
 555 *           zero shift
 556 *
 557             MU = ONE
 558             NU = ZERO
 559 *
 560          ELSE
 561 *
 562 *           Compute shifts for B11 and B21 and use the lesser
 563 *
 564             CALL DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11,
 565      $                  DUMMY )
 566             CALL DLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21,
 567      $                  DUMMY )
 568 *
 569             IF( SIGMA11 .LE. SIGMA21 ) THEN
 570                MU = SIGMA11
 571                NU = SQRT( ONE - MU**2 )
 572                IF( MU .LT. THRESH ) THEN
 573                   MU = ZERO
 574                   NU = ONE
 575                END IF
 576             ELSE
 577                NU = SIGMA21
 578                MU = SQRT( 1.0 - NU**2 )
 579                IF( NU .LT. THRESH ) THEN
 580                   MU = ONE
 581                   NU = ZERO
 582                END IF
 583             END IF
 584          END IF
 585 *
 586 *        Rotate to produce bulges in B11 and B21
 587 *
 588          IF( MU .LE. NU ) THEN
 589             CALL DLARTGS( B11D(IMIN), B11E(IMIN), MU,
 590      $                    WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1) )
 591          ELSE
 592             CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU,
 593      $                    WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1) )
 594          END IF
 595 *
 596          TEMP = WORK(IV1TCS+IMIN-1)*B11D(IMIN) +
 597      $          WORK(IV1TSN+IMIN-1)*B11E(IMIN)
 598          B11E(IMIN) = WORK(IV1TCS+IMIN-1)*B11E(IMIN) -
 599      $                WORK(IV1TSN+IMIN-1)*B11D(IMIN)
 600          B11D(IMIN) = TEMP
 601          B11BULGE = WORK(IV1TSN+IMIN-1)*B11D(IMIN+1)
 602          B11D(IMIN+1) = WORK(IV1TCS+IMIN-1)*B11D(IMIN+1)
 603          TEMP = WORK(IV1TCS+IMIN-1)*B21D(IMIN) +
 604      $          WORK(IV1TSN+IMIN-1)*B21E(IMIN)
 605          B21E(IMIN) = WORK(IV1TCS+IMIN-1)*B21E(IMIN) -
 606      $                WORK(IV1TSN+IMIN-1)*B21D(IMIN)
 607          B21D(IMIN) = TEMP
 608          B21BULGE = WORK(IV1TSN+IMIN-1)*B21D(IMIN+1)
 609          B21D(IMIN+1) = WORK(IV1TCS+IMIN-1)*B21D(IMIN+1)
 610 *
 611 *        Compute THETA(IMIN)
 612 *
 613          THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ),
 614      $                   SQRT( B11D(IMIN)**2+B11BULGE**2 ) )
 615 *
 616 *        Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
 617 *
 618          IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
 619             CALL DLARTGP( B11BULGE, B11D(IMIN), WORK(IU1SN+IMIN-1),
 620      $                    WORK(IU1CS+IMIN-1), R )
 621          ELSE IF( MU .LE. NU ) THEN
 622             CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU,
 623      $                    WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) )
 624          ELSE
 625             CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
 626      $                    WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) )
 627          END IF
 628          IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
 629             CALL DLARTGP( B21BULGE, B21D(IMIN), WORK(IU2SN+IMIN-1),
 630      $                    WORK(IU2CS+IMIN-1), R )
 631          ELSE IF( NU .LT. MU ) THEN
 632             CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU,
 633      $                    WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1) )
 634          ELSE
 635             CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU,
 636      $                    WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1) )
 637          END IF
 638          WORK(IU2CS+IMIN-1) = -WORK(IU2CS+IMIN-1)
 639          WORK(IU2SN+IMIN-1) = -WORK(IU2SN+IMIN-1)
 640 *
 641          TEMP = WORK(IU1CS+IMIN-1)*B11E(IMIN) +
 642      $          WORK(IU1SN+IMIN-1)*B11D(IMIN+1)
 643          B11D(IMIN+1) = WORK(IU1CS+IMIN-1)*B11D(IMIN+1) -
 644      $                  WORK(IU1SN+IMIN-1)*B11E(IMIN)
 645          B11E(IMIN) = TEMP
 646          IF( IMAX .GT. IMIN+1 ) THEN
 647             B11BULGE = WORK(IU1SN+IMIN-1)*B11E(IMIN+1)
 648             B11E(IMIN+1) = WORK(IU1CS+IMIN-1)*B11E(IMIN+1)
 649          END IF
 650          TEMP = WORK(IU1CS+IMIN-1)*B12D(IMIN) +
 651      $          WORK(IU1SN+IMIN-1)*B12E(IMIN)
 652          B12E(IMIN) = WORK(IU1CS+IMIN-1)*B12E(IMIN) -
 653      $                WORK(IU1SN+IMIN-1)*B12D(IMIN)
 654          B12D(IMIN) = TEMP
 655          B12BULGE = WORK(IU1SN+IMIN-1)*B12D(IMIN+1)
 656          B12D(IMIN+1) = WORK(IU1CS+IMIN-1)*B12D(IMIN+1)
 657          TEMP = WORK(IU2CS+IMIN-1)*B21E(IMIN) +
 658      $          WORK(IU2SN+IMIN-1)*B21D(IMIN+1)
 659          B21D(IMIN+1) = WORK(IU2CS+IMIN-1)*B21D(IMIN+1) -
 660      $                  WORK(IU2SN+IMIN-1)*B21E(IMIN)
 661          B21E(IMIN) = TEMP
 662          IF( IMAX .GT. IMIN+1 ) THEN
 663             B21BULGE = WORK(IU2SN+IMIN-1)*B21E(IMIN+1)
 664             B21E(IMIN+1) = WORK(IU2CS+IMIN-1)*B21E(IMIN+1)
 665          END IF
 666          TEMP = WORK(IU2CS+IMIN-1)*B22D(IMIN) +
 667      $          WORK(IU2SN+IMIN-1)*B22E(IMIN)
 668          B22E(IMIN) = WORK(IU2CS+IMIN-1)*B22E(IMIN) -
 669      $                WORK(IU2SN+IMIN-1)*B22D(IMIN)
 670          B22D(IMIN) = TEMP
 671          B22BULGE = WORK(IU2SN+IMIN-1)*B22D(IMIN+1)
 672          B22D(IMIN+1) = WORK(IU2CS+IMIN-1)*B22D(IMIN+1)
 673 *
 674 *        Inner loop: chase bulges from B11(IMIN,IMIN+2),
 675 *        B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
 676 *        bottom-right
 677 *
 678          DO I = IMIN+1, IMAX-1
 679 *
 680 *           Compute PHI(I-1)
 681 *
 682             X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1)
 683             X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE
 684             Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1)
 685             Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE
 686 *
 687             PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) )
 688 *
 689 *           Determine if there are bulges to chase or if a new direct
 690 *           summand has been reached
 691 *
 692             RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
 693             RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
 694             RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
 695             RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
 696 *
 697 *           If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
 698 *           B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
 699 *           chasing by applying the original shift again.
 700 *
 701             IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN
 702                CALL DLARTGP( X2, X1, WORK(IV1TSN+I-1), WORK(IV1TCS+I-1),
 703      $                       R )
 704             ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN
 705                CALL DLARTGP( B11BULGE, B11E(I-1), WORK(IV1TSN+I-1),
 706      $                       WORK(IV1TCS+I-1), R )
 707             ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN
 708                CALL DLARTGP( B21BULGE, B21E(I-1), WORK(IV1TSN+I-1),
 709      $                       WORK(IV1TCS+I-1), R )
 710             ELSE IF( MU .LE. NU ) THEN
 711                CALL DLARTGS( B11D(I), B11E(I), MU, WORK(IV1TCS+I-1),
 712      $                       WORK(IV1TSN+I-1) )
 713             ELSE
 714                CALL DLARTGS( B21D(I), B21E(I), NU, WORK(IV1TCS+I-1),
 715      $                       WORK(IV1TSN+I-1) )
 716             END IF
 717             WORK(IV1TCS+I-1) = -WORK(IV1TCS+I-1)
 718             WORK(IV1TSN+I-1) = -WORK(IV1TSN+I-1)
 719             IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
 720                CALL DLARTGP( Y2, Y1, WORK(IV2TSN+I-1-1),
 721      $                       WORK(IV2TCS+I-1-1), R )
 722             ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
 723                CALL DLARTGP( B12BULGE, B12D(I-1), WORK(IV2TSN+I-1-1),
 724      $                       WORK(IV2TCS+I-1-1), R )
 725             ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
 726                CALL DLARTGP( B22BULGE, B22D(I-1), WORK(IV2TSN+I-1-1),
 727      $                       WORK(IV2TCS+I-1-1), R )
 728             ELSE IF( NU .LT. MU ) THEN
 729                CALL DLARTGS( B12E(I-1), B12D(I), NU, WORK(IV2TCS+I-1-1),
 730      $                       WORK(IV2TSN+I-1-1) )
 731             ELSE
 732                CALL DLARTGS( B22E(I-1), B22D(I), MU, WORK(IV2TCS+I-1-1),
 733      $                       WORK(IV2TSN+I-1-1) )
 734             END IF
 735 *
 736             TEMP = WORK(IV1TCS+I-1)*B11D(I) + WORK(IV1TSN+I-1)*B11E(I)
 737             B11E(I) = WORK(IV1TCS+I-1)*B11E(I) -
 738      $                WORK(IV1TSN+I-1)*B11D(I)
 739             B11D(I) = TEMP
 740             B11BULGE = WORK(IV1TSN+I-1)*B11D(I+1)
 741             B11D(I+1) = WORK(IV1TCS+I-1)*B11D(I+1)
 742             TEMP = WORK(IV1TCS+I-1)*B21D(I) + WORK(IV1TSN+I-1)*B21E(I)
 743             B21E(I) = WORK(IV1TCS+I-1)*B21E(I) -
 744      $                WORK(IV1TSN+I-1)*B21D(I)
 745             B21D(I) = TEMP
 746             B21BULGE = WORK(IV1TSN+I-1)*B21D(I+1)
 747             B21D(I+1) = WORK(IV1TCS+I-1)*B21D(I+1)
 748             TEMP = WORK(IV2TCS+I-1-1)*B12E(I-1) +
 749      $             WORK(IV2TSN+I-1-1)*B12D(I)
 750             B12D(I) = WORK(IV2TCS+I-1-1)*B12D(I) -
 751      $                WORK(IV2TSN+I-1-1)*B12E(I-1)
 752             B12E(I-1) = TEMP
 753             B12BULGE = WORK(IV2TSN+I-1-1)*B12E(I)
 754             B12E(I) = WORK(IV2TCS+I-1-1)*B12E(I)
 755             TEMP = WORK(IV2TCS+I-1-1)*B22E(I-1) +
 756      $             WORK(IV2TSN+I-1-1)*B22D(I)
 757             B22D(I) = WORK(IV2TCS+I-1-1)*B22D(I) -
 758      $                WORK(IV2TSN+I-1-1)*B22E(I-1)
 759             B22E(I-1) = TEMP
 760             B22BULGE = WORK(IV2TSN+I-1-1)*B22E(I)
 761             B22E(I) = WORK(IV2TCS+I-1-1)*B22E(I)
 762 *
 763 *           Compute THETA(I)
 764 *
 765             X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1)
 766             X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE
 767             Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1)
 768             Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE
 769 *
 770             THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) )
 771 *
 772 *           Determine if there are bulges to chase or if a new direct
 773 *           summand has been reached
 774 *
 775             RESTART11 =   B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
 776             RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
 777             RESTART21 =   B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
 778             RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
 779 *
 780 *           If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
 781 *           B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
 782 *           chasing by applying the original shift again.
 783 *
 784             IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN
 785                CALL DLARTGP( X2, X1, WORK(IU1SN+I-1), WORK(IU1CS+I-1),
 786      $                       R )
 787             ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN
 788                CALL DLARTGP( B11BULGE, B11D(I), WORK(IU1SN+I-1),
 789      $                       WORK(IU1CS+I-1), R )
 790             ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN
 791                CALL DLARTGP( B12BULGE, B12E(I-1), WORK(IU1SN+I-1),
 792      $                       WORK(IU1CS+I-1), R )
 793             ELSE IF( MU .LE. NU ) THEN
 794                CALL DLARTGS( B11E(I), B11D(I+1), MU, WORK(IU1CS+I-1),
 795      $                       WORK(IU1SN+I-1) )
 796             ELSE
 797                CALL DLARTGS( B12D(I), B12E(I), NU, WORK(IU1CS+I-1),
 798      $                       WORK(IU1SN+I-1) )
 799             END IF
 800             IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN
 801                CALL DLARTGP( Y2, Y1, WORK(IU2SN+I-1), WORK(IU2CS+I-1),
 802      $                       R )
 803             ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN
 804                CALL DLARTGP( B21BULGE, B21D(I), WORK(IU2SN+I-1),
 805      $                       WORK(IU2CS+I-1), R )
 806             ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN
 807                CALL DLARTGP( B22BULGE, B22E(I-1), WORK(IU2SN+I-1),
 808      $                       WORK(IU2CS+I-1), R )
 809             ELSE IF( NU .LT. MU ) THEN
 810                CALL DLARTGS( B21E(I), B21E(I+1), NU, WORK(IU2CS+I-1),
 811      $                       WORK(IU2SN+I-1) )
 812             ELSE
 813                CALL DLARTGS( B22D(I), B22E(I), MU, WORK(IU2CS+I-1),
 814      $                       WORK(IU2SN+I-1) )
 815             END IF
 816             WORK(IU2CS+I-1) = -WORK(IU2CS+I-1)
 817             WORK(IU2SN+I-1) = -WORK(IU2SN+I-1)
 818 *
 819             TEMP = WORK(IU1CS+I-1)*B11E(I) + WORK(IU1SN+I-1)*B11D(I+1)
 820             B11D(I+1) = WORK(IU1CS+I-1)*B11D(I+1) -
 821      $                  WORK(IU1SN+I-1)*B11E(I)
 822             B11E(I) = TEMP
 823             IF( I .LT. IMAX - 1 ) THEN
 824                B11BULGE = WORK(IU1SN+I-1)*B11E(I+1)
 825                B11E(I+1) = WORK(IU1CS+I-1)*B11E(I+1)
 826             END IF
 827             TEMP = WORK(IU2CS+I-1)*B21E(I) + WORK(IU2SN+I-1)*B21D(I+1)
 828             B21D(I+1) = WORK(IU2CS+I-1)*B21D(I+1) -
 829      $                  WORK(IU2SN+I-1)*B21E(I)
 830             B21E(I) = TEMP
 831             IF( I .LT. IMAX - 1 ) THEN
 832                B21BULGE = WORK(IU2SN+I-1)*B21E(I+1)
 833                B21E(I+1) = WORK(IU2CS+I-1)*B21E(I+1)
 834             END IF
 835             TEMP = WORK(IU1CS+I-1)*B12D(I) + WORK(IU1SN+I-1)*B12E(I)
 836             B12E(I) = WORK(IU1CS+I-1)*B12E(I) - WORK(IU1SN+I-1)*B12D(I)
 837             B12D(I) = TEMP
 838             B12BULGE = WORK(IU1SN+I-1)*B12D(I+1)
 839             B12D(I+1) = WORK(IU1CS+I-1)*B12D(I+1)
 840             TEMP = WORK(IU2CS+I-1)*B22D(I) + WORK(IU2SN+I-1)*B22E(I)
 841             B22E(I) = WORK(IU2CS+I-1)*B22E(I) - WORK(IU2SN+I-1)*B22D(I)
 842             B22D(I) = TEMP
 843             B22BULGE = WORK(IU2SN+I-1)*B22D(I+1)
 844             B22D(I+1) = WORK(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 DLARTGP( Y2, Y1, WORK(IV2TSN+IMAX-1-1),
 865      $                    WORK(IV2TCS+IMAX-1-1), R )
 866          ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
 867             CALL DLARTGP( B12BULGE, B12D(IMAX-1), WORK(IV2TSN+IMAX-1-1),
 868      $                    WORK(IV2TCS+IMAX-1-1), R )
 869          ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
 870             CALL DLARTGP( B22BULGE, B22D(IMAX-1), WORK(IV2TSN+IMAX-1-1),
 871      $                    WORK(IV2TCS+IMAX-1-1), R )
 872          ELSE IF( NU .LT. MU ) THEN
 873             CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU,
 874      $                    WORK(IV2TCS+IMAX-1-1), WORK(IV2TSN+IMAX-1-1) )
 875          ELSE
 876             CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU,
 877      $                    WORK(IV2TCS+IMAX-1-1), WORK(IV2TSN+IMAX-1-1) )
 878          END IF
 879 *
 880          TEMP = WORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) +
 881      $          WORK(IV2TSN+IMAX-1-1)*B12D(IMAX)
 882          B12D(IMAX) = WORK(IV2TCS+IMAX-1-1)*B12D(IMAX) -
 883      $                WORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1)
 884          B12E(IMAX-1) = TEMP
 885          TEMP = WORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) +
 886      $          WORK(IV2TSN+IMAX-1-1)*B22D(IMAX)
 887          B22D(IMAX) = WORK(IV2TCS+IMAX-1-1)*B22D(IMAX) -
 888      $                WORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1)
 889          B22E(IMAX-1) = TEMP
 890 *
 891 *        Update singular vectors
 892 *
 893          IF( WANTU1 ) THEN
 894             IF( COLMAJOR ) THEN
 895                CALL DLASR( 'R', 'V', 'F', P, IMAX-IMIN+1,
 896      $                     WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1),
 897      $                     U1(1,IMIN), LDU1 )
 898             ELSE
 899                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, P,
 900      $                     WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1),
 901      $                     U1(IMIN,1), LDU1 )
 902             END IF
 903          END IF
 904          IF( WANTU2 ) THEN
 905             IF( COLMAJOR ) THEN
 906                CALL DLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1,
 907      $                     WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1),
 908      $                     U2(1,IMIN), LDU2 )
 909             ELSE
 910                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P,
 911      $                     WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1),
 912      $                     U2(IMIN,1), LDU2 )
 913             END IF
 914          END IF
 915          IF( WANTV1T ) THEN
 916             IF( COLMAJOR ) THEN
 917                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q,
 918      $                     WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1),
 919      $                     V1T(IMIN,1), LDV1T )
 920             ELSE
 921                CALL DLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1,
 922      $                     WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1),
 923      $                     V1T(1,IMIN), LDV1T )
 924             END IF
 925          END IF
 926          IF( WANTV2T ) THEN
 927             IF( COLMAJOR ) THEN
 928                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q,
 929      $                     WORK(IV2TCS+IMIN-1), WORK(IV2TSN+IMIN-1),
 930      $                     V2T(IMIN,1), LDV2T )
 931             ELSE
 932                CALL DLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1,
 933      $                     WORK(IV2TCS+IMIN-1), WORK(IV2TSN+IMIN-1),
 934      $                     V2T(1,IMIN), LDV2T )
 935             END IF
 936          END IF
 937 *
 938 *        Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
 939 *
 940          IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN
 941             B11D(IMAX) = -B11D(IMAX)
 942             B21D(IMAX) = -B21D(IMAX)
 943             IF( WANTV1T ) THEN
 944                IF( COLMAJOR ) THEN
 945                   CALL DSCAL( Q, NEGONE, V1T(IMAX,1), LDV1T )
 946                ELSE
 947                   CALL DSCAL( Q, NEGONE, V1T(1,IMAX), 1 )
 948                END IF
 949             END IF
 950          END IF
 951 *
 952 *        Compute THETA(IMAX)
 953 *
 954          X1 = COS(PHI(IMAX-1))*B11D(IMAX) +
 955      $        SIN(PHI(IMAX-1))*B12E(IMAX-1)
 956          Y1 = COS(PHI(IMAX-1))*B21D(IMAX) +
 957      $        SIN(PHI(IMAX-1))*B22E(IMAX-1)
 958 *
 959          THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) )
 960 *
 961 *        Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
 962 *        and B22(IMAX,IMAX-1)
 963 *
 964          IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN
 965             B12D(IMAX) = -B12D(IMAX)
 966             IF( WANTU1 ) THEN
 967                IF( COLMAJOR ) THEN
 968                   CALL DSCAL( P, NEGONE, U1(1,IMAX), 1 )
 969                ELSE
 970                   CALL DSCAL( P, NEGONE, U1(IMAX,1), LDU1 )
 971                END IF
 972             END IF
 973          END IF
 974          IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN
 975             B22D(IMAX) = -B22D(IMAX)
 976             IF( WANTU2 ) THEN
 977                IF( COLMAJOR ) THEN
 978                   CALL DSCAL( M-P, NEGONE, U2(1,IMAX), 1 )
 979                ELSE
 980                   CALL DSCAL( M-P, NEGONE, U2(IMAX,1), LDU2 )
 981                END IF
 982             END IF
 983          END IF
 984 *
 985 *        Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
 986 *
 987          IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
 988             IF( WANTV2T ) THEN
 989                IF( COLMAJOR ) THEN
 990                   CALL DSCAL( M-Q, NEGONE, V2T(IMAX,1), LDV2T )
 991                ELSE
 992                   CALL DSCAL( M-Q, NEGONE, V2T(1,IMAX), 1 )
 993                END IF
 994             END IF
 995          END IF
 996 *
 997 *        Test for negligible sines or cosines
 998 *
 999          DO I = IMIN, IMAX
1000             IF( THETA(I) .LT. THRESH ) THEN
1001                THETA(I) = ZERO
1002             ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
1003                THETA(I) = PIOVER2
1004             END IF
1005          END DO
1006          DO I = IMIN, IMAX-1
1007             IF( PHI(I) .LT. THRESH ) THEN
1008                PHI(I) = ZERO
1009             ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
1010                PHI(I) = PIOVER2
1011             END IF
1012          END DO
1013 *
1014 *        Deflate
1015 *
1016          IF (IMAX .GT. 1) THEN
1017             DO WHILE( PHI(IMAX-1) .EQ. ZERO )
1018                IMAX = IMAX - 1
1019                IF (IMAX .LE. 1) EXIT
1020             END DO
1021          END IF
1022          IF( IMIN .GT. IMAX - 1 )
1023      $      IMIN = IMAX - 1
1024          IF (IMIN .GT. 1) THEN
1025             DO WHILE (PHI(IMIN-1) .NE. ZERO)
1026                 IMIN = IMIN - 1
1027                 IF (IMIN .LE. 1) EXIT
1028             END DO
1029          END IF
1030 *
1031 *        Repeat main iteration loop
1032 *
1033       END DO
1034 *
1035 *     Postprocessing: order THETA from least to greatest
1036 *
1037       DO I = 1, Q
1038 *
1039          MINI = I
1040          THETAMIN = THETA(I)
1041          DO J = I+1, Q
1042             IF( THETA(J) .LT. THETAMIN ) THEN
1043                MINI = J
1044                THETAMIN = THETA(J)
1045             END IF
1046          END DO
1047 *
1048          IF( MINI .NE. I ) THEN
1049             THETA(MINI) = THETA(I)
1050             THETA(I) = THETAMIN
1051             IF( COLMAJOR ) THEN
1052                IF( WANTU1 )
1053      $            CALL DSWAP( P, U1(1,I), 1, U1(1,MINI), 1 )
1054                IF( WANTU2 )
1055      $            CALL DSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 )
1056                IF( WANTV1T )
1057      $            CALL DSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T )
1058                IF( WANTV2T )
1059      $            CALL DSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1),
1060      $               LDV2T )
1061             ELSE
1062                IF( WANTU1 )
1063      $            CALL DSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 )
1064                IF( WANTU2 )
1065      $            CALL DSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 )
1066                IF( WANTV1T )
1067      $            CALL DSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 )
1068                IF( WANTV2T )
1069      $            CALL DSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 )
1070             END IF
1071          END IF
1072 *
1073       END DO
1074 *
1075       RETURN
1076 *
1077 *     End of DBBCSD
1078 *
1079       END
1080