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