LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dbbcsd.f
Go to the documentation of this file.
00001 *> \brief \b DBBCSD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DBBCSD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbbcsd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbbcsd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbbcsd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DBBCSD( 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, WORK, LWORK, INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
00028 *       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       DOUBLE PRECISION   B11D( * ), B11E( * ), B12D( * ), B12E( * ),
00032 *      $                   B21D( * ), B21E( * ), B22D( * ), B22E( * ),
00033 *      $                   PHI( * ), THETA( * ), WORK( * )
00034 *       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00035 *      $                   V2T( LDV2T, * )
00036 *       ..
00037 *  
00038 *
00039 *> \par Purpose:
00040 *  =============
00041 *>
00042 *> \verbatim
00043 *>
00044 *> DBBCSD computes the CS decomposition of an orthogonal 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 |    ]**T
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 DORCSD 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 orthogonal 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 orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDV1T,Q)
00180 *>          On entry, a LDV1T-by-Q matrix. On exit, V1T is premultiplied
00181 *>          by the 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 DOUBLE PRECISION array, dimenison (LDV2T,M-Q)
00194 *>          On entry, a LDV2T-by-(M-Q) matrix. On exit, V2T is
00195 *>          premultiplied by the 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 DBBCSD converges, B11D contains the cosines of THETA(1),
00210 *>          ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B11E contains zeros. If DBBCSD 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 DBBCSD converges, B12D contains the negative sines of
00227 *>          THETA(1), ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B12E contains zeros. If DBBCSD 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] WORK
00275 *> \verbatim
00276 *>          WORK 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] LWORK
00281 *> \verbatim
00282 *>          LWORK is INTEGER
00283 *>          The dimension of the array WORK. LWORK >= MAX(1,8*Q).
00284 *>
00285 *>          If LWORK = -1, then a workspace query is assumed; the
00286 *>          routine only calculates the optimal size of the WORK array,
00287 *>          returns this value as the first entry of the work array, and
00288 *>          no error message related to LWORK 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 DBBCSD 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 doubleOTHERcomputational
00328 *
00329 *  =====================================================================
00330       SUBROUTINE DBBCSD( 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, WORK, LWORK, 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, LWORK, M, P, Q
00343 *     ..
00344 *     .. Array Arguments ..
00345       DOUBLE PRECISION   B11D( * ), B11E( * ), B12D( * ), B12E( * ),
00346      $                   B21D( * ), B21E( * ), B22D( * ), B22E( * ),
00347      $                   PHI( * ), THETA( * ), WORK( * )
00348       DOUBLE PRECISION   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       DOUBLE PRECISION   NEGONECOMPLEX
00362       PARAMETER          ( NEGONECOMPLEX = -1.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      $                   LWORKMIN, LWORKOPT, 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 Subroutines ..
00377       EXTERNAL           DLASR, DSCAL, DSWAP, DLARTGP, DLARTGS, DLAS2,
00378      $                   XERBLA
00379 *     ..
00380 *     .. External Functions ..
00381       DOUBLE PRECISION   DLAMCH
00382       LOGICAL            LSAME
00383       EXTERNAL           LSAME, DLAMCH
00384 *     ..
00385 *     .. Intrinsic Functions ..
00386       INTRINSIC          ABS, ATAN2, COS, MAX, MIN, SIN, SQRT
00387 *     ..
00388 *     .. Executable Statements ..
00389 *
00390 *     Test input arguments
00391 *
00392       INFO = 0
00393       LQUERY = LWORK .EQ. -1
00394       WANTU1 = LSAME( JOBU1, 'Y' )
00395       WANTU2 = LSAME( JOBU2, 'Y' )
00396       WANTV1T = LSAME( JOBV1T, 'Y' )
00397       WANTV2T = LSAME( JOBV2T, 'Y' )
00398       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00399 *
00400       IF( M .LT. 0 ) THEN
00401          INFO = -6
00402       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00403          INFO = -7
00404       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
00405          INFO = -8
00406       ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN
00407          INFO = -8
00408       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
00409          INFO = -12
00410       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
00411          INFO = -14
00412       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
00413          INFO = -16
00414       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
00415          INFO = -18
00416       END IF
00417 *
00418 *     Quick return if Q = 0
00419 *
00420       IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN
00421          LWORKMIN = 1
00422          WORK(1) = LWORKMIN
00423          RETURN
00424       END IF
00425 *
00426 *     Compute workspace
00427 *
00428       IF( INFO .EQ. 0 ) THEN
00429          IU1CS = 1
00430          IU1SN = IU1CS + Q
00431          IU2CS = IU1SN + Q
00432          IU2SN = IU2CS + Q
00433          IV1TCS = IU2SN + Q
00434          IV1TSN = IV1TCS + Q
00435          IV2TCS = IV1TSN + Q
00436          IV2TSN = IV2TCS + Q
00437          LWORKOPT = IV2TSN + Q - 1
00438          LWORKMIN = LWORKOPT
00439          WORK(1) = LWORKOPT
00440          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
00441             INFO = -28
00442          END IF
00443       END IF
00444 *
00445       IF( INFO .NE. 0 ) THEN
00446          CALL XERBLA( 'DBBCSD', -INFO )
00447          RETURN
00448       ELSE IF( LQUERY ) THEN
00449          RETURN
00450       END IF
00451 *
00452 *     Get machine constants
00453 *
00454       EPS = DLAMCH( 'Epsilon' )
00455       UNFL = DLAMCH( 'Safe minimum' )
00456       TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) )
00457       TOL = TOLMUL*EPS
00458       THRESH = MAX( TOL, MAXITR*Q*Q*UNFL )
00459 *
00460 *     Test for negligible sines or cosines
00461 *
00462       DO I = 1, Q
00463          IF( THETA(I) .LT. THRESH ) THEN
00464             THETA(I) = ZERO
00465          ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
00466             THETA(I) = PIOVER2
00467          END IF
00468       END DO
00469       DO I = 1, Q-1
00470          IF( PHI(I) .LT. THRESH ) THEN
00471             PHI(I) = ZERO
00472          ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
00473             PHI(I) = PIOVER2
00474          END IF
00475       END DO
00476 *
00477 *     Initial deflation
00478 *
00479       IMAX = Q
00480       DO WHILE( ( IMAX .GT. 1 ) .AND. ( PHI(IMAX-1) .EQ. ZERO ) )
00481          IMAX = IMAX - 1
00482       END DO
00483       IMIN = IMAX - 1
00484       IF  ( IMIN .GT. 1 ) THEN
00485          DO WHILE( PHI(IMIN-1) .NE. ZERO )
00486             IMIN = IMIN - 1
00487             IF  ( IMIN .LE. 1 ) EXIT
00488          END DO
00489       END IF
00490 *
00491 *     Initialize iteration counter
00492 *
00493       MAXIT = MAXITR*Q*Q
00494       ITER = 0
00495 *
00496 *     Begin main iteration loop
00497 *
00498       DO WHILE( IMAX .GT. 1 )
00499 *
00500 *        Compute the matrix entries
00501 *
00502          B11D(IMIN) = COS( THETA(IMIN) )
00503          B21D(IMIN) = -SIN( THETA(IMIN) )
00504          DO I = IMIN, IMAX - 1
00505             B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) )
00506             B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) )
00507             B12D(I) = SIN( THETA(I) ) * COS( PHI(I) )
00508             B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) )
00509             B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) )
00510             B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) )
00511             B22D(I) = COS( THETA(I) ) * COS( PHI(I) )
00512             B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) )
00513          END DO
00514          B12D(IMAX) = SIN( THETA(IMAX) )
00515          B22D(IMAX) = COS( THETA(IMAX) )
00516 *
00517 *        Abort if not converging; otherwise, increment ITER
00518 *
00519          IF( ITER .GT. MAXIT ) THEN
00520             INFO = 0
00521             DO I = 1, Q
00522                IF( PHI(I) .NE. ZERO )
00523      $            INFO = INFO + 1
00524             END DO
00525             RETURN
00526          END IF
00527 *
00528          ITER = ITER + IMAX - IMIN
00529 *
00530 *        Compute shifts
00531 *
00532          THETAMAX = THETA(IMIN)
00533          THETAMIN = THETA(IMIN)
00534          DO I = IMIN+1, IMAX
00535             IF( THETA(I) > THETAMAX )
00536      $         THETAMAX = THETA(I)
00537             IF( THETA(I) < THETAMIN )
00538      $         THETAMIN = THETA(I)
00539          END DO
00540 *
00541          IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN
00542 *
00543 *           Zero on diagonals of B11 and B22; induce deflation with a
00544 *           zero shift
00545 *
00546             MU = ZERO
00547             NU = ONE
00548 *
00549          ELSE IF( THETAMIN .LT. THRESH ) THEN
00550 *
00551 *           Zero on diagonals of B12 and B22; induce deflation with a
00552 *           zero shift
00553 *
00554             MU = ONE
00555             NU = ZERO
00556 *
00557          ELSE
00558 *
00559 *           Compute shifts for B11 and B21 and use the lesser
00560 *
00561             CALL DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11,
00562      $                  DUMMY )
00563             CALL DLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21,
00564      $                  DUMMY )
00565 *
00566             IF( SIGMA11 .LE. SIGMA21 ) THEN
00567                MU = SIGMA11
00568                NU = SQRT( ONE - MU**2 )
00569                IF( MU .LT. THRESH ) THEN
00570                   MU = ZERO
00571                   NU = ONE
00572                END IF
00573             ELSE
00574                NU = SIGMA21
00575                MU = SQRT( 1.0 - NU**2 )
00576                IF( NU .LT. THRESH ) THEN
00577                   MU = ONE
00578                   NU = ZERO
00579                END IF
00580             END IF
00581          END IF
00582 *
00583 *        Rotate to produce bulges in B11 and B21
00584 *
00585          IF( MU .LE. NU ) THEN
00586             CALL DLARTGS( B11D(IMIN), B11E(IMIN), MU,
00587      $                    WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1) )
00588          ELSE
00589             CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU,
00590      $                    WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1) )
00591          END IF
00592 *
00593          TEMP = WORK(IV1TCS+IMIN-1)*B11D(IMIN) +
00594      $          WORK(IV1TSN+IMIN-1)*B11E(IMIN)
00595          B11E(IMIN) = WORK(IV1TCS+IMIN-1)*B11E(IMIN) -
00596      $                WORK(IV1TSN+IMIN-1)*B11D(IMIN)
00597          B11D(IMIN) = TEMP
00598          B11BULGE = WORK(IV1TSN+IMIN-1)*B11D(IMIN+1)
00599          B11D(IMIN+1) = WORK(IV1TCS+IMIN-1)*B11D(IMIN+1)
00600          TEMP = WORK(IV1TCS+IMIN-1)*B21D(IMIN) +
00601      $          WORK(IV1TSN+IMIN-1)*B21E(IMIN)
00602          B21E(IMIN) = WORK(IV1TCS+IMIN-1)*B21E(IMIN) -
00603      $                WORK(IV1TSN+IMIN-1)*B21D(IMIN)
00604          B21D(IMIN) = TEMP
00605          B21BULGE = WORK(IV1TSN+IMIN-1)*B21D(IMIN+1)
00606          B21D(IMIN+1) = WORK(IV1TCS+IMIN-1)*B21D(IMIN+1)
00607 *
00608 *        Compute THETA(IMIN)
00609 *
00610          THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ),
00611      $                   SQRT( B11D(IMIN)**2+B11BULGE**2 ) )
00612 *
00613 *        Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
00614 *
00615          IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
00616             CALL DLARTGP( B11BULGE, B11D(IMIN), WORK(IU1SN+IMIN-1),
00617      $                    WORK(IU1CS+IMIN-1), R )
00618          ELSE IF( MU .LE. NU ) THEN
00619             CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU,
00620      $                    WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) )
00621          ELSE
00622             CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
00623      $                    WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) )
00624          END IF
00625          IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
00626             CALL DLARTGP( B21BULGE, B21D(IMIN), WORK(IU2SN+IMIN-1),
00627      $                    WORK(IU2CS+IMIN-1), R )
00628          ELSE IF( NU .LT. MU ) THEN
00629             CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU,
00630      $                    WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1) )
00631          ELSE
00632             CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU,
00633      $                    WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1) )
00634          END IF
00635          WORK(IU2CS+IMIN-1) = -WORK(IU2CS+IMIN-1)
00636          WORK(IU2SN+IMIN-1) = -WORK(IU2SN+IMIN-1)
00637 *
00638          TEMP = WORK(IU1CS+IMIN-1)*B11E(IMIN) +
00639      $          WORK(IU1SN+IMIN-1)*B11D(IMIN+1)
00640          B11D(IMIN+1) = WORK(IU1CS+IMIN-1)*B11D(IMIN+1) -
00641      $                  WORK(IU1SN+IMIN-1)*B11E(IMIN)
00642          B11E(IMIN) = TEMP
00643          IF( IMAX .GT. IMIN+1 ) THEN
00644             B11BULGE = WORK(IU1SN+IMIN-1)*B11E(IMIN+1)
00645             B11E(IMIN+1) = WORK(IU1CS+IMIN-1)*B11E(IMIN+1)
00646          END IF
00647          TEMP = WORK(IU1CS+IMIN-1)*B12D(IMIN) +
00648      $          WORK(IU1SN+IMIN-1)*B12E(IMIN)
00649          B12E(IMIN) = WORK(IU1CS+IMIN-1)*B12E(IMIN) -
00650      $                WORK(IU1SN+IMIN-1)*B12D(IMIN)
00651          B12D(IMIN) = TEMP
00652          B12BULGE = WORK(IU1SN+IMIN-1)*B12D(IMIN+1)
00653          B12D(IMIN+1) = WORK(IU1CS+IMIN-1)*B12D(IMIN+1)
00654          TEMP = WORK(IU2CS+IMIN-1)*B21E(IMIN) +
00655      $          WORK(IU2SN+IMIN-1)*B21D(IMIN+1)
00656          B21D(IMIN+1) = WORK(IU2CS+IMIN-1)*B21D(IMIN+1) -
00657      $                  WORK(IU2SN+IMIN-1)*B21E(IMIN)
00658          B21E(IMIN) = TEMP
00659          IF( IMAX .GT. IMIN+1 ) THEN
00660             B21BULGE = WORK(IU2SN+IMIN-1)*B21E(IMIN+1)
00661             B21E(IMIN+1) = WORK(IU2CS+IMIN-1)*B21E(IMIN+1)
00662          END IF
00663          TEMP = WORK(IU2CS+IMIN-1)*B22D(IMIN) +
00664      $          WORK(IU2SN+IMIN-1)*B22E(IMIN)
00665          B22E(IMIN) = WORK(IU2CS+IMIN-1)*B22E(IMIN) -
00666      $                WORK(IU2SN+IMIN-1)*B22D(IMIN)
00667          B22D(IMIN) = TEMP
00668          B22BULGE = WORK(IU2SN+IMIN-1)*B22D(IMIN+1)
00669          B22D(IMIN+1) = WORK(IU2CS+IMIN-1)*B22D(IMIN+1)
00670 *
00671 *        Inner loop: chase bulges from B11(IMIN,IMIN+2),
00672 *        B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
00673 *        bottom-right
00674 *
00675          DO I = IMIN+1, IMAX-1
00676 *
00677 *           Compute PHI(I-1)
00678 *
00679             X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1)
00680             X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE
00681             Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1)
00682             Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE
00683 *
00684             PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) )
00685 *
00686 *           Determine if there are bulges to chase or if a new direct
00687 *           summand has been reached
00688 *
00689             RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
00690             RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
00691             RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
00692             RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
00693 *
00694 *           If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
00695 *           B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
00696 *           chasing by applying the original shift again.
00697 *
00698             IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN
00699                CALL DLARTGP( X2, X1, WORK(IV1TSN+I-1), WORK(IV1TCS+I-1),
00700      $                       R )
00701             ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN
00702                CALL DLARTGP( B11BULGE, B11E(I-1), WORK(IV1TSN+I-1),
00703      $                       WORK(IV1TCS+I-1), R )
00704             ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN
00705                CALL DLARTGP( B21BULGE, B21E(I-1), WORK(IV1TSN+I-1),
00706      $                       WORK(IV1TCS+I-1), R )
00707             ELSE IF( MU .LE. NU ) THEN
00708                CALL DLARTGS( B11D(I), B11E(I), MU, WORK(IV1TCS+I-1),
00709      $                       WORK(IV1TSN+I-1) )
00710             ELSE
00711                CALL DLARTGS( B21D(I), B21E(I), NU, WORK(IV1TCS+I-1),
00712      $                       WORK(IV1TSN+I-1) )
00713             END IF
00714             WORK(IV1TCS+I-1) = -WORK(IV1TCS+I-1)
00715             WORK(IV1TSN+I-1) = -WORK(IV1TSN+I-1)
00716             IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
00717                CALL DLARTGP( Y2, Y1, WORK(IV2TSN+I-1-1),
00718      $                       WORK(IV2TCS+I-1-1), R )
00719             ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
00720                CALL DLARTGP( B12BULGE, B12D(I-1), WORK(IV2TSN+I-1-1),
00721      $                       WORK(IV2TCS+I-1-1), R )
00722             ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
00723                CALL DLARTGP( B22BULGE, B22D(I-1), WORK(IV2TSN+I-1-1),
00724      $                       WORK(IV2TCS+I-1-1), R )
00725             ELSE IF( NU .LT. MU ) THEN
00726                CALL DLARTGS( B12E(I-1), B12D(I), NU, WORK(IV2TCS+I-1-1),
00727      $                       WORK(IV2TSN+I-1-1) )
00728             ELSE
00729                CALL DLARTGS( B22E(I-1), B22D(I), MU, WORK(IV2TCS+I-1-1),
00730      $                       WORK(IV2TSN+I-1-1) )
00731             END IF
00732 *
00733             TEMP = WORK(IV1TCS+I-1)*B11D(I) + WORK(IV1TSN+I-1)*B11E(I)
00734             B11E(I) = WORK(IV1TCS+I-1)*B11E(I) -
00735      $                WORK(IV1TSN+I-1)*B11D(I)
00736             B11D(I) = TEMP
00737             B11BULGE = WORK(IV1TSN+I-1)*B11D(I+1)
00738             B11D(I+1) = WORK(IV1TCS+I-1)*B11D(I+1)
00739             TEMP = WORK(IV1TCS+I-1)*B21D(I) + WORK(IV1TSN+I-1)*B21E(I)
00740             B21E(I) = WORK(IV1TCS+I-1)*B21E(I) -
00741      $                WORK(IV1TSN+I-1)*B21D(I)
00742             B21D(I) = TEMP
00743             B21BULGE = WORK(IV1TSN+I-1)*B21D(I+1)
00744             B21D(I+1) = WORK(IV1TCS+I-1)*B21D(I+1)
00745             TEMP = WORK(IV2TCS+I-1-1)*B12E(I-1) +
00746      $             WORK(IV2TSN+I-1-1)*B12D(I)
00747             B12D(I) = WORK(IV2TCS+I-1-1)*B12D(I) -
00748      $                WORK(IV2TSN+I-1-1)*B12E(I-1)
00749             B12E(I-1) = TEMP
00750             B12BULGE = WORK(IV2TSN+I-1-1)*B12E(I)
00751             B12E(I) = WORK(IV2TCS+I-1-1)*B12E(I)
00752             TEMP = WORK(IV2TCS+I-1-1)*B22E(I-1) +
00753      $             WORK(IV2TSN+I-1-1)*B22D(I)
00754             B22D(I) = WORK(IV2TCS+I-1-1)*B22D(I) -
00755      $                WORK(IV2TSN+I-1-1)*B22E(I-1)
00756             B22E(I-1) = TEMP
00757             B22BULGE = WORK(IV2TSN+I-1-1)*B22E(I)
00758             B22E(I) = WORK(IV2TCS+I-1-1)*B22E(I)
00759 *
00760 *           Compute THETA(I)
00761 *
00762             X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1)
00763             X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE
00764             Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1)
00765             Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE
00766 *
00767             THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) )
00768 *
00769 *           Determine if there are bulges to chase or if a new direct
00770 *           summand has been reached
00771 *
00772             RESTART11 =   B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
00773             RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
00774             RESTART21 =   B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
00775             RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
00776 *
00777 *           If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
00778 *           B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
00779 *           chasing by applying the original shift again.
00780 *
00781             IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN
00782                CALL DLARTGP( X2, X1, WORK(IU1SN+I-1), WORK(IU1CS+I-1),
00783      $                       R )
00784             ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN
00785                CALL DLARTGP( B11BULGE, B11D(I), WORK(IU1SN+I-1),
00786      $                       WORK(IU1CS+I-1), R )
00787             ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN
00788                CALL DLARTGP( B12BULGE, B12E(I-1), WORK(IU1SN+I-1),
00789      $                       WORK(IU1CS+I-1), R )
00790             ELSE IF( MU .LE. NU ) THEN
00791                CALL DLARTGS( B11E(I), B11D(I+1), MU, WORK(IU1CS+I-1),
00792      $                       WORK(IU1SN+I-1) )
00793             ELSE
00794                CALL DLARTGS( B12D(I), B12E(I), NU, WORK(IU1CS+I-1),
00795      $                       WORK(IU1SN+I-1) )
00796             END IF
00797             IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN
00798                CALL DLARTGP( Y2, Y1, WORK(IU2SN+I-1), WORK(IU2CS+I-1),
00799      $                       R )
00800             ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN
00801                CALL DLARTGP( B21BULGE, B21D(I), WORK(IU2SN+I-1),
00802      $                       WORK(IU2CS+I-1), R )
00803             ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN
00804                CALL DLARTGP( B22BULGE, B22E(I-1), WORK(IU2SN+I-1),
00805      $                       WORK(IU2CS+I-1), R )
00806             ELSE IF( NU .LT. MU ) THEN
00807                CALL DLARTGS( B21E(I), B21E(I+1), NU, WORK(IU2CS+I-1),
00808      $                       WORK(IU2SN+I-1) )
00809             ELSE
00810                CALL DLARTGS( B22D(I), B22E(I), MU, WORK(IU2CS+I-1),
00811      $                       WORK(IU2SN+I-1) )
00812             END IF
00813             WORK(IU2CS+I-1) = -WORK(IU2CS+I-1)
00814             WORK(IU2SN+I-1) = -WORK(IU2SN+I-1)
00815 *
00816             TEMP = WORK(IU1CS+I-1)*B11E(I) + WORK(IU1SN+I-1)*B11D(I+1)
00817             B11D(I+1) = WORK(IU1CS+I-1)*B11D(I+1) -
00818      $                  WORK(IU1SN+I-1)*B11E(I)
00819             B11E(I) = TEMP
00820             IF( I .LT. IMAX - 1 ) THEN
00821                B11BULGE = WORK(IU1SN+I-1)*B11E(I+1)
00822                B11E(I+1) = WORK(IU1CS+I-1)*B11E(I+1)
00823             END IF
00824             TEMP = WORK(IU2CS+I-1)*B21E(I) + WORK(IU2SN+I-1)*B21D(I+1)
00825             B21D(I+1) = WORK(IU2CS+I-1)*B21D(I+1) -
00826      $                  WORK(IU2SN+I-1)*B21E(I)
00827             B21E(I) = TEMP
00828             IF( I .LT. IMAX - 1 ) THEN
00829                B21BULGE = WORK(IU2SN+I-1)*B21E(I+1)
00830                B21E(I+1) = WORK(IU2CS+I-1)*B21E(I+1)
00831             END IF
00832             TEMP = WORK(IU1CS+I-1)*B12D(I) + WORK(IU1SN+I-1)*B12E(I)
00833             B12E(I) = WORK(IU1CS+I-1)*B12E(I) - WORK(IU1SN+I-1)*B12D(I)
00834             B12D(I) = TEMP
00835             B12BULGE = WORK(IU1SN+I-1)*B12D(I+1)
00836             B12D(I+1) = WORK(IU1CS+I-1)*B12D(I+1)
00837             TEMP = WORK(IU2CS+I-1)*B22D(I) + WORK(IU2SN+I-1)*B22E(I)
00838             B22E(I) = WORK(IU2CS+I-1)*B22E(I) - WORK(IU2SN+I-1)*B22D(I)
00839             B22D(I) = TEMP
00840             B22BULGE = WORK(IU2SN+I-1)*B22D(I+1)
00841             B22D(I+1) = WORK(IU2CS+I-1)*B22D(I+1)
00842 *
00843          END DO
00844 *
00845 *        Compute PHI(IMAX-1)
00846 *
00847          X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) +
00848      $        COS(THETA(IMAX-1))*B21E(IMAX-1)
00849          Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) +
00850      $        COS(THETA(IMAX-1))*B22D(IMAX-1)
00851          Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE
00852 *
00853          PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) )
00854 *
00855 *        Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
00856 *
00857          RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
00858          RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
00859 *
00860          IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
00861             CALL DLARTGP( Y2, Y1, WORK(IV2TSN+IMAX-1-1),
00862      $                    WORK(IV2TCS+IMAX-1-1), R )
00863          ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
00864             CALL DLARTGP( B12BULGE, B12D(IMAX-1), WORK(IV2TSN+IMAX-1-1),
00865      $                    WORK(IV2TCS+IMAX-1-1), R )
00866          ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
00867             CALL DLARTGP( B22BULGE, B22D(IMAX-1), WORK(IV2TSN+IMAX-1-1),
00868      $                    WORK(IV2TCS+IMAX-1-1), R )
00869          ELSE IF( NU .LT. MU ) THEN
00870             CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU,
00871      $                    WORK(IV2TCS+IMAX-1-1), WORK(IV2TSN+IMAX-1-1) )
00872          ELSE
00873             CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU,
00874      $                    WORK(IV2TCS+IMAX-1-1), WORK(IV2TSN+IMAX-1-1) )
00875          END IF
00876 *
00877          TEMP = WORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) +
00878      $          WORK(IV2TSN+IMAX-1-1)*B12D(IMAX)
00879          B12D(IMAX) = WORK(IV2TCS+IMAX-1-1)*B12D(IMAX) -
00880      $                WORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1)
00881          B12E(IMAX-1) = TEMP
00882          TEMP = WORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) +
00883      $          WORK(IV2TSN+IMAX-1-1)*B22D(IMAX)
00884          B22D(IMAX) = WORK(IV2TCS+IMAX-1-1)*B22D(IMAX) -
00885      $                WORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1)
00886          B22E(IMAX-1) = TEMP
00887 *
00888 *        Update singular vectors
00889 *
00890          IF( WANTU1 ) THEN
00891             IF( COLMAJOR ) THEN
00892                CALL DLASR( 'R', 'V', 'F', P, IMAX-IMIN+1,
00893      $                     WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1),
00894      $                     U1(1,IMIN), LDU1 )
00895             ELSE
00896                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, P,
00897      $                     WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1),
00898      $                     U1(IMIN,1), LDU1 )
00899             END IF
00900          END IF
00901          IF( WANTU2 ) THEN
00902             IF( COLMAJOR ) THEN
00903                CALL DLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1,
00904      $                     WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1),
00905      $                     U2(1,IMIN), LDU2 )
00906             ELSE
00907                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P,
00908      $                     WORK(IU2CS+IMIN-1), WORK(IU2SN+IMIN-1),
00909      $                     U2(IMIN,1), LDU2 )
00910             END IF
00911          END IF
00912          IF( WANTV1T ) THEN
00913             IF( COLMAJOR ) THEN
00914                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q,
00915      $                     WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1),
00916      $                     V1T(IMIN,1), LDV1T )
00917             ELSE
00918                CALL DLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1,
00919      $                     WORK(IV1TCS+IMIN-1), WORK(IV1TSN+IMIN-1),
00920      $                     V1T(1,IMIN), LDV1T )
00921             END IF
00922          END IF
00923          IF( WANTV2T ) THEN
00924             IF( COLMAJOR ) THEN
00925                CALL DLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q,
00926      $                     WORK(IV2TCS+IMIN-1), WORK(IV2TSN+IMIN-1),
00927      $                     V2T(IMIN,1), LDV2T )
00928             ELSE
00929                CALL DLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1,
00930      $                     WORK(IV2TCS+IMIN-1), WORK(IV2TSN+IMIN-1),
00931      $                     V2T(1,IMIN), LDV2T )
00932             END IF
00933          END IF
00934 *
00935 *        Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
00936 *
00937          IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN
00938             B11D(IMAX) = -B11D(IMAX)
00939             B21D(IMAX) = -B21D(IMAX)
00940             IF( WANTV1T ) THEN
00941                IF( COLMAJOR ) THEN
00942                   CALL DSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
00943                ELSE
00944                   CALL DSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
00945                END IF
00946             END IF
00947          END IF
00948 *
00949 *        Compute THETA(IMAX)
00950 *
00951          X1 = COS(PHI(IMAX-1))*B11D(IMAX) +
00952      $        SIN(PHI(IMAX-1))*B12E(IMAX-1)
00953          Y1 = COS(PHI(IMAX-1))*B21D(IMAX) +
00954      $        SIN(PHI(IMAX-1))*B22E(IMAX-1)
00955 *
00956          THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) )
00957 *
00958 *        Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
00959 *        and B22(IMAX,IMAX-1)
00960 *
00961          IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN
00962             B12D(IMAX) = -B12D(IMAX)
00963             IF( WANTU1 ) THEN
00964                IF( COLMAJOR ) THEN
00965                   CALL DSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
00966                ELSE
00967                   CALL DSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
00968                END IF
00969             END IF
00970          END IF
00971          IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN
00972             B22D(IMAX) = -B22D(IMAX)
00973             IF( WANTU2 ) THEN
00974                IF( COLMAJOR ) THEN
00975                   CALL DSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
00976                ELSE
00977                   CALL DSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
00978                END IF
00979             END IF
00980          END IF
00981 *
00982 *        Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
00983 *
00984          IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
00985             IF( WANTV2T ) THEN
00986                IF( COLMAJOR ) THEN
00987                   CALL DSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
00988                ELSE
00989                   CALL DSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
00990                END IF
00991             END IF
00992          END IF
00993 *
00994 *        Test for negligible sines or cosines
00995 *
00996          DO I = IMIN, IMAX
00997             IF( THETA(I) .LT. THRESH ) THEN
00998                THETA(I) = ZERO
00999             ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
01000                THETA(I) = PIOVER2
01001             END IF
01002          END DO
01003          DO I = IMIN, IMAX-1
01004             IF( PHI(I) .LT. THRESH ) THEN
01005                PHI(I) = ZERO
01006             ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
01007                PHI(I) = PIOVER2
01008             END IF
01009          END DO
01010 *
01011 *        Deflate
01012 *
01013          IF (IMAX .GT. 1) THEN
01014             DO WHILE( PHI(IMAX-1) .EQ. ZERO )
01015                IMAX = IMAX - 1
01016                IF (IMAX .LE. 1) EXIT
01017             END DO
01018          END IF
01019          IF( IMIN .GT. IMAX - 1 )
01020      $      IMIN = IMAX - 1
01021          IF (IMIN .GT. 1) THEN
01022             DO WHILE (PHI(IMIN-1) .NE. ZERO)
01023                 IMIN = IMIN - 1
01024                 IF (IMIN .LE. 1) EXIT
01025             END DO
01026          END IF
01027 *
01028 *        Repeat main iteration loop
01029 *
01030       END DO
01031 *
01032 *     Postprocessing: order THETA from least to greatest
01033 *
01034       DO I = 1, Q
01035 *
01036          MINI = I
01037          THETAMIN = THETA(I)
01038          DO J = I+1, Q
01039             IF( THETA(J) .LT. THETAMIN ) THEN
01040                MINI = J
01041                THETAMIN = THETA(J)
01042             END IF
01043          END DO
01044 *
01045          IF( MINI .NE. I ) THEN
01046             THETA(MINI) = THETA(I)
01047             THETA(I) = THETAMIN
01048             IF( COLMAJOR ) THEN
01049                IF( WANTU1 )
01050      $            CALL DSWAP( P, U1(1,I), 1, U1(1,MINI), 1 )
01051                IF( WANTU2 )
01052      $            CALL DSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 )
01053                IF( WANTV1T )
01054      $            CALL DSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T )
01055                IF( WANTV2T )
01056      $            CALL DSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1),
01057      $               LDV2T )
01058             ELSE
01059                IF( WANTU1 )
01060      $            CALL DSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 )
01061                IF( WANTU2 )
01062      $            CALL DSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 )
01063                IF( WANTV1T )
01064      $            CALL DSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 )
01065                IF( WANTV2T )
01066      $            CALL DSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 )
01067             END IF
01068          END IF
01069 *
01070       END DO
01071 *
01072       RETURN
01073 *
01074 *     End of DBBCSD
01075 *
01076       END
01077 
 All Files Functions