dbdsvdx.f

*> \brief \b DBDSVDX
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DBDSVDX + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*     SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
*    $                    NS, S, Z, LDZ, WORK, IWORK, INFO )
*
*     .. Scalar Arguments ..
*      CHARACTER          JOBZ, RANGE, UPLO
*      INTEGER            IL, INFO, IU, LDZ, N, NS
*      DOUBLE PRECISION   VL, VU
*     ..
*     .. Array Arguments ..
*      INTEGER            IWORK( * )
*      DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),
*                         Z( LDZ, * )
*       ..
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>  DBDSVDX computes the singular value decomposition (SVD) of a real
*>  N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,
*>  where S is a diagonal matrix with non-negative diagonal elements
*>  (the singular values of B), and U and VT are orthogonal matrices
*>  of left and right singular vectors, respectively.
*>
*>  Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]
*>  and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the
*>  singular value decompositon of B through the eigenvalues and
*>  eigenvectors of the N*2-by-N*2 tridiagonal matrix
*>
*>        |  0  d_1                |
*>        | d_1  0  e_1            |
*>  TGK = |     e_1  0  d_2        |
*>        |         d_2  .   .     |
*>        |              .   .   . |
*>
*>  If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then
*>  (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /
*>  sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and
*>  P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].
*>
*>  Given a TGK matrix, one can either a) compute -s,-v and change signs
*>  so that the singular values (and corresponding vectors) are already in
*>  descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder
*>  the values (and corresponding vectors). DBDSVDX implements a) by
*>  calling DSTEVX (bisection plus inverse iteration, to be replaced
*>  with a version of the Multiple Relative Robust Representation
*>  algorithm. (See P. Willems and B. Lang, A framework for the MR^3
*>  algorithm: theory and implementation, SIAM J. Sci. Comput.,
*>  35:740-766, 2013.)
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  B is upper bidiagonal;
*>          = 'L':  B is lower bidiagonal.
*> \endverbatim
*>
*> \param[in] JOBZ
*> \verbatim
*>          JOBZ is CHARACTER*1
*>          = 'N':  Compute singular values only;
*>          = 'V':  Compute singular values and singular vectors.
*> \endverbatim
*>
*> \param[in] RANGE
*> \verbatim
*>          RANGE is CHARACTER*1
*>          = 'A': all singular values will be found.
*>          = 'V': all singular values in the half-open interval [VL,VU)
*>                 will be found.
*>          = 'I': the IL-th through IU-th singular values will be found.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the bidiagonal matrix.  N >= 0.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>          The n diagonal elements of the bidiagonal matrix B.
*> \endverbatim
*>
*> \param[in] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (max(1,N-1))
*>          The (n-1) superdiagonal elements of the bidiagonal matrix
*>          B in elements 1 to N-1.
*> \endverbatim
*>
*> \param[in] VL
*> \verbatim
*>         VL is DOUBLE PRECISION
*>          If RANGE='V', the lower bound of the interval to
*>          be searched for singular values. VU > VL.
*>          Not referenced if RANGE = 'A' or 'I'.
*> \endverbatim
*>
*> \param[in] VU
*> \verbatim
*>         VU is DOUBLE PRECISION
*>          If RANGE='V', the upper bound of the interval to
*>          be searched for singular values. VU > VL.
*>          Not referenced if RANGE = 'A' or 'I'.
*> \endverbatim
*>
*> \param[in] IL
*> \verbatim
*>          IL is INTEGER
*>          If RANGE='I', the index of the
*>          smallest singular value to be returned.
*>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
*>          Not referenced if RANGE = 'A' or 'V'.
*> \endverbatim
*>
*> \param[in] IU
*> \verbatim
*>          IU is INTEGER
*>          If RANGE='I', the index of the
*>          largest singular value to be returned.
*>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
*>          Not referenced if RANGE = 'A' or 'V'.
*> \endverbatim
*>
*> \param[out] NS
*> \verbatim
*>          NS is INTEGER
*>          The total number of singular values found.  0 <= NS <= N.
*>          If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1.
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          S is DOUBLE PRECISION array, dimension (N)
*>          The first NS elements contain the selected singular values in
*>          ascending order.
*> \endverbatim
*>
*> \param[out] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension (2*N,K)
*>          If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
*>          contain the singular vectors of the matrix B corresponding to
*>          the selected singular values, with U in rows 1 to N and V
*>          in rows N+1 to N*2, i.e.
*>          Z = [ U ]
*>              [ V ]
*>          If JOBZ = 'N', then Z is not referenced.
*>          Note: The user must ensure that at least K = NS+1 columns are
*>          supplied in the array Z; if RANGE = 'V', the exact value of
*>          NS is not known in advance and an upper bound must be used.
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>          The leading dimension of the array Z. LDZ >= 1, and if
*>          JOBZ = 'V', LDZ >= max(2,N*2).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (14*N)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (12*N)
*>          If JOBZ = 'V', then if INFO = 0, the first NS elements of
*>          IWORK are zero. If INFO > 0, then IWORK contains the indices
*>          of the eigenvectors that failed to converge in DSTEVX.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*>          > 0:  if INFO = i, then i eigenvectors failed to converge
*>                   in DSTEVX. The indices of the eigenvectors
*>                   (as returned by DSTEVX) are stored in the
*>                   array IWORK.
*>                if INFO = N*2 + 1, an internal error occurred.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup bdsvdx
*
*  =====================================================================
      SUBROUTINE dbdsvdx( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
     $                    NS, S, Z, LDZ, WORK, IWORK, INFO)
*
*  -- LAPACK driver routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ, RANGE, UPLO
      INTEGER            IL, INFO, IU, LDZ, N, NS
      DOUBLE PRECISION   VL, VU
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),
     $                   z( ldz, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TEN, HNDRD, MEIGTH
      parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0,
     $                     hndrd = 100.0d0, meigth = -0.1250d0 )
      DOUBLE PRECISION   FUDGE
      parameter( fudge = 2.0d0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          RNGVX
      LOGICAL            ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
      INTEGER            I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
     $                   ietgk, iifail, iiwork, iltgk, irowu, irowv,
     $                   irowz, isbeg, isplt, itemp, iutgk, j, k,
     $                   ntgk, nru, nrv, nsl
      DOUBLE PRECISION   ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
     $                   smin, sqrt2, thresh, tol, ulp,
     $                   vltgk, vutgk, zjtji
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            IDAMAX
      DOUBLE PRECISION   DDOT, DLAMCH, DNRM2
      EXTERNAL           idamax, lsame, daxpy, ddot, dlamch, dnrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           dstevx, dcopy, dlaset, dscal, dswap, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, sign, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      allsv = lsame( range, 'A' )
      valsv = lsame( range, 'V' )
      indsv = lsame( range, 'I' )
      wantz = lsame( jobz, 'V' )
      lower = lsame( uplo, 'L' )
*
      info = 0
      IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
         info = -1
      ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
         info = -2
      ELSE IF( .NOT.( allsv .OR. valsv .OR. indsv ) ) THEN
         info = -3
      ELSE IF( n.LT.0 ) THEN
         info = -4
      ELSE IF( n.GT.0 ) THEN
         IF( valsv ) THEN
            IF( vl.LT.zero ) THEN
               info = -7
            ELSE IF( vu.LE.vl ) THEN
               info = -8
            END IF
         ELSE IF( indsv ) THEN
            IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
               info = -9
            ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
               info = -10
            END IF
         END IF
      END IF
      IF( info.EQ.0 ) THEN
         IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n*2 ) ) info = -14
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DBDSVDX', -info )
         RETURN
      END IF
*
*     Quick return if possible (N.LE.1)
*
      ns = 0
      IF( n.EQ.0 ) RETURN
*
      IF( n.EQ.1 ) THEN
         IF( allsv .OR. indsv ) THEN
            ns = 1
            s( 1 ) = abs( d( 1 ) )
         ELSE
            IF( vl.LT.abs( d( 1 ) ) .AND. vu.GE.abs( d( 1 ) ) ) THEN
               ns = 1
               s( 1 ) = abs( d( 1 ) )
            END IF
         END IF
         IF( wantz ) THEN
            z( 1, 1 ) = sign( one, d( 1 ) )
            z( 2, 1 ) = one
         END IF
         RETURN
      END IF
*
      abstol = 2*dlamch( 'Safe Minimum' )
      ulp = dlamch( 'Precision' )
      eps = dlamch( 'Epsilon' )
      sqrt2 = sqrt( 2.0d0 )
      ortol = sqrt( ulp )
*
*     Criterion for splitting is taken from DBDSQR when singular
*     values are computed to relative accuracy TOL. (See J. Demmel and
*     W. Kahan, Accurate singular values of bidiagonal matrices, SIAM
*     J. Sci. and Stat. Comput., 11:873–912, 1990.)
*
      tol = max( ten, min( hndrd, eps**meigth ) )*eps
*
*     Compute approximate maximum, minimum singular values.
*
      i = idamax( n, d, 1 )
      smax = abs( d( i ) )
      i = idamax( n-1, e, 1 )
      smax = max( smax, abs( e( i ) ) )
*
*     Compute threshold for neglecting D's and E's.
*
      smin = abs( d( 1 ) )
      IF( smin.NE.zero ) THEN
         mu = smin
         DO i = 2, n
            mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
            smin = min( smin, mu )
            IF( smin.EQ.zero ) EXIT
         END DO
      END IF
      smin = smin / sqrt( dble( n ) )
      thresh = tol*smin
*
*     Check for zeros in D and E (splits), i.e. submatrices.
*
      DO i = 1, n-1
         IF( abs( d( i ) ).LE.thresh ) d( i ) = zero
         IF( abs( e( i ) ).LE.thresh ) e( i ) = zero
      END DO
      IF( abs( d( n ) ).LE.thresh ) d( n ) = zero
*
*     Pointers for arrays used by DSTEVX.
*
      idtgk = 1
      ietgk = idtgk + n*2
      itemp = ietgk + n*2
      iifail = 1
      iiwork = iifail + n*2
*
*     Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode.
*     VL,VU or IL,IU are redefined to conform to implementation a)
*     described in the leading comments.
*
      iltgk = 0
      iutgk = 0
      vltgk = zero
      vutgk = zero
*
      IF( allsv ) THEN
*
*        All singular values will be found. We aim at -s (see
*        leading comments) with RNGVX = 'I'. IL and IU are set
*        later (as ILTGK and IUTGK) according to the dimension
*        of the active submatrix.
*
         rngvx = 'I'
         IF( wantz ) CALL dlaset( 'F', n*2, n+1, zero, zero, z, ldz )
      ELSE IF( valsv ) THEN
*
*        Find singular values in a half-open interval. We aim
*        at -s (see leading comments) and we swap VL and VU
*        (as VUTGK and VLTGK), changing their signs.
*
         rngvx = 'V'
         vltgk = -vu
         vutgk = -vl
         work( idtgk:idtgk+2*n-1 ) = zero
         CALL dcopy( n, d, 1, work( ietgk ), 2 )
         CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
         CALL dstevx( 'N', 'V', n*2, work( idtgk ), work( ietgk ),
     $                vltgk, vutgk, iltgk, iltgk, abstol, ns, s,
     $                z, ldz, work( itemp ), iwork( iiwork ),
     $                iwork( iifail ), info )
         IF( ns.EQ.0 ) THEN
            RETURN
         ELSE
            IF( wantz ) CALL dlaset( 'F', n*2, ns, zero, zero, z, ldz )
         END IF
      ELSE IF( indsv ) THEN
*
*        Find the IL-th through the IU-th singular values. We aim
*        at -s (see leading comments) and indices are mapped into
*        values, therefore mimicking DSTEBZ, where
*
*        GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
*        GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
*
         iltgk = il
         iutgk = iu
         rngvx = 'V'
         work( idtgk:idtgk+2*n-1 ) = zero
         CALL dcopy( n, d, 1, work( ietgk ), 2 )
         CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
         CALL dstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
     $                vltgk, vltgk, iltgk, iltgk, abstol, ns, s,
     $                z, ldz, work( itemp ), iwork( iiwork ),
     $                iwork( iifail ), info )
         vltgk = s( 1 ) - fudge*smax*ulp*n
         work( idtgk:idtgk+2*n-1 ) = zero
         CALL dcopy( n, d, 1, work( ietgk ), 2 )
         CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
         CALL dstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
     $                vutgk, vutgk, iutgk, iutgk, abstol, ns, s,
     $                z, ldz, work( itemp ), iwork( iiwork ),
     $                iwork( iifail ), info )
         vutgk = s( 1 ) + fudge*smax*ulp*n
         vutgk = min( vutgk, zero )
*
*        If VLTGK=VUTGK, DSTEVX returns an error message,
*        so if needed we change VUTGK slightly.
*
         IF( vltgk.EQ.vutgk ) vltgk = vltgk - tol
*
         IF( wantz ) CALL dlaset( 'F', n*2, iu-il+1, zero, zero, z, ldz)
      END IF
*
*     Initialize variables and pointers for S, Z, and WORK.
*
*     NRU, NRV: number of rows in U and V for the active submatrix
*     IDBEG, ISBEG: offsets for the entries of D and S
*     IROWZ, ICOLZ: offsets for the rows and columns of Z
*     IROWU, IROWV: offsets for the rows of U and V
*
      ns = 0
      nru = 0
      nrv = 0
      idbeg = 1
      isbeg = 1
      irowz = 1
      icolz = 1
      irowu = 2
      irowv = 1
      split = .false.
      sveq0 = .false.
*
*     Form the tridiagonal TGK matrix.
*
      s( 1:n ) = zero
      work( ietgk+2*n-1 ) = zero
      work( idtgk:idtgk+2*n-1 ) = zero
      CALL dcopy( n, d, 1, work( ietgk ), 2 )
      CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
*
*
*     Check for splits in two levels, outer level
*     in E and inner level in D.
*
      DO ieptr = 2, n*2, 2
         IF( work( ietgk+ieptr-1 ).EQ.zero ) THEN
*
*           Split in E (this piece of B is square) or bottom
*           of the (input bidiagonal) matrix.
*
            isplt = idbeg
            idend = ieptr - 1
            DO idptr = idbeg, idend, 2
               IF( work( ietgk+idptr-1 ).EQ.zero ) THEN
*
*                 Split in D (rectangular submatrix). Set the number
*                 of rows in U and V (NRU and NRV) accordingly.
*
                  IF( idptr.EQ.idbeg ) THEN
*
*                    D=0 at the top.
*
                     sveq0 = .true.
                     IF( idbeg.EQ.idend) THEN
                        nru = 1
                        nrv = 1
                     END IF
                  ELSE IF( idptr.EQ.idend ) THEN
*
*                    D=0 at the bottom.
*
                     sveq0 = .true.
                     nru = (idend-isplt)/2 + 1
                     nrv = nru
                     IF( isplt.NE.idbeg ) THEN
                        nru = nru + 1
                     END IF
                  ELSE
                     IF( isplt.EQ.idbeg ) THEN
*
*                       Split: top rectangular submatrix.
*
                        nru = (idptr-idbeg)/2
                        nrv = nru + 1
                     ELSE
*
*                       Split: middle square submatrix.
*
                        nru = (idptr-isplt)/2 + 1
                        nrv = nru
                     END IF
                  END IF
               ELSE IF( idptr.EQ.idend ) THEN
*
*                 Last entry of D in the active submatrix.
*
                  IF( isplt.EQ.idbeg ) THEN
*
*                    No split (trivial case).
*
                     nru = (idend-idbeg)/2 + 1
                     nrv = nru
                  ELSE
*
*                    Split: bottom rectangular submatrix.
*
                     nrv = (idend-isplt)/2 + 1
                     nru = nrv + 1
                  END IF
               END IF
*
               ntgk = nru + nrv
*
               IF( ntgk.GT.0 ) THEN
*
*                 Compute eigenvalues/vectors of the active
*                 submatrix according to RANGE:
*                 if RANGE='A' (ALLSV) then RNGVX = 'I'
*                 if RANGE='V' (VALSV) then RNGVX = 'V'
*                 if RANGE='I' (INDSV) then RNGVX = 'V'
*
                  iltgk = 1
                  iutgk = ntgk / 2
                  IF( allsv .OR. vutgk.EQ.zero ) THEN
                     IF( sveq0 .OR.
     $                   smin.LT.eps .OR.
     $                   mod(ntgk,2).GT.0 ) THEN
*                        Special case: eigenvalue equal to zero or very
*                        small, additional eigenvector is needed.
                         iutgk = iutgk + 1
                     END IF
                  END IF
*
*                 Workspace needed by DSTEVX:
*                 WORK( ITEMP: ): 2*5*NTGK
*                 IWORK( 1: ): 2*6*NTGK
*
                  CALL dstevx( jobz, rngvx, ntgk, work( idtgk+isplt-1 ),
     $                         work( ietgk+isplt-1 ), vltgk, vutgk,
     $                         iltgk, iutgk, abstol, nsl, s( isbeg ),
     $                         z( irowz,icolz ), ldz, work( itemp ),
     $                         iwork( iiwork ), iwork( iifail ),
     $                         info )
                  IF( info.NE.0 ) THEN
*                    Exit with the error code from DSTEVX.
                     RETURN
                  END IF
                  emin = abs( maxval( s( isbeg:isbeg+nsl-1 ) ) )
*
                  IF( nsl.GT.0 .AND. wantz ) THEN
*
*                    Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:),
*                    changing the sign of v as discussed in the leading
*                    comments. The norms of u and v may be (slightly)
*                    different from 1/sqrt(2) if the corresponding
*                    eigenvalues are very small or too close. We check
*                    those norms and, if needed, reorthogonalize the
*                    vectors.
*
                     IF( nsl.GT.1 .AND.
     $                   vutgk.EQ.zero .AND.
     $                   mod(ntgk,2).EQ.0 .AND.
     $                   emin.EQ.0 .AND. .NOT.split ) THEN
*
*                       D=0 at the top or bottom of the active submatrix:
*                       one eigenvalue is equal to zero; concatenate the
*                       eigenvectors corresponding to the two smallest
*                       eigenvalues.
*
                        z( irowz:irowz+ntgk-1,icolz+nsl-2 ) =
     $                  z( irowz:irowz+ntgk-1,icolz+nsl-2 ) +
     $                  z( irowz:irowz+ntgk-1,icolz+nsl-1 )
                        z( irowz:irowz+ntgk-1,icolz+nsl-1 ) =
     $                  zero
*                       IF( IUTGK*2.GT.NTGK ) THEN
*                          Eigenvalue equal to zero or very small.
*                          NSL = NSL - 1
*                       END IF
                     END IF
*
                     DO i = 0, min( nsl-1, nru-1 )
                        nrmu = dnrm2( nru, z( irowu, icolz+i ), 2 )
                        IF( nrmu.EQ.zero ) THEN
                           info = n*2 + 1
                           RETURN
                        END IF
                        CALL dscal( nru, one/nrmu,
     $                              z( irowu,icolz+i ), 2 )
                        IF( nrmu.NE.one .AND.
     $                      abs( nrmu-ortol )*sqrt2.GT.one )
     $                      THEN
                           DO j = 0, i-1
                              zjtji = -ddot( nru, z( irowu, icolz+j ),
     $                                       2, z( irowu, icolz+i ), 2 )
                              CALL daxpy( nru, zjtji,
     $                                    z( irowu, icolz+j ), 2,
     $                                    z( irowu, icolz+i ), 2 )
                           END DO
                           nrmu = dnrm2( nru, z( irowu, icolz+i ), 2 )
                           CALL dscal( nru, one/nrmu,
     $                                 z( irowu,icolz+i ), 2 )
                        END IF
                     END DO
                     DO i = 0, min( nsl-1, nrv-1 )
                        nrmv = dnrm2( nrv, z( irowv, icolz+i ), 2 )
                        IF( nrmv.EQ.zero ) THEN
                           info = n*2 + 1
                           RETURN
                        END IF
                        CALL dscal( nrv, -one/nrmv,
     $                              z( irowv,icolz+i ), 2 )
                        IF( nrmv.NE.one .AND.
     $                      abs( nrmv-ortol )*sqrt2.GT.one )
     $                      THEN
                           DO j = 0, i-1
                              zjtji = -ddot( nrv, z( irowv, icolz+j ),
     $                                       2, z( irowv, icolz+i ), 2 )
                              CALL daxpy( nru, zjtji,
     $                                    z( irowv, icolz+j ), 2,
     $                                    z( irowv, icolz+i ), 2 )
                           END DO
                           nrmv = dnrm2( nrv, z( irowv, icolz+i ), 2 )
                           CALL dscal( nrv, one/nrmv,
     $                                 z( irowv,icolz+i ), 2 )
                        END IF
                     END DO
                     IF( vutgk.EQ.zero .AND.
     $                   idptr.LT.idend .AND.
     $                   mod(ntgk,2).GT.0 ) THEN
*
*                       D=0 in the middle of the active submatrix (one
*                       eigenvalue is equal to zero): save the corresponding
*                       eigenvector for later use (when bottom of the
*                       active submatrix is reached).
*
                        split = .true.
                        z( irowz:irowz+ntgk-1,n+1 ) =
     $                     z( irowz:irowz+ntgk-1,ns+nsl )
                        z( irowz:irowz+ntgk-1,ns+nsl ) =
     $                     zero
                     END IF
                  END IF !** WANTZ **!
*
                  nsl = min( nsl, nru )
                  sveq0 = .false.
*
*                 Absolute values of the eigenvalues of TGK.
*
                  DO i = 0, nsl-1
                     s( isbeg+i ) = abs( s( isbeg+i ) )
                  END DO
*
*                 Update pointers for TGK, S and Z.
*
                  isbeg = isbeg + nsl
                  irowz = irowz + ntgk
                  icolz = icolz + nsl
                  irowu = irowz
                  irowv = irowz + 1
                  isplt = idptr + 1
                  ns = ns + nsl
                  nru = 0
                  nrv = 0
               END IF !** NTGK.GT.0 **!
               IF( irowz.LT.n*2 .AND. wantz ) THEN
                  z( 1:irowz-1, icolz ) = zero
               END IF
            END DO !** IDPTR loop **!
            IF( split .AND. wantz ) THEN
*
*              Bring back eigenvector corresponding
*              to eigenvalue equal to zero.
*
               z( idbeg:idend-ntgk+1,isbeg-1 ) =
     $         z( idbeg:idend-ntgk+1,isbeg-1 ) +
     $         z( idbeg:idend-ntgk+1,n+1 )
               z( idbeg:idend-ntgk+1,n+1 ) = 0
            END IF
            irowv = irowv - 1
            irowu = irowu + 1
            idbeg = ieptr + 1
            sveq0 = .false.
            split = .false.
         END IF !** Check for split in E **!
      END DO !** IEPTR loop **!
*
*     Sort the singular values into decreasing order (insertion sort on
*     singular values, but only one transposition per singular vector)
*
      DO i = 1, ns-1
         k = 1
         smin = s( 1 )
         DO j = 2, ns + 1 - i
            IF( s( j ).LE.smin ) THEN
               k = j
               smin = s( j )
            END IF
         END DO
         IF( k.NE.ns+1-i ) THEN
            s( k ) = s( ns+1-i )
            s( ns+1-i ) = smin
            IF( wantz ) CALL dswap( n*2, z( 1,k ), 1, z( 1,ns+1-i ), 1 )
         END IF
      END DO
*
*     If RANGE=I, check for singular values/vectors to be discarded.
*
      IF( indsv ) THEN
         k = iu - il + 1
         IF( k.LT.ns ) THEN
            s( k+1:ns ) = zero
            IF( wantz ) z( 1:n*2,k+1:ns ) = zero
            ns = k
         END IF
      END IF
*
*     Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).
*     If B is a lower diagonal, swap U and V.
*
      IF( wantz ) THEN
      DO i = 1, ns
         CALL dcopy( n*2, z( 1,i ), 1, work, 1 )
         IF( lower ) THEN
            CALL dcopy( n, work( 2 ), 2, z( n+1,i ), 1 )
            CALL dcopy( n, work( 1 ), 2, z( 1  ,i ), 1 )
         ELSE
            CALL dcopy( n, work( 2 ), 2, z( 1  ,i ), 1 )
            CALL dcopy( n, work( 1 ), 2, z( n+1,i ), 1 )
         END IF
      END DO
      END IF
*
      RETURN
*
*     End of DBDSVDX
*
      END