slasq2.f

Go to the documentation of this file.

*> \brief \b SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download SLASQ2 + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq2.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq2.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq2.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE SLASQ2( N, Z, INFO )
* 
*       .. Scalar Arguments ..
*       INTEGER            INFO, N
*       ..
*       .. Array Arguments ..
*       REAL               Z( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SLASQ2 computes all the eigenvalues of the symmetric positive 
*> definite tridiagonal matrix associated with the qd array Z to high
*> relative accuracy are computed to high relative accuracy, in the
*> absence of denormalization, underflow and overflow.
*>
*> To see the relation of Z to the tridiagonal matrix, let L be a
*> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
*> let U be an upper bidiagonal matrix with 1's above and diagonal
*> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
*> symmetric tridiagonal to which it is similar.
*>
*> Note : SLASQ2 defines a logical variable, IEEE, which is true
*> on machines which follow ieee-754 floating-point standard in their
*> handling of infinities and NaNs, and false otherwise. This variable
*> is passed to SLASQ3.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>        The number of rows and columns in the matrix. N >= 0.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is REAL array, dimension ( 4*N )
*>        On entry Z holds the qd array. On exit, entries 1 to N hold
*>        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
*>        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
*>        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
*>        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
*>        shifts that failed.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>        = 0: successful exit
*>        < 0: if the i-th argument is a scalar and had an illegal
*>             value, then INFO = -i, if the i-th argument is an
*>             array and the j-entry had an illegal value, then
*>             INFO = -(i*100+j)
*>        > 0: the algorithm failed
*>              = 1, a split was marked by a positive value in E
*>              = 2, current block of Z not diagonalized after 100*N
*>                   iterations (in inner while loop).  On exit Z holds
*>                   a qd array with the same eigenvalues as the given Z.
*>              = 3, termination criterion of outer while loop not met 
*>                   (program created more than N unreduced blocks)
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date September 2012
*
*> \ingroup auxOTHERcomputational
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  Local Variables: I0:N0 defines a current unreduced segment of Z.
*>  The shifts are accumulated in SIGMA. Iteration count is in ITER.
*>  Ping-pong is controlled by PP (alternates between 0 and 1).
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE slasq2( N, Z, INFO )
*
*  -- LAPACK computational routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      INTEGER            info, n
*     ..
*     .. Array Arguments ..
      REAL               z( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               cbias
      parameter( cbias = 1.50e0 )
      REAL               zero, half, one, two, four, hundrd
      parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0,
     $                     two = 2.0e0, four = 4.0e0, hundrd = 100.0e0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ieee
      INTEGER            i0, i4, iinfo, ipn4, iter, iwhila, iwhilb, k,
     $                   kmin, n0, nbig, ndiv, nfail, pp, splt, ttype,
     $                   i1, n1
      REAL               d, dee, deemin, desig, dmin, dmin1, dmin2, dn,
     $                   dn1, dn2, e, emax, emin, eps, g, oldemn, qmax,
     $                   qmin, s, safmin, sigma, t, tau, temp, tol,
     $                   tol2, trace, zmax, tempe, tempq
*     ..
*     .. External Subroutines ..
      EXTERNAL           slasq3, slasrt, xerbla
*     ..
*     .. External Functions ..
      INTEGER            ilaenv
      REAL               slamch
      EXTERNAL           ilaenv, slamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, REAL, sqrt
*     ..
*     .. Executable Statements ..
*      
*     Test the input arguments.
*     (in case SLASQ2 is not called by SLASQ1)
*
      info = 0
      eps = slamch( 'Precision' )
      safmin = slamch( 'Safe minimum' )
      tol = eps*hundrd
      tol2 = tol**2
*
      IF( n.LT.0 ) THEN
         info = -1
         CALL xerbla( 'SLASQ2', 1 )
         return
      ELSE IF( n.EQ.0 ) THEN
         return
      ELSE IF( n.EQ.1 ) THEN
*
*        1-by-1 case.
*
         IF( z( 1 ).LT.zero ) THEN
            info = -201
            CALL xerbla( 'SLASQ2', 2 )
         END IF
         return
      ELSE IF( n.EQ.2 ) THEN
*
*        2-by-2 case.
*
         IF( z( 2 ).LT.zero .OR. z( 3 ).LT.zero ) THEN
            info = -2
            CALL xerbla( 'SLASQ2', 2 )
            return
         ELSE IF( z( 3 ).GT.z( 1 ) ) THEN
            d = z( 3 )
            z( 3 ) = z( 1 )
            z( 1 ) = d
         END IF
         z( 5 ) = z( 1 ) + z( 2 ) + z( 3 )
         IF( z( 2 ).GT.z( 3 )*tol2 ) THEN
            t = half*( ( z( 1 )-z( 3 ) )+z( 2 ) ) 
            s = z( 3 )*( z( 2 ) / t )
            IF( s.LE.t ) THEN
               s = z( 3 )*( z( 2 ) / ( t*( one+sqrt( one+s / t ) ) ) )
            ELSE
               s = z( 3 )*( z( 2 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
            END IF
            t = z( 1 ) + ( s+z( 2 ) )
            z( 3 ) = z( 3 )*( z( 1 ) / t )
            z( 1 ) = t
         END IF
         z( 2 ) = z( 3 )
         z( 6 ) = z( 2 ) + z( 1 )
         return
      END IF
*
*     Check for negative data and compute sums of q's and e's.
*
      z( 2*n ) = zero
      emin = z( 2 )
      qmax = zero
      zmax = zero
      d = zero
      e = zero
*
      DO 10 k = 1, 2*( n-1 ), 2
         IF( z( k ).LT.zero ) THEN
            info = -( 200+k )
            CALL xerbla( 'SLASQ2', 2 )
            return
         ELSE IF( z( k+1 ).LT.zero ) THEN
            info = -( 200+k+1 )
            CALL xerbla( 'SLASQ2', 2 )
            return
         END IF
         d = d + z( k )
         e = e + z( k+1 )
         qmax = max( qmax, z( k ) )
         emin = min( emin, z( k+1 ) )
         zmax = max( qmax, zmax, z( k+1 ) )
   10 continue
      IF( z( 2*n-1 ).LT.zero ) THEN
         info = -( 200+2*n-1 )
         CALL xerbla( 'SLASQ2', 2 )
         return
      END IF
      d = d + z( 2*n-1 )
      qmax = max( qmax, z( 2*n-1 ) )
      zmax = max( qmax, zmax )
*
*     Check for diagonality.
*
      IF( e.EQ.zero ) THEN
         DO 20 k = 2, n
            z( k ) = z( 2*k-1 )
   20    continue
         CALL slasrt( 'D', n, z, iinfo )
         z( 2*n-1 ) = d
         return
      END IF
*
      trace = d + e
*
*     Check for zero data.
*
      IF( trace.EQ.zero ) THEN
         z( 2*n-1 ) = zero
         return
      END IF
*         
*     Check whether the machine is IEEE conformable.
*         
*     IEEE = ILAENV( 10, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
*    $       ILAENV( 11, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1      
*
*     [11/15/2008] The case IEEE=.TRUE. has a problem in single precision with
*     some the test matrices of type 16. The double precision code is fine.
*
      ieee = .false.
*         
*     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
*
      DO 30 k = 2*n, 2, -2
         z( 2*k ) = zero 
         z( 2*k-1 ) = z( k ) 
         z( 2*k-2 ) = zero 
         z( 2*k-3 ) = z( k-1 ) 
   30 continue
*
      i0 = 1
      n0 = n
*
*     Reverse the qd-array, if warranted.
*
      IF( cbias*z( 4*i0-3 ).LT.z( 4*n0-3 ) ) THEN
         ipn4 = 4*( i0+n0 )
         DO 40 i4 = 4*i0, 2*( i0+n0-1 ), 4
            temp = z( i4-3 )
            z( i4-3 ) = z( ipn4-i4-3 )
            z( ipn4-i4-3 ) = temp
            temp = z( i4-1 )
            z( i4-1 ) = z( ipn4-i4-5 )
            z( ipn4-i4-5 ) = temp
   40    continue
      END IF
*
*     Initial split checking via dqd and Li's test.
*
      pp = 0
*
      DO 80 k = 1, 2
*
         d = z( 4*n0+pp-3 )
         DO 50 i4 = 4*( n0-1 ) + pp, 4*i0 + pp, -4
            IF( z( i4-1 ).LE.tol2*d ) THEN
               z( i4-1 ) = -zero
               d = z( i4-3 )
            ELSE
               d = z( i4-3 )*( d / ( d+z( i4-1 ) ) )
            END IF
   50    continue
*
*        dqd maps Z to ZZ plus Li's test.
*
         emin = z( 4*i0+pp+1 )
         d = z( 4*i0+pp-3 )
         DO 60 i4 = 4*i0 + pp, 4*( n0-1 ) + pp, 4
            z( i4-2*pp-2 ) = d + z( i4-1 )
            IF( z( i4-1 ).LE.tol2*d ) THEN
               z( i4-1 ) = -zero
               z( i4-2*pp-2 ) = d
               z( i4-2*pp ) = zero
               d = z( i4+1 )
            ELSE IF( safmin*z( i4+1 ).LT.z( i4-2*pp-2 ) .AND.
     $               safmin*z( i4-2*pp-2 ).LT.z( i4+1 ) ) THEN
               temp = z( i4+1 ) / z( i4-2*pp-2 )
               z( i4-2*pp ) = z( i4-1 )*temp
               d = d*temp
            ELSE
               z( i4-2*pp ) = z( i4+1 )*( z( i4-1 ) / z( i4-2*pp-2 ) )
               d = z( i4+1 )*( d / z( i4-2*pp-2 ) )
            END IF
            emin = min( emin, z( i4-2*pp ) )
   60    continue 
         z( 4*n0-pp-2 ) = d
*
*        Now find qmax.
*
         qmax = z( 4*i0-pp-2 )
         DO 70 i4 = 4*i0 - pp + 2, 4*n0 - pp - 2, 4
            qmax = max( qmax, z( i4 ) )
   70    continue
*
*        Prepare for the next iteration on K.
*
         pp = 1 - pp
   80 continue
*
*     Initialise variables to pass to SLASQ3.
*
      ttype = 0
      dmin1 = zero
      dmin2 = zero
      dn    = zero
      dn1   = zero
      dn2   = zero
      g     = zero
      tau   = zero
*
      iter = 2
      nfail = 0
      ndiv = 2*( n0-i0 )
*
      DO 160 iwhila = 1, n + 1
         IF( n0.LT.1 ) 
     $      go to 170
*
*        While array unfinished do 
*
*        E(N0) holds the value of SIGMA when submatrix in I0:N0
*        splits from the rest of the array, but is negated.
*      
         desig = zero
         IF( n0.EQ.n ) THEN
            sigma = zero
         ELSE
            sigma = -z( 4*n0-1 )
         END IF
         IF( sigma.LT.zero ) THEN
            info = 1
            return
         END IF
*
*        Find last unreduced submatrix's top index I0, find QMAX and
*        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
*
         emax = zero 
         IF( n0.GT.i0 ) THEN
            emin = abs( z( 4*n0-5 ) )
         ELSE
            emin = zero
         END IF
         qmin = z( 4*n0-3 )
         qmax = qmin
         DO 90 i4 = 4*n0, 8, -4
            IF( z( i4-5 ).LE.zero )
     $         go to 100
            IF( qmin.GE.four*emax ) THEN
               qmin = min( qmin, z( i4-3 ) )
               emax = max( emax, z( i4-5 ) )
            END IF
            qmax = max( qmax, z( i4-7 )+z( i4-5 ) )
            emin = min( emin, z( i4-5 ) )
   90    continue
         i4 = 4 
*
  100    continue
         i0 = i4 / 4
         pp = 0
*
         IF( n0-i0.GT.1 ) THEN
            dee = z( 4*i0-3 )
            deemin = dee
            kmin = i0
            DO 110 i4 = 4*i0+1, 4*n0-3, 4
               dee = z( i4 )*( dee /( dee+z( i4-2 ) ) )
               IF( dee.LE.deemin ) THEN
                  deemin = dee
                  kmin = ( i4+3 )/4
               END IF
  110       continue
            IF( (kmin-i0)*2.LT.n0-kmin .AND. 
     $         deemin.LE.half*z(4*n0-3) ) THEN
               ipn4 = 4*( i0+n0 )
               pp = 2
               DO 120 i4 = 4*i0, 2*( i0+n0-1 ), 4
                  temp = z( i4-3 )
                  z( i4-3 ) = z( ipn4-i4-3 )
                  z( ipn4-i4-3 ) = temp
                  temp = z( i4-2 )
                  z( i4-2 ) = z( ipn4-i4-2 )
                  z( ipn4-i4-2 ) = temp
                  temp = z( i4-1 )
                  z( i4-1 ) = z( ipn4-i4-5 )
                  z( ipn4-i4-5 ) = temp
                  temp = z( i4 )
                  z( i4 ) = z( ipn4-i4-4 )
                  z( ipn4-i4-4 ) = temp
  120          continue
            END IF
         END IF
*
*        Put -(initial shift) into DMIN.
*
         dmin = -max( zero, qmin-two*sqrt( qmin )*sqrt( emax ) )
*
*        Now I0:N0 is unreduced. 
*        PP = 0 for ping, PP = 1 for pong.
*        PP = 2 indicates that flipping was applied to the Z array and
*               and that the tests for deflation upon entry in SLASQ3 
*               should not be performed.
*
         nbig = 100*( n0-i0+1 )
         DO 140 iwhilb = 1, nbig
            IF( i0.GT.n0 ) 
     $         go to 150
*
*           While submatrix unfinished take a good dqds step.
*
            CALL slasq3( i0, n0, z, pp, dmin, sigma, desig, qmax, nfail,
     $                   iter, ndiv, ieee, ttype, dmin1, dmin2, dn, dn1,
     $                   dn2, g, tau )
*
            pp = 1 - pp
*
*           When EMIN is very small check for splits.
*
            IF( pp.EQ.0 .AND. n0-i0.GE.3 ) THEN
               IF( z( 4*n0 ).LE.tol2*qmax .OR.
     $             z( 4*n0-1 ).LE.tol2*sigma ) THEN
                  splt = i0 - 1
                  qmax = z( 4*i0-3 )
                  emin = z( 4*i0-1 )
                  oldemn = z( 4*i0 )
                  DO 130 i4 = 4*i0, 4*( n0-3 ), 4
                     IF( z( i4 ).LE.tol2*z( i4-3 ) .OR.
     $                   z( i4-1 ).LE.tol2*sigma ) THEN
                        z( i4-1 ) = -sigma
                        splt = i4 / 4
                        qmax = zero
                        emin = z( i4+3 )
                        oldemn = z( i4+4 )
                     ELSE
                        qmax = max( qmax, z( i4+1 ) )
                        emin = min( emin, z( i4-1 ) )
                        oldemn = min( oldemn, z( i4 ) )
                     END IF
  130             continue
                  z( 4*n0-1 ) = emin
                  z( 4*n0 ) = oldemn
                  i0 = splt + 1
               END IF
            END IF
*
  140    continue
*
         info = 2
*       
*        Maximum number of iterations exceeded, restore the shift 
*        SIGMA and place the new d's and e's in a qd array.
*        This might need to be done for several blocks
*
         i1 = i0
         n1 = n0
 145     continue
         tempq = z( 4*i0-3 )
         z( 4*i0-3 ) = z( 4*i0-3 ) + sigma
         DO k = i0+1, n0
            tempe = z( 4*k-5 )
            z( 4*k-5 ) = z( 4*k-5 ) * (tempq / z( 4*k-7 ))
            tempq = z( 4*k-3 )
            z( 4*k-3 ) = z( 4*k-3 ) + sigma + tempe - z( 4*k-5 )
         END DO
*
*        Prepare to do this on the previous block if there is one
*
         IF( i1.GT.1 ) THEN
            n1 = i1-1
            DO WHILE( ( i1.GE.2 ) .AND. ( z(4*i1-5).GE.zero ) )
               i1 = i1 - 1
            END DO
            IF( i1.GE.1 ) THEN
               sigma = -z(4*n1-1)
               go to 145
            END IF
         END IF

         DO k = 1, n
            z( 2*k-1 ) = z( 4*k-3 )
*
*        Only the block 1..N0 is unfinished.  The rest of the e's
*        must be essentially zero, although sometimes other data
*        has been stored in them.
*
            IF( k.LT.n0 ) THEN
               z( 2*k ) = z( 4*k-1 )
            ELSE
               z( 2*k ) = 0
            END IF
         END DO
         return
*
*        end IWHILB
*
  150    continue
*
  160 continue
*
      info = 3
      return
*
*     end IWHILA   
*
  170 continue
*      
*     Move q's to the front.
*      
      DO 180 k = 2, n
         z( k ) = z( 4*k-3 )
  180 continue
*      
*     Sort and compute sum of eigenvalues.
*
      CALL slasrt( 'D', n, z, iinfo )
*
      e = zero
      DO 190 k = n, 1, -1
         e = e + z( k )
  190 continue
*
*     Store trace, sum(eigenvalues) and information on performance.
*
      z( 2*n+1 ) = trace 
      z( 2*n+2 ) = e
      z( 2*n+3 ) = REAL( iter )
      z( 2*n+4 ) = REAL( NDIV ) / REAL( n**2 )
      z( 2*n+5 ) = hundrd*nfail / REAL( iter )
      return
*
*     End of SLASQ2
*
      END