dlasd7.f

*> \brief \b DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLASD7 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd7.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd7.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd7.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
*                          VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
*                          PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
*                          C, S, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
*      $                   NR, SQRE
*       DOUBLE PRECISION   ALPHA, BETA, C, S
*       ..
*       .. Array Arguments ..
*       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
*      $                   IDXQ( * ), PERM( * )
*       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
*      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
*      $                   ZW( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLASD7 merges the two sets of singular values together into a single
*> sorted set. Then it tries to deflate the size of the problem. There
*> are two ways in which deflation can occur:  when two or more singular
*> values are close together or if there is a tiny entry in the Z
*> vector. For each such occurrence the order of the related
*> secular equation problem is reduced by one.
*>
*> DLASD7 is called from DLASD6.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] ICOMPQ
*> \verbatim
*>          ICOMPQ is INTEGER
*>          Specifies whether singular vectors are to be computed
*>          in compact form, as follows:
*>          = 0: Compute singular values only.
*>          = 1: Compute singular vectors of upper
*>               bidiagonal matrix in compact form.
*> \endverbatim
*>
*> \param[in] NL
*> \verbatim
*>          NL is INTEGER
*>         The row dimension of the upper block. NL >= 1.
*> \endverbatim
*>
*> \param[in] NR
*> \verbatim
*>          NR is INTEGER
*>         The row dimension of the lower block. NR >= 1.
*> \endverbatim
*>
*> \param[in] SQRE
*> \verbatim
*>          SQRE is INTEGER
*>         = 0: the lower block is an NR-by-NR square matrix.
*>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
*>
*>         The bidiagonal matrix has
*>         N = NL + NR + 1 rows and
*>         M = N + SQRE >= N columns.
*> \endverbatim
*>
*> \param[out] K
*> \verbatim
*>          K is INTEGER
*>         Contains the dimension of the non-deflated matrix, this is
*>         the order of the related secular equation. 1 <= K <=N.
*> \endverbatim
*>
*> \param[in,out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension ( N )
*>         On entry D contains the singular values of the two submatrices
*>         to be combined. On exit D contains the trailing (N-K) updated
*>         singular values (those which were deflated) sorted into
*>         increasing order.
*> \endverbatim
*>
*> \param[out] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension ( M )
*>         On exit Z contains the updating row vector in the secular
*>         equation.
*> \endverbatim
*>
*> \param[out] ZW
*> \verbatim
*>          ZW is DOUBLE PRECISION array, dimension ( M )
*>         Workspace for Z.
*> \endverbatim
*>
*> \param[in,out] VF
*> \verbatim
*>          VF is DOUBLE PRECISION array, dimension ( M )
*>         On entry, VF(1:NL+1) contains the first components of all
*>         right singular vectors of the upper block; and VF(NL+2:M)
*>         contains the first components of all right singular vectors
*>         of the lower block. On exit, VF contains the first components
*>         of all right singular vectors of the bidiagonal matrix.
*> \endverbatim
*>
*> \param[out] VFW
*> \verbatim
*>          VFW is DOUBLE PRECISION array, dimension ( M )
*>         Workspace for VF.
*> \endverbatim
*>
*> \param[in,out] VL
*> \verbatim
*>          VL is DOUBLE PRECISION array, dimension ( M )
*>         On entry, VL(1:NL+1) contains the  last components of all
*>         right singular vectors of the upper block; and VL(NL+2:M)
*>         contains the last components of all right singular vectors
*>         of the lower block. On exit, VL contains the last components
*>         of all right singular vectors of the bidiagonal matrix.
*> \endverbatim
*>
*> \param[out] VLW
*> \verbatim
*>          VLW is DOUBLE PRECISION array, dimension ( M )
*>         Workspace for VL.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*>          ALPHA is DOUBLE PRECISION
*>         Contains the diagonal element associated with the added row.
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*>          BETA is DOUBLE PRECISION
*>         Contains the off-diagonal element associated with the added
*>         row.
*> \endverbatim
*>
*> \param[out] DSIGMA
*> \verbatim
*>          DSIGMA is DOUBLE PRECISION array, dimension ( N )
*>         Contains a copy of the diagonal elements (K-1 singular values
*>         and one zero) in the secular equation.
*> \endverbatim
*>
*> \param[out] IDX
*> \verbatim
*>          IDX is INTEGER array, dimension ( N )
*>         This will contain the permutation used to sort the contents of
*>         D into ascending order.
*> \endverbatim
*>
*> \param[out] IDXP
*> \verbatim
*>          IDXP is INTEGER array, dimension ( N )
*>         This will contain the permutation used to place deflated
*>         values of D at the end of the array. On output IDXP(2:K)
*>         points to the nondeflated D-values and IDXP(K+1:N)
*>         points to the deflated singular values.
*> \endverbatim
*>
*> \param[in] IDXQ
*> \verbatim
*>          IDXQ is INTEGER array, dimension ( N )
*>         This contains the permutation which separately sorts the two
*>         sub-problems in D into ascending order.  Note that entries in
*>         the first half of this permutation must first be moved one
*>         position backward; and entries in the second half
*>         must first have NL+1 added to their values.
*> \endverbatim
*>
*> \param[out] PERM
*> \verbatim
*>          PERM is INTEGER array, dimension ( N )
*>         The permutations (from deflation and sorting) to be applied
*>         to each singular block. Not referenced if ICOMPQ = 0.
*> \endverbatim
*>
*> \param[out] GIVPTR
*> \verbatim
*>          GIVPTR is INTEGER
*>         The number of Givens rotations which took place in this
*>         subproblem. Not referenced if ICOMPQ = 0.
*> \endverbatim
*>
*> \param[out] GIVCOL
*> \verbatim
*>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
*>         Each pair of numbers indicates a pair of columns to take place
*>         in a Givens rotation. Not referenced if ICOMPQ = 0.
*> \endverbatim
*>
*> \param[in] LDGCOL
*> \verbatim
*>          LDGCOL is INTEGER
*>         The leading dimension of GIVCOL, must be at least N.
*> \endverbatim
*>
*> \param[out] GIVNUM
*> \verbatim
*>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
*>         Each number indicates the C or S value to be used in the
*>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
*> \endverbatim
*>
*> \param[in] LDGNUM
*> \verbatim
*>          LDGNUM is INTEGER
*>         The leading dimension of GIVNUM, must be at least N.
*> \endverbatim
*>
*> \param[out] C
*> \verbatim
*>          C is DOUBLE PRECISION
*>         C contains garbage if SQRE =0 and the C-value of a Givens
*>         rotation related to the right null space if SQRE = 1.
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          S is DOUBLE PRECISION
*>         S contains garbage if SQRE =0 and the S-value of a Givens
*>         rotation related to the right null space if SQRE = 1.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>         = 0:  successful exit.
*>         < 0:  if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup lasd7
*
*> \par Contributors:
*  ==================
*>
*>     Ming Gu and Huan Ren, Computer Science Division, University of
*>     California at Berkeley, USA
*>
*  =====================================================================
      SUBROUTINE dlasd7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
     $                   VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
     $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
     $                   C, S, INFO )
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
     $                   nr, sqre
      DOUBLE PRECISION   ALPHA, BETA, C, S
*     ..
*     .. Array Arguments ..
      INTEGER            GIVCOL( ldgcol, * ), IDX( * ), IDXP( * ),
     $                   idxq( * ), perm( * )
      DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( ldgnum, * ),
     $                   vf( * ), vfw( * ), vl( * ), vlw( * ), z( * ),
     $                   zw( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
      parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
     $                   eight = 8.0d+0 )
*     ..
*     .. Local Scalars ..
*
      INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
     $                   nlp1, nlp2
      DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlamrg, drot, xerbla
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLAPY2
      EXTERNAL           dlamch, dlapy2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      n = nl + nr + 1
      m = n + sqre
*
      IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
         info = -1
      ELSE IF( nl.LT.1 ) THEN
         info = -2
      ELSE IF( nr.LT.1 ) THEN
         info = -3
      ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
         info = -4
      ELSE IF( ldgcol.LT.n ) THEN
         info = -22
      ELSE IF( ldgnum.LT.n ) THEN
         info = -24
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DLASD7', -info )
         RETURN
      END IF
*
      nlp1 = nl + 1
      nlp2 = nl + 2
      IF( icompq.EQ.1 ) THEN
         givptr = 0
      END IF
*
*     Generate the first part of the vector Z and move the singular
*     values in the first part of D one position backward.
*
      z1 = alpha*vl( nlp1 )
      vl( nlp1 ) = zero
      tau = vf( nlp1 )
      DO 10 i = nl, 1, -1
         z( i+1 ) = alpha*vl( i )
         vl( i ) = zero
         vf( i+1 ) = vf( i )
         d( i+1 ) = d( i )
         idxq( i+1 ) = idxq( i ) + 1
   10 CONTINUE
      vf( 1 ) = tau
*
*     Generate the second part of the vector Z.
*
      DO 20 i = nlp2, m
         z( i ) = beta*vf( i )
         vf( i ) = zero
   20 CONTINUE
*
*     Sort the singular values into increasing order
*
      DO 30 i = nlp2, n
         idxq( i ) = idxq( i ) + nlp1
   30 CONTINUE
*
*     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
*
      DO 40 i = 2, n
         dsigma( i ) = d( idxq( i ) )
         zw( i ) = z( idxq( i ) )
         vfw( i ) = vf( idxq( i ) )
         vlw( i ) = vl( idxq( i ) )
   40 CONTINUE
*
      CALL dlamrg( nl, nr, dsigma( 2 ), 1, 1, idx( 2 ) )
*
      DO 50 i = 2, n
         idxi = 1 + idx( i )
         d( i ) = dsigma( idxi )
         z( i ) = zw( idxi )
         vf( i ) = vfw( idxi )
         vl( i ) = vlw( idxi )
   50 CONTINUE
*
*     Calculate the allowable deflation tolerance
*
      eps = dlamch( 'Epsilon' )
      tol = max( abs( alpha ), abs( beta ) )
      tol = eight*eight*eps*max( abs( d( n ) ), tol )
*
*     There are 2 kinds of deflation -- first a value in the z-vector
*     is small, second two (or more) singular values are very close
*     together (their difference is small).
*
*     If the value in the z-vector is small, we simply permute the
*     array so that the corresponding singular value is moved to the
*     end.
*
*     If two values in the D-vector are close, we perform a two-sided
*     rotation designed to make one of the corresponding z-vector
*     entries zero, and then permute the array so that the deflated
*     singular value is moved to the end.
*
*     If there are multiple singular values then the problem deflates.
*     Here the number of equal singular values are found.  As each equal
*     singular value is found, an elementary reflector is computed to
*     rotate the corresponding singular subspace so that the
*     corresponding components of Z are zero in this new basis.
*
      k = 1
      k2 = n + 1
      DO 60 j = 2, n
         IF( abs( z( j ) ).LE.tol ) THEN
*
*           Deflate due to small z component.
*
            k2 = k2 - 1
            idxp( k2 ) = j
            IF( j.EQ.n )
     $         GO TO 100
         ELSE
            jprev = j
            GO TO 70
         END IF
   60 CONTINUE
   70 CONTINUE
      j = jprev
   80 CONTINUE
      j = j + 1
      IF( j.GT.n )
     $   GO TO 90
      IF( abs( z( j ) ).LE.tol ) THEN
*
*        Deflate due to small z component.
*
         k2 = k2 - 1
         idxp( k2 ) = j
      ELSE
*
*        Check if singular values are close enough to allow deflation.
*
         IF( abs( d( j )-d( jprev ) ).LE.tol ) THEN
*
*           Deflation is possible.
*
            s = z( jprev )
            c = z( j )
*
*           Find sqrt(a**2+b**2) without overflow or
*           destructive underflow.
*
            tau = dlapy2( c, s )
            z( j ) = tau
            z( jprev ) = zero
            c = c / tau
            s = -s / tau
*
*           Record the appropriate Givens rotation
*
            IF( icompq.EQ.1 ) THEN
               givptr = givptr + 1
               idxjp = idxq( idx( jprev )+1 )
               idxj = idxq( idx( j )+1 )
               IF( idxjp.LE.nlp1 ) THEN
                  idxjp = idxjp - 1
               END IF
               IF( idxj.LE.nlp1 ) THEN
                  idxj = idxj - 1
               END IF
               givcol( givptr, 2 ) = idxjp
               givcol( givptr, 1 ) = idxj
               givnum( givptr, 2 ) = c
               givnum( givptr, 1 ) = s
            END IF
            CALL drot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
            CALL drot( 1, vl( jprev ), 1, vl( j ), 1, c, s )
            k2 = k2 - 1
            idxp( k2 ) = jprev
            jprev = j
         ELSE
            k = k + 1
            zw( k ) = z( jprev )
            dsigma( k ) = d( jprev )
            idxp( k ) = jprev
            jprev = j
         END IF
      END IF
      GO TO 80
   90 CONTINUE
*
*     Record the last singular value.
*
      k = k + 1
      zw( k ) = z( jprev )
      dsigma( k ) = d( jprev )
      idxp( k ) = jprev
*
  100 CONTINUE
*
*     Sort the singular values into DSIGMA. The singular values which
*     were not deflated go into the first K slots of DSIGMA, except
*     that DSIGMA(1) is treated separately.
*
      DO 110 j = 2, n
         jp = idxp( j )
         dsigma( j ) = d( jp )
         vfw( j ) = vf( jp )
         vlw( j ) = vl( jp )
  110 CONTINUE
      IF( icompq.EQ.1 ) THEN
         DO 120 j = 2, n
            jp = idxp( j )
            perm( j ) = idxq( idx( jp )+1 )
            IF( perm( j ).LE.nlp1 ) THEN
               perm( j ) = perm( j ) - 1
            END IF
  120    CONTINUE
      END IF
*
*     The deflated singular values go back into the last N - K slots of
*     D.
*
      CALL dcopy( n-k, dsigma( k+1 ), 1, d( k+1 ), 1 )
*
*     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
*     VL(M).
*
      dsigma( 1 ) = zero
      hlftol = tol / two
      IF( abs( dsigma( 2 ) ).LE.hlftol )
     $   dsigma( 2 ) = hlftol
      IF( m.GT.n ) THEN
         z( 1 ) = dlapy2( z1, z( m ) )
         IF( z( 1 ).LE.tol ) THEN
            c = one
            s = zero
            z( 1 ) = tol
         ELSE
            c = z1 / z( 1 )
            s = -z( m ) / z( 1 )
         END IF
         CALL drot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
         CALL drot( 1, vl( m ), 1, vl( 1 ), 1, c, s )
      ELSE
         IF( abs( z1 ).LE.tol ) THEN
            z( 1 ) = tol
         ELSE
            z( 1 ) = z1
         END IF
      END IF
*
*     Restore Z, VF, and VL.
*
      CALL dcopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
      CALL dcopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
      CALL dcopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
*
      RETURN
*
*     End of DLASD7
*
      END