dlasq3.f

*> \brief \b DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLASQ3 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
*                          ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
*                          DN2, G, TAU )
*
*       .. Scalar Arguments ..
*       LOGICAL            IEEE
*       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
*       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
*      $                   QMAX, SIGMA, TAU
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   Z( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
*> In case of failure it changes shifts, and tries again until output
*> is positive.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] I0
*> \verbatim
*>          I0 is INTEGER
*>         First index.
*> \endverbatim
*>
*> \param[in,out] N0
*> \verbatim
*>          N0 is INTEGER
*>         Last index.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension ( 4*N0 )
*>         Z holds the qd array.
*> \endverbatim
*>
*> \param[in,out] PP
*> \verbatim
*>          PP is INTEGER
*>         PP=0 for ping, PP=1 for pong.
*>         PP=2 indicates that flipping was applied to the Z array
*>         and that the initial tests for deflation should not be
*>         performed.
*> \endverbatim
*>
*> \param[out] DMIN
*> \verbatim
*>          DMIN is DOUBLE PRECISION
*>         Minimum value of d.
*> \endverbatim
*>
*> \param[out] SIGMA
*> \verbatim
*>          SIGMA is DOUBLE PRECISION
*>         Sum of shifts used in current segment.
*> \endverbatim
*>
*> \param[in,out] DESIG
*> \verbatim
*>          DESIG is DOUBLE PRECISION
*>         Lower order part of SIGMA
*> \endverbatim
*>
*> \param[in] QMAX
*> \verbatim
*>          QMAX is DOUBLE PRECISION
*>         Maximum value of q.
*> \endverbatim
*>
*> \param[in,out] NFAIL
*> \verbatim
*>          NFAIL is INTEGER
*>         Increment NFAIL by 1 each time the shift was too big.
*> \endverbatim
*>
*> \param[in,out] ITER
*> \verbatim
*>          ITER is INTEGER
*>         Increment ITER by 1 for each iteration.
*> \endverbatim
*>
*> \param[in,out] NDIV
*> \verbatim
*>          NDIV is INTEGER
*>         Increment NDIV by 1 for each division.
*> \endverbatim
*>
*> \param[in] IEEE
*> \verbatim
*>          IEEE is LOGICAL
*>         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
*> \endverbatim
*>
*> \param[in,out] TTYPE
*> \verbatim
*>          TTYPE is INTEGER
*>         Shift type.
*> \endverbatim
*>
*> \param[in,out] DMIN1
*> \verbatim
*>          DMIN1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DMIN2
*> \verbatim
*>          DMIN2 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DN
*> \verbatim
*>          DN is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DN1
*> \verbatim
*>          DN1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DN2
*> \verbatim
*>          DN2 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] G
*> \verbatim
*>          G is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] TAU
*> \verbatim
*>          TAU is DOUBLE PRECISION
*>
*>         These are passed as arguments in order to save their values
*>         between calls to DLASQ3.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup lasq3
*
*  =====================================================================
      SUBROUTINE dlasq3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
     $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
     $                   DN2, G, TAU )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            IEEE
      INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
      DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
     $                   qmax, sigma, tau
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   Z( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   CBIAS
      parameter( cbias = 1.50d0 )
      DOUBLE PRECISION   ZERO, QURTR, HALF, ONE, TWO, HUNDRD
      parameter( zero = 0.0d0, qurtr = 0.250d0, half = 0.5d0,
     $                     one = 1.0d0, two = 2.0d0, hundrd = 100.0d0 )
*     ..
*     .. Local Scalars ..
      INTEGER            IPN4, J4, N0IN, NN, TTYPE
      DOUBLE PRECISION   EPS, S, T, TEMP, TOL, TOL2
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlasq4, dlasq5, dlasq6
*     ..
*     .. External Function ..
      DOUBLE PRECISION   DLAMCH
      LOGICAL            DISNAN
      EXTERNAL           disnan, dlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, sqrt
*     ..
*     .. Executable Statements ..
*
      n0in = n0
      eps = dlamch( 'Precision' )
      tol = eps*hundrd
      tol2 = tol**2
*
*     Check for deflation.
*
   10 CONTINUE
*
      IF( n0.LT.i0 )
     $   RETURN
      IF( n0.EQ.i0 )
     $   GO TO 20
      nn = 4*n0 + pp
      IF( n0.EQ.( i0+1 ) )
     $   GO TO 40
*
*     Check whether E(N0-1) is negligible, 1 eigenvalue.
*
      IF( z( nn-5 ).GT.tol2*( sigma+z( nn-3 ) ) .AND.
     $    z( nn-2*pp-4 ).GT.tol2*z( nn-7 ) )
     $   GO TO 30
*
   20 CONTINUE
*
      z( 4*n0-3 ) = z( 4*n0+pp-3 ) + sigma
      n0 = n0 - 1
      GO TO 10
*
*     Check  whether E(N0-2) is negligible, 2 eigenvalues.
*
   30 CONTINUE
*
      IF( z( nn-9 ).GT.tol2*sigma .AND.
     $    z( nn-2*pp-8 ).GT.tol2*z( nn-11 ) )
     $   GO TO 50
*
   40 CONTINUE
*
      IF( z( nn-3 ).GT.z( nn-7 ) ) THEN
         s = z( nn-3 )
         z( nn-3 ) = z( nn-7 )
         z( nn-7 ) = s
      END IF
      t = half*( ( z( nn-7 )-z( nn-3 ) )+z( nn-5 ) )
      IF( z( nn-5 ).GT.z( nn-3 )*tol2.AND.t.NE.zero ) THEN
         s = z( nn-3 )*( z( nn-5 ) / t )
         IF( s.LE.t ) THEN
            s = z( nn-3 )*( z( nn-5 ) /
     $          ( t*( one+sqrt( one+s / t ) ) ) )
         ELSE
            s = z( nn-3 )*( z( nn-5 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
         END IF
         t = z( nn-7 ) + ( s+z( nn-5 ) )
         z( nn-3 ) = z( nn-3 )*( z( nn-7 ) / t )
         z( nn-7 ) = t
      END IF
      z( 4*n0-7 ) = z( nn-7 ) + sigma
      z( 4*n0-3 ) = z( nn-3 ) + sigma
      n0 = n0 - 2
      GO TO 10
*
   50 CONTINUE
      IF( pp.EQ.2 )
     $   pp = 0
*
*     Reverse the qd-array, if warranted.
*
      IF( dmin.LE.zero .OR. n0.LT.n0in ) THEN
         IF( cbias*z( 4*i0+pp-3 ).LT.z( 4*n0+pp-3 ) ) THEN
            ipn4 = 4*( i0+n0 )
            DO 60 j4 = 4*i0, 2*( i0+n0-1 ), 4
               temp = z( j4-3 )
               z( j4-3 ) = z( ipn4-j4-3 )
               z( ipn4-j4-3 ) = temp
               temp = z( j4-2 )
               z( j4-2 ) = z( ipn4-j4-2 )
               z( ipn4-j4-2 ) = temp
               temp = z( j4-1 )
               z( j4-1 ) = z( ipn4-j4-5 )
               z( ipn4-j4-5 ) = temp
               temp = z( j4 )
               z( j4 ) = z( ipn4-j4-4 )
               z( ipn4-j4-4 ) = temp
   60       CONTINUE
            IF( n0-i0.LE.4 ) THEN
               z( 4*n0+pp-1 ) = z( 4*i0+pp-1 )
               z( 4*n0-pp ) = z( 4*i0-pp )
            END IF
            dmin2 = min( dmin2, z( 4*n0+pp-1 ) )
            z( 4*n0+pp-1 ) = min( z( 4*n0+pp-1 ), z( 4*i0+pp-1 ),
     $                            z( 4*i0+pp+3 ) )
            z( 4*n0-pp ) = min( z( 4*n0-pp ), z( 4*i0-pp ),
     $                          z( 4*i0-pp+4 ) )
            qmax = max( qmax, z( 4*i0+pp-3 ), z( 4*i0+pp+1 ) )
            dmin = -zero
         END IF
      END IF
*
*     Choose a shift.
*
      CALL dlasq4( i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1,
     $             dn2, tau, ttype, g )
*
*     Call dqds until DMIN > 0.
*
   70 CONTINUE
*
      CALL dlasq5( i0, n0, z, pp, tau, sigma, dmin, dmin1, dmin2, dn,
     $             dn1, dn2, ieee, eps )
*
      ndiv = ndiv + ( n0-i0+2 )
      iter = iter + 1
*
*     Check status.
*
      IF( dmin.GE.zero .AND. dmin1.GE.zero ) THEN
*
*        Success.
*
         GO TO 90
*
      ELSE IF( dmin.LT.zero .AND. dmin1.GT.zero .AND.
     $         z( 4*( n0-1 )-pp ).LT.tol*( sigma+dn1 ) .AND.
     $         abs( dn ).LT.tol*sigma ) THEN
*
*        Convergence hidden by negative DN.
*
         z( 4*( n0-1 )-pp+2 ) = zero
         dmin = zero
         GO TO 90
      ELSE IF( dmin.LT.zero ) THEN
*
*        TAU too big. Select new TAU and try again.
*
         nfail = nfail + 1
         IF( ttype.LT.-22 ) THEN
*
*           Failed twice. Play it safe.
*
            tau = zero
         ELSE IF( dmin1.GT.zero ) THEN
*
*           Late failure. Gives excellent shift.
*
            tau = ( tau+dmin )*( one-two*eps )
            ttype = ttype - 11
         ELSE
*
*           Early failure. Divide by 4.
*
            tau = qurtr*tau
            ttype = ttype - 12
         END IF
         GO TO 70
      ELSE IF( disnan( dmin ) ) THEN
*
*        NaN.
*
         IF( tau.EQ.zero ) THEN
            GO TO 80
         ELSE
            tau = zero
            GO TO 70
         END IF
      ELSE
*
*        Possible underflow. Play it safe.
*
         GO TO 80
      END IF
*
*     Risk of underflow.
*
   80 CONTINUE
      CALL dlasq6( i0, n0, z, pp, dmin, dmin1, dmin2, dn, dn1, dn2 )
      ndiv = ndiv + ( n0-i0+2 )
      iter = iter + 1
      tau = zero
*
   90 CONTINUE
      IF( tau.LT.sigma ) THEN
         desig = desig + tau
         t = sigma + desig
         desig = desig - ( t-sigma )
      ELSE
         t = sigma + tau
         desig = sigma - ( t-tau ) + desig
      END IF
      sigma = t
*
      RETURN
*
*     End of DLASQ3
*
      END