dbdsvdx()

subroutine dbdsvdx	(	character	uplo,
		character	jobz,
		character	range,
		integer	n,
		double precision, dimension( * )	d,
		double precision, dimension( * )	e,
		double precision	vl,
		double precision	vu,
		integer	il,
		integer	iu,
		integer	ns,
		double precision, dimension( * )	s,
		double precision, dimension( ldz, * )	z,
		integer	ldz,
		double precision, dimension( * )	work,
		integer, dimension( * )	iwork,
		integer	info )

DBDSVDX

Download DBDSVDX + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!>  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 decomposition 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.)
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> = 'U': B is upper bidiagonal; !> = 'L': B is lower bidiagonal. !>
[in]	JOBZ	!> JOBZ is CHARACTER*1 !> = 'N': Compute singular values only; !> = 'V': Compute singular values and singular vectors. !>
[in]	RANGE	!> 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. !>
[in]	N	!> N is INTEGER !> The order of the bidiagonal matrix. N >= 0. !>
[in]	D	!> D is DOUBLE PRECISION array, dimension (N) !> The n diagonal elements of the bidiagonal matrix B. !>
[in]	E	!> 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. !>
[in]	VL	!> 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'. !>
[in]	VU	!> 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'. !>
[in]	IL	!> 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'. !>
[in]	IU	!> 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'. !>
[out]	NS	!> 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. !>
[out]	S	!> S is DOUBLE PRECISION array, dimension (N) !> The first NS elements contain the selected singular values in !> ascending order. !>
[out]	Z	!> 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. !>
[in]	LDZ	!> LDZ is INTEGER !> The leading dimension of the array Z. LDZ >= 1, and if !> JOBZ = 'V', LDZ >= max(2,N*2). !>
[out]	WORK	!> WORK is DOUBLE PRECISION array, dimension (14*N) !>
[out]	IWORK	!> 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. !>
[out]	INFO	!> 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. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 222 of file dbdsvdx.f.

*
*  -- 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
*

Here is the call graph for this function:

Here is the caller graph for this function: