SRC/dbdsvdx.f

   1 *> \brief \b DBDSVDX
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DBDSVDX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *     SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
  22 *    $                    NS, S, Z, LDZ, WORK, IWORK, INFO )
  23 *
  24 *     .. Scalar Arguments ..
  25 *      CHARACTER          JOBZ, RANGE, UPLO
  26 *      INTEGER            IL, INFO, IU, LDZ, N, NS
  27 *      DOUBLE PRECISION   VL, VU
  28 *     ..
  29 *     .. Array Arguments ..
  30 *      INTEGER            IWORK( * )
  31 *      DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),
  32 *                         Z( LDZ, * )
  33 *       ..
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *>  DBDSVDX computes the singular value decomposition (SVD) of a real
  41 *>  N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,
  42 *>  where S is a diagonal matrix with non-negative diagonal elements
  43 *>  (the singular values of B), and U and VT are orthogonal matrices
  44 *>  of left and right singular vectors, respectively.
  45 *>
  46 *>  Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]
  47 *>  and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the
  48 *>  singular value decompositon of B through the eigenvalues and
  49 *>  eigenvectors of the N*2-by-N*2 tridiagonal matrix
  50 *>
  51 *>        |  0  d_1                |
  52 *>        | d_1  0  e_1            |
  53 *>  TGK = |     e_1  0  d_2        |
  54 *>        |         d_2  .   .     |
  55 *>        |              .   .   . |
  56 *>
  57 *>  If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then
  58 *>  (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /
  59 *>  sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and
  60 *>  P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].
  61 *>
  62 *>  Given a TGK matrix, one can either a) compute -s,-v and change signs
  63 *>  so that the singular values (and corresponding vectors) are already in
  64 *>  descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder
  65 *>  the values (and corresponding vectors). DBDSVDX implements a) by
  66 *>  calling DSTEVX (bisection plus inverse iteration, to be replaced
  67 *>  with a version of the Multiple Relative Robust Representation
  68 *>  algorithm. (See P. Willems and B. Lang, A framework for the MR^3
  69 *>  algorithm: theory and implementation, SIAM J. Sci. Comput.,
  70 *>  35:740-766, 2013.)
  71 *> \endverbatim
  72 *
  73 *  Arguments:
  74 *  ==========
  75 *
  76 *> \param[in] UPLO
  77 *> \verbatim
  78 *>          UPLO is CHARACTER*1
  79 *>          = 'U':  B is upper bidiagonal;
  80 *>          = 'L':  B is lower bidiagonal.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] JOBZ
  84 *> \verbatim
  85 *>          JOBZ is CHARACTER*1
  86 *>          = 'N':  Compute singular values only;
  87 *>          = 'V':  Compute singular values and singular vectors.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] RANGE
  91 *> \verbatim
  92 *>          RANGE is CHARACTER*1
  93 *>          = 'A': all singular values will be found.
  94 *>          = 'V': all singular values in the half-open interval [VL,VU)
  95 *>                 will be found.
  96 *>          = 'I': the IL-th through IU-th singular values will be found.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] N
 100 *> \verbatim
 101 *>          N is INTEGER
 102 *>          The order of the bidiagonal matrix.  N >= 0.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] D
 106 *> \verbatim
 107 *>          D is DOUBLE PRECISION array, dimension (N)
 108 *>          The n diagonal elements of the bidiagonal matrix B.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] E
 112 *> \verbatim
 113 *>          E is DOUBLE PRECISION array, dimension (max(1,N-1))
 114 *>          The (n-1) superdiagonal elements of the bidiagonal matrix
 115 *>          B in elements 1 to N-1.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] VL
 119 *> \verbatim
 120 *>         VL is DOUBLE PRECISION
 121 *>          If RANGE='V', the lower bound of the interval to
 122 *>          be searched for singular values. VU > VL.
 123 *>          Not referenced if RANGE = 'A' or 'I'.
 124 *> \endverbatim
 125 *>
 126 *> \param[in] VU
 127 *> \verbatim
 128 *>         VU is DOUBLE PRECISION
 129 *>          If RANGE='V', the upper bound of the interval to
 130 *>          be searched for singular values. VU > VL.
 131 *>          Not referenced if RANGE = 'A' or 'I'.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] IL
 135 *> \verbatim
 136 *>          IL is INTEGER
 137 *>          If RANGE='I', the index of the
 138 *>          smallest singular value to be returned.
 139 *>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
 140 *>          Not referenced if RANGE = 'A' or 'V'.
 141 *> \endverbatim
 142 *>
 143 *> \param[in] IU
 144 *> \verbatim
 145 *>          IU is INTEGER
 146 *>          If RANGE='I', the index of the
 147 *>          largest singular value to be returned.
 148 *>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
 149 *>          Not referenced if RANGE = 'A' or 'V'.
 150 *> \endverbatim
 151 *>
 152 *> \param[out] NS
 153 *> \verbatim
 154 *>          NS is INTEGER
 155 *>          The total number of singular values found.  0 <= NS <= N.
 156 *>          If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1.
 157 *> \endverbatim
 158 *>
 159 *> \param[out] S
 160 *> \verbatim
 161 *>          S is DOUBLE PRECISION array, dimension (N)
 162 *>          The first NS elements contain the selected singular values in
 163 *>          ascending order.
 164 *> \endverbatim
 165 *>
 166 *> \param[out] Z
 167 *> \verbatim
 168 *>          Z is DOUBLE PRECISION array, dimension (2*N,K) )
 169 *>          If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
 170 *>          contain the singular vectors of the matrix B corresponding to
 171 *>          the selected singular values, with U in rows 1 to N and V
 172 *>          in rows N+1 to N*2, i.e.
 173 *>          Z = [ U ]
 174 *>              [ V ]
 175 *>          If JOBZ = 'N', then Z is not referenced.
 176 *>          Note: The user must ensure that at least K = NS+1 columns are
 177 *>          supplied in the array Z; if RANGE = 'V', the exact value of
 178 *>          NS is not known in advance and an upper bound must be used.
 179 *> \endverbatim
 180 *>
 181 *> \param[in] LDZ
 182 *> \verbatim
 183 *>          LDZ is INTEGER
 184 *>          The leading dimension of the array Z. LDZ >= 1, and if
 185 *>          JOBZ = 'V', LDZ >= max(2,N*2).
 186 *> \endverbatim
 187 *>
 188 *> \param[out] WORK
 189 *> \verbatim
 190 *>          WORK is DOUBLE PRECISION array, dimension (14*N)
 191 *> \endverbatim
 192 *>
 193 *> \param[out] IWORK
 194 *> \verbatim
 195 *>          IWORK is INTEGER array, dimension (12*N)
 196 *>          If JOBZ = 'V', then if INFO = 0, the first NS elements of
 197 *>          IWORK are zero. If INFO > 0, then IWORK contains the indices
 198 *>          of the eigenvectors that failed to converge in DSTEVX.
 199 *> \endverbatim
 200 *>
 201 *> \param[out] INFO
 202 *> \verbatim
 203 *>          INFO is INTEGER
 204 *>          = 0:  successful exit
 205 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 206 *>          > 0:  if INFO = i, then i eigenvectors failed to converge
 207 *>                   in DSTEVX. The indices of the eigenvectors
 208 *>                   (as returned by DSTEVX) are stored in the
 209 *>                   array IWORK.
 210 *>                if INFO = N*2 + 1, an internal error occurred.
 211 *> \endverbatim
 212 *
 213 *  Authors:
 214 *  ========
 215 *
 216 *> \author Univ. of Tennessee
 217 *> \author Univ. of California Berkeley
 218 *> \author Univ. of Colorado Denver
 219 *> \author NAG Ltd.
 220 *
 221 *> \date June 2016
 222 *
 223 *> \ingroup doubleOTHEReigen
 224 *
 225 *  =====================================================================
 226       SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 227      $                    NS, S, Z, LDZ, WORK, IWORK, INFO)
 228 *
 229 *  -- LAPACK driver routine (version 3.6.1) --
 230 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 231 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 232 *     November 2016
 233 *
 234 *     .. Scalar Arguments ..
 235       CHARACTER          JOBZ, RANGE, UPLO
 236       INTEGER            IL, INFO, IU, LDZ, N, NS
 237       DOUBLE PRECISION   VL, VU
 238 *     ..
 239 *     .. Array Arguments ..
 240       INTEGER            IWORK( * )
 241       DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),
 242      $                   Z( LDZ, * )
 243 *     ..
 244 *
 245 *  =====================================================================
 246 *
 247 *     .. Parameters ..
 248       DOUBLE PRECISION   ZERO, ONE, TEN, HNDRD, MEIGTH
 249       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0,
 250      $                     HNDRD = 100.0D0, MEIGTH = -0.1250D0 )
 251       DOUBLE PRECISION   FUDGE
 252       PARAMETER          ( FUDGE = 2.0D0 )
 253 *     ..
 254 *     .. Local Scalars ..
 255       CHARACTER          RNGVX
 256       LOGICAL            ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
 257       INTEGER            I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
 258      $                   IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV,
 259      $                   IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K,
 260      $                   NTGK, NRU, NRV, NSL
 261       DOUBLE PRECISION   ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
 262      $                   SMIN, SQRT2, THRESH, TOL, ULP,
 263      $                   VLTGK, VUTGK, ZJTJI
 264 *     ..
 265 *     .. External Functions ..
 266       LOGICAL            LSAME
 267       INTEGER            IDAMAX
 268       DOUBLE PRECISION   DDOT, DLAMCH, DNRM2
 269       EXTERNAL           IDAMAX, LSAME, DAXPY, DDOT, DLAMCH, DNRM2
 270 *     ..
 271 *     .. External Subroutines ..
 272       EXTERNAL           DCOPY, DLASET, DSCAL, DSWAP
 273 *     ..
 274 *     .. Intrinsic Functions ..
 275       INTRINSIC          ABS, DBLE, SIGN, SQRT
 276 *     ..
 277 *     .. Executable Statements ..
 278 *
 279 *     Test the input parameters.
 280 *
 281       ALLSV = LSAME( RANGE, 'A' )
 282       VALSV = LSAME( RANGE, 'V' )
 283       INDSV = LSAME( RANGE, 'I' )
 284       WANTZ = LSAME( JOBZ, 'V' )
 285       LOWER = LSAME( UPLO, 'L' )
 286 *
 287       INFO = 0
 288       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
 289          INFO = -1
 290       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 291          INFO = -2
 292       ELSE IF( .NOT.( ALLSV .OR. VALSV .OR. INDSV ) ) THEN
 293          INFO = -3
 294       ELSE IF( N.LT.0 ) THEN
 295          INFO = -4
 296       ELSE IF( N.GT.0 ) THEN
 297          IF( VALSV ) THEN
 298             IF( VL.LT.ZERO ) THEN
 299                INFO = -7
 300             ELSE IF( VU.LE.VL ) THEN
 301                INFO = -8
 302             END IF
 303          ELSE IF( INDSV ) THEN
 304             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
 305                INFO = -9
 306             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
 307                INFO = -10
 308             END IF
 309          END IF
 310       END IF
 311       IF( INFO.EQ.0 ) THEN
 312          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N*2 ) ) INFO = -14
 313       END IF
 314 *
 315       IF( INFO.NE.0 ) THEN
 316          CALL XERBLA( 'DBDSVDX', -INFO )
 317          RETURN
 318       END IF
 319 *
 320 *     Quick return if possible (N.LE.1)
 321 *
 322       NS = 0
 323       IF( N.EQ.0 ) RETURN
 324 *
 325       IF( N.EQ.1 ) THEN
 326          IF( ALLSV .OR. INDSV ) THEN
 327             NS = 1
 328             S( 1 ) = ABS( D( 1 ) )
 329          ELSE
 330             IF( VL.LT.ABS( D( 1 ) ) .AND. VU.GE.ABS( D( 1 ) ) ) THEN
 331                NS = 1
 332                S( 1 ) = ABS( D( 1 ) )
 333             END IF
 334          END IF
 335          IF( WANTZ ) THEN
 336             Z( 1, 1 ) = SIGN( ONE, D( 1 ) )
 337             Z( 2, 1 ) = ONE
 338          END IF
 339          RETURN
 340       END IF
 341 *
 342       ABSTOL = 2*DLAMCH( 'Safe Minimum' )
 343       ULP = DLAMCH( 'Precision' )
 344       EPS = DLAMCH( 'Epsilon' )
 345       SQRT2 = SQRT( 2.0D0 )
 346       ORTOL = SQRT( ULP )
 347 *
 348 *     Criterion for splitting is taken from DBDSQR when singular
 349 *     values are computed to relative accuracy TOL. (See J. Demmel and
 350 *     W. Kahan, Accurate singular values of bidiagonal matrices, SIAM
 351 *     J. Sci. and Stat. Comput., 11:873–912, 1990.)
 352 *
 353       TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS
 354 *
 355 *     Compute approximate maximum, minimum singular values.
 356 *
 357       I = IDAMAX( N, D, 1 )
 358       SMAX = ABS( D( I ) )
 359       I = IDAMAX( N-1, E, 1 )
 360       SMAX = MAX( SMAX, ABS( E( I ) ) )
 361 *
 362 *     Compute threshold for neglecting D's and E's.
 363 *
 364       SMIN = ABS( D( 1 ) )
 365       IF( SMIN.NE.ZERO ) THEN
 366          MU = SMIN
 367          DO I = 2, N
 368             MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
 369             SMIN = MIN( SMIN, MU )
 370             IF( SMIN.EQ.ZERO ) EXIT
 371          END DO
 372       END IF
 373       SMIN = SMIN / SQRT( DBLE( N ) )
 374       THRESH = TOL*SMIN
 375 *
 376 *     Check for zeros in D and E (splits), i.e. submatrices.
 377 *
 378       DO I = 1, N-1
 379          IF( ABS( D( I ) ).LE.THRESH ) D( I ) = ZERO
 380          IF( ABS( E( I ) ).LE.THRESH ) E( I ) = ZERO
 381       END DO
 382       IF( ABS( D( N ) ).LE.THRESH ) D( N ) = ZERO
 383 *
 384 *     Pointers for arrays used by DSTEVX.
 385 *
 386       IDTGK = 1
 387       IETGK = IDTGK + N*2
 388       ITEMP = IETGK + N*2
 389       IIFAIL = 1
 390       IIWORK = IIFAIL + N*2
 391 *
 392 *     Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode.
 393 *     VL,VU or IL,IU are redefined to conform to implementation a)
 394 *     described in the leading comments.
 395 *
 396       ILTGK = 0
 397       IUTGK = 0
 398       VLTGK = ZERO
 399       VUTGK = ZERO
 400 *
 401       IF( ALLSV ) THEN
 402 *
 403 *        All singular values will be found. We aim at -s (see
 404 *        leading comments) with RNGVX = 'I'. IL and IU are set
 405 *        later (as ILTGK and IUTGK) according to the dimension
 406 *        of the active submatrix.
 407 *
 408          RNGVX = 'I'
 409          IF( WANTZ ) CALL DLASET( 'F', N*2, N+1, ZERO, ZERO, Z, LDZ )
 410       ELSE IF( VALSV ) THEN
 411 *
 412 *        Find singular values in a half-open interval. We aim
 413 *        at -s (see leading comments) and we swap VL and VU
 414 *        (as VUTGK and VLTGK), changing their signs.
 415 *
 416          RNGVX = 'V'
 417          VLTGK = -VU
 418          VUTGK = -VL
 419          WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
 420          CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
 421          CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
 422          CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ),
 423      $                VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S,
 424      $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
 425      $                IWORK( IIFAIL ), INFO )
 426          IF( NS.EQ.0 ) THEN
 427             RETURN
 428          ELSE
 429             IF( WANTZ ) CALL DLASET( 'F', N*2, NS, ZERO, ZERO, Z, LDZ )
 430          END IF
 431       ELSE IF( INDSV ) THEN
 432 *
 433 *        Find the IL-th through the IU-th singular values. We aim
 434 *        at -s (see leading comments) and indices are mapped into
 435 *        values, therefore mimicking DSTEBZ, where
 436 *
 437 *        GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
 438 *        GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
 439 *
 440          ILTGK = IL
 441          IUTGK = IU
 442          RNGVX = 'V'
 443          WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
 444          CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
 445          CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
 446          CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),
 447      $                VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S,
 448      $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
 449      $                IWORK( IIFAIL ), INFO )
 450          VLTGK = S( 1 ) - FUDGE*SMAX*ULP*N
 451          WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
 452          CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
 453          CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
 454          CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),
 455      $                VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S,
 456      $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
 457      $                IWORK( IIFAIL ), INFO )
 458          VUTGK = S( 1 ) + FUDGE*SMAX*ULP*N
 459          VUTGK = MIN( VUTGK, ZERO )
 460 *
 461 *        If VLTGK=VUTGK, DSTEVX returns an error message,
 462 *        so if needed we change VUTGK slightly.
 463 *
 464          IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL
 465 *
 466          IF( WANTZ ) CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ)
 467       END IF
 468 *
 469 *     Initialize variables and pointers for S, Z, and WORK.
 470 *
 471 *     NRU, NRV: number of rows in U and V for the active submatrix
 472 *     IDBEG, ISBEG: offsets for the entries of D and S
 473 *     IROWZ, ICOLZ: offsets for the rows and columns of Z
 474 *     IROWU, IROWV: offsets for the rows of U and V
 475 *
 476       NS = 0
 477       NRU = 0
 478       NRV = 0
 479       IDBEG = 1
 480       ISBEG = 1
 481       IROWZ = 1
 482       ICOLZ = 1
 483       IROWU = 2
 484       IROWV = 1
 485       SPLIT = .FALSE.
 486       SVEQ0 = .FALSE.
 487 *
 488 *     Form the tridiagonal TGK matrix.
 489 *
 490       S( 1:N ) = ZERO
 491       WORK( IETGK+2*N-1 ) = ZERO
 492       WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
 493       CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
 494       CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
 495 *
 496 *
 497 *     Check for splits in two levels, outer level
 498 *     in E and inner level in D.
 499 *
 500       DO IEPTR = 2, N*2, 2
 501          IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN
 502 *
 503 *           Split in E (this piece of B is square) or bottom
 504 *           of the (input bidiagonal) matrix.
 505 *
 506             ISPLT = IDBEG
 507             IDEND = IEPTR - 1
 508             DO IDPTR = IDBEG, IDEND, 2
 509                IF( WORK( IETGK+IDPTR-1 ).EQ.ZERO ) THEN
 510 *
 511 *                 Split in D (rectangular submatrix). Set the number
 512 *                 of rows in U and V (NRU and NRV) accordingly.
 513 *
 514                   IF( IDPTR.EQ.IDBEG ) THEN
 515 *
 516 *                    D=0 at the top.
 517 *
 518                      SVEQ0 = .TRUE.
 519                      IF( IDBEG.EQ.IDEND) THEN
 520                         NRU = 1
 521                         NRV = 1
 522                      END IF
 523                   ELSE IF( IDPTR.EQ.IDEND ) THEN
 524 *
 525 *                    D=0 at the bottom.
 526 *
 527                      SVEQ0 = .TRUE.
 528                      NRU = (IDEND-ISPLT)/2 + 1
 529                      NRV = NRU
 530                      IF( ISPLT.NE.IDBEG ) THEN
 531                         NRU = NRU + 1
 532                      END IF
 533                   ELSE
 534                      IF( ISPLT.EQ.IDBEG ) THEN
 535 *
 536 *                       Split: top rectangular submatrix.
 537 *
 538                         NRU = (IDPTR-IDBEG)/2
 539                         NRV = NRU + 1
 540                      ELSE
 541 *
 542 *                       Split: middle square submatrix.
 543 *
 544                         NRU = (IDPTR-ISPLT)/2 + 1
 545                         NRV = NRU
 546                      END IF
 547                   END IF
 548                ELSE IF( IDPTR.EQ.IDEND ) THEN
 549 *
 550 *                 Last entry of D in the active submatrix.
 551 *
 552                   IF( ISPLT.EQ.IDBEG ) THEN
 553 *
 554 *                    No split (trivial case).
 555 *
 556                      NRU = (IDEND-IDBEG)/2 + 1
 557                      NRV = NRU
 558                   ELSE
 559 *
 560 *                    Split: bottom rectangular submatrix.
 561 *
 562                      NRV = (IDEND-ISPLT)/2 + 1
 563                      NRU = NRV + 1
 564                   END IF
 565                END IF
 566 *
 567                NTGK = NRU + NRV
 568 *
 569                IF( NTGK.GT.0 ) THEN
 570 *
 571 *                 Compute eigenvalues/vectors of the active
 572 *                 submatrix according to RANGE:
 573 *                 if RANGE='A' (ALLSV) then RNGVX = 'I'
 574 *                 if RANGE='V' (VALSV) then RNGVX = 'V'
 575 *                 if RANGE='I' (INDSV) then RNGVX = 'V'
 576 *
 577                   ILTGK = 1
 578                   IUTGK = NTGK / 2
 579                   IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN
 580                      IF( SVEQ0 .OR.
 581      $                   SMIN.LT.EPS .OR.
 582      $                   MOD(NTGK,2).GT.0 ) THEN
 583 *                        Special case: eigenvalue equal to zero or very
 584 *                        small, additional eigenvector is needed.
 585                          IUTGK = IUTGK + 1
 586                      END IF
 587                   END IF
 588 *
 589 *                 Workspace needed by DSTEVX:
 590 *                 WORK( ITEMP: ): 2*5*NTGK
 591 *                 IWORK( 1: ): 2*6*NTGK
 592 *
 593                   CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ),
 594      $                         WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK,
 595      $                         ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ),
 596      $                         Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ),
 597      $                         IWORK( IIWORK ), IWORK( IIFAIL ),
 598      $                         INFO )
 599                   IF( INFO.NE.0 ) THEN
 600 *                    Exit with the error code from DSTEVX.
 601                      RETURN
 602                   END IF
 603                   EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) )
 604 *
 605                   IF( NSL.GT.0 .AND. WANTZ ) THEN
 606 *
 607 *                    Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:),
 608 *                    changing the sign of v as discussed in the leading
 609 *                    comments. The norms of u and v may be (slightly)
 610 *                    different from 1/sqrt(2) if the corresponding
 611 *                    eigenvalues are very small or too close. We check
 612 *                    those norms and, if needed, reorthogonalize the
 613 *                    vectors.
 614 *
 615                      IF( NSL.GT.1 .AND.
 616      $                   VUTGK.EQ.ZERO .AND.
 617      $                   MOD(NTGK,2).EQ.0 .AND.
 618      $                   EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN
 619 *
 620 *                       D=0 at the top or bottom of the active submatrix:
 621 *                       one eigenvalue is equal to zero; concatenate the
 622 *                       eigenvectors corresponding to the two smallest
 623 *                       eigenvalues.
 624 *
 625                         Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) =
 626      $                  Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) +
 627      $                  Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 )
 628                         Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) =
 629      $                  ZERO
 630 *                       IF( IUTGK*2.GT.NTGK ) THEN
 631 *                          Eigenvalue equal to zero or very small.
 632 *                          NSL = NSL - 1
 633 *                       END IF
 634                      END IF
 635 *
 636                      DO I = 0, MIN( NSL-1, NRU-1 )
 637                         NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 )
 638                         IF( NRMU.EQ.ZERO ) THEN
 639                            INFO = N*2 + 1
 640                            RETURN
 641                         END IF
 642                         CALL DSCAL( NRU, ONE/NRMU,
 643      $                              Z( IROWU,ICOLZ+I ), 2 )
 644                         IF( NRMU.NE.ONE .AND.
 645      $                      ABS( NRMU-ORTOL )*SQRT2.GT.ONE )
 646      $                      THEN
 647                            DO J = 0, I-1
 648                               ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ),
 649      $                                       2, Z( IROWU, ICOLZ+I ), 2 )
 650                               CALL DAXPY( NRU, ZJTJI,
 651      $                                    Z( IROWU, ICOLZ+J ), 2,
 652      $                                    Z( IROWU, ICOLZ+I ), 2 )
 653                            END DO
 654                            NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 )
 655                            CALL DSCAL( NRU, ONE/NRMU,
 656      $                                 Z( IROWU,ICOLZ+I ), 2 )
 657                         END IF
 658                      END DO
 659                      DO I = 0, MIN( NSL-1, NRV-1 )
 660                         NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 )
 661                         IF( NRMV.EQ.ZERO ) THEN
 662                            INFO = N*2 + 1
 663                            RETURN
 664                         END IF
 665                         CALL DSCAL( NRV, -ONE/NRMV,
 666      $                              Z( IROWV,ICOLZ+I ), 2 )
 667                         IF( NRMV.NE.ONE .AND.
 668      $                      ABS( NRMV-ORTOL )*SQRT2.GT.ONE )
 669      $                      THEN
 670                            DO J = 0, I-1
 671                               ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ),
 672      $                                       2, Z( IROWV, ICOLZ+I ), 2 )
 673                               CALL DAXPY( NRU, ZJTJI,
 674      $                                    Z( IROWV, ICOLZ+J ), 2,
 675      $                                    Z( IROWV, ICOLZ+I ), 2 )
 676                            END DO
 677                            NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 )
 678                            CALL DSCAL( NRV, ONE/NRMV,
 679      $                                 Z( IROWV,ICOLZ+I ), 2 )
 680                         END IF
 681                      END DO
 682                      IF( VUTGK.EQ.ZERO .AND.
 683      $                   IDPTR.LT.IDEND .AND.
 684      $                   MOD(NTGK,2).GT.0 ) THEN
 685 *
 686 *                       D=0 in the middle of the active submatrix (one
 687 *                       eigenvalue is equal to zero): save the corresponding
 688 *                       eigenvector for later use (when bottom of the
 689 *                       active submatrix is reached).
 690 *
 691                         SPLIT = .TRUE.
 692                         Z( IROWZ:IROWZ+NTGK-1,N+1 ) =
 693      $                     Z( IROWZ:IROWZ+NTGK-1,NS+NSL )
 694                         Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) =
 695      $                     ZERO
 696                      END IF
 697                   END IF !** WANTZ **!
 698 *
 699                   NSL = MIN( NSL, NRU )
 700                   SVEQ0 = .FALSE.
 701 *
 702 *                 Absolute values of the eigenvalues of TGK.
 703 *
 704                   DO I = 0, NSL-1
 705                      S( ISBEG+I ) = ABS( S( ISBEG+I ) )
 706                   END DO
 707 *
 708 *                 Update pointers for TGK, S and Z.
 709 *
 710                   ISBEG = ISBEG + NSL
 711                   IROWZ = IROWZ + NTGK
 712                   ICOLZ = ICOLZ + NSL
 713                   IROWU = IROWZ
 714                   IROWV = IROWZ + 1
 715                   ISPLT = IDPTR + 1
 716                   NS = NS + NSL
 717                   NRU = 0
 718                   NRV = 0
 719                END IF !** NTGK.GT.0 **!
 720                IF( IROWZ.LT.N*2 .AND. WANTZ ) THEN
 721                   Z( 1:IROWZ-1, ICOLZ ) = ZERO
 722                END IF
 723             END DO !** IDPTR loop **!
 724             IF( SPLIT .AND. WANTZ ) THEN
 725 *
 726 *              Bring back eigenvector corresponding
 727 *              to eigenvalue equal to zero.
 728 *
 729                Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) =
 730      $         Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) +
 731      $         Z( IDBEG:IDEND-NTGK+1,N+1 )
 732                Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0
 733             END IF
 734             IROWV = IROWV - 1
 735             IROWU = IROWU + 1
 736             IDBEG = IEPTR + 1
 737             SVEQ0 = .FALSE.
 738             SPLIT = .FALSE.
 739          END IF !** Check for split in E **!
 740       END DO !** IEPTR loop **!
 741 *
 742 *     Sort the singular values into decreasing order (insertion sort on
 743 *     singular values, but only one transposition per singular vector)
 744 *
 745       DO I = 1, NS-1
 746          K = 1
 747          SMIN = S( 1 )
 748          DO J = 2, NS + 1 - I
 749             IF( S( J ).LE.SMIN ) THEN
 750                K = J
 751                SMIN = S( J )
 752             END IF
 753          END DO
 754          IF( K.NE.NS+1-I ) THEN
 755             S( K ) = S( NS+1-I )
 756             S( NS+1-I ) = SMIN
 757             IF( WANTZ ) CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 )
 758          END IF
 759       END DO
 760 *
 761 *     If RANGE=I, check for singular values/vectors to be discarded.
 762 *
 763       IF( INDSV ) THEN
 764          K = IU - IL + 1
 765          IF( K.LT.NS ) THEN
 766             S( K+1:NS ) = ZERO
 767             IF( WANTZ ) Z( 1:N*2,K+1:NS ) = ZERO
 768             NS = K
 769          END IF
 770       END IF
 771 *
 772 *     Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).
 773 *     If B is a lower diagonal, swap U and V.
 774 *
 775       IF( WANTZ ) THEN
 776       DO I = 1, NS
 777          CALL DCOPY( N*2, Z( 1,I ), 1, WORK, 1 )
 778          IF( LOWER ) THEN
 779             CALL DCOPY( N, WORK( 2 ), 2, Z( N+1,I ), 1 )
 780             CALL DCOPY( N, WORK( 1 ), 2, Z( 1  ,I ), 1 )
 781          ELSE
 782             CALL DCOPY( N, WORK( 2 ), 2, Z( 1  ,I ), 1 )
 783             CALL DCOPY( N, WORK( 1 ), 2, Z( N+1,I ), 1 )
 784          END IF
 785       END DO
 786       END IF
 787 *
 788       RETURN
 789 *
 790 *     End of DBDSVDX
 791 *
 792       END