subroutine csteqr	(	character	COMPZ,
		integer	N,
		real, dimension( * )	D,
		real, dimension( * )	E,
		complex, dimension( ldz, * )	Z,
		integer	LDZ,
		real, dimension( * )	WORK,
		integer	INFO
	)

CSTEQR

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

Purpose:

 CSTEQR computes all eigenvalues and, optionally, eigenvectors of a
 symmetric tridiagonal matrix using the implicit QL or QR method.
 The eigenvectors of a full or band complex Hermitian matrix can also
 be found if CHETRD or CHPTRD or CHBTRD has been used to reduce this
 matrix to tridiagonal form.

Parameters

[in]	COMPZ	COMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'V': Compute eigenvalues and eigenvectors of the original Hermitian matrix. On entry, Z must contain the unitary matrix used to reduce the original matrix to tridiagonal form. = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix.
[in]	N	N is INTEGER The order of the matrix. N >= 0.
[in,out]	D	D is REAL array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.
[in,out]	E	E is REAL array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.
[in,out]	Z	Z is COMPLEX array, dimension (LDZ, N) On entry, if COMPZ = 'V', then Z contains the unitary matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original Hermitian matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.
[in]	LDZ	LDZ is INTEGER The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N).
[out]	WORK	WORK is REAL array, dimension (max(1,2*N-2)) If COMPZ = 'N', then WORK is not referenced.
[out]	INFO	INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is unitarily similar to the original matrix.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Date

November 2011

Definition at line 134 of file csteqr.f.

*
*  -- LAPACK computational routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      CHARACTER          compz
      INTEGER            info, ldz, n
*     ..
*     .. Array Arguments ..
      REAL               d( * ), e( * ), work( * )
      COMPLEX            z( ldz, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               zero, one, two, three
      parameter                ( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
     $                   three = 3.0e0 )
      COMPLEX            czero, cone
      parameter                ( czero = ( 0.0e0, 0.0e0 ),
     $                   cone = ( 1.0e0, 0.0e0 ) )
      INTEGER            maxit
      parameter                ( maxit = 30 )
*     ..
*     .. Local Scalars ..
      INTEGER            i, icompz, ii, iscale, j, jtot, k, l, l1, lend,
     $                   lendm1, lendp1, lendsv, lm1, lsv, m, mm, mm1,
     $                   nm1, nmaxit
      REAL               anorm, b, c, eps, eps2, f, g, p, r, rt1, rt2,
     $                   s, safmax, safmin, ssfmax, ssfmin, tst
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      REAL               slamch, slanst, slapy2
      EXTERNAL           lsame, slamch, slanst, slapy2
*     ..
*     .. External Subroutines ..
      EXTERNAL           claset, clasr, cswap, slae2, slaev2, slartg,
     $                   slascl, slasrt, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, sign, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
*
      IF( lsame( compz, 'N' ) ) THEN
         icompz = 0
      ELSE IF( lsame( compz, 'V' ) ) THEN
         icompz = 1
      ELSE IF( lsame( compz, 'I' ) ) THEN
         icompz = 2
      ELSE
         icompz = -1
      END IF
      IF( icompz.LT.0 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
     $         n ) ) ) THEN
         info = -6
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CSTEQR', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
*
      IF( n.EQ.1 ) THEN
         IF( icompz.EQ.2 )
     $      z( 1, 1 ) = cone
         RETURN
      END IF
*
*     Determine the unit roundoff and over/underflow thresholds.
*
      eps = slamch( 'E' )
      eps2 = eps**2
      safmin = slamch( 'S' )
      safmax = one / safmin
      ssfmax = sqrt( safmax ) / three
      ssfmin = sqrt( safmin ) / eps2
*
*     Compute the eigenvalues and eigenvectors of the tridiagonal
*     matrix.
*
      IF( icompz.EQ.2 )
     $   CALL claset( 'Full', n, n, czero, cone, z, ldz )
*
      nmaxit = n*maxit
      jtot = 0
*
*     Determine where the matrix splits and choose QL or QR iteration
*     for each block, according to whether top or bottom diagonal
*     element is smaller.
*
      l1 = 1
      nm1 = n - 1
*
   10 CONTINUE
      IF( l1.GT.n )
     $   GO TO 160
      IF( l1.GT.1 )
     $   e( l1-1 ) = zero
      IF( l1.LE.nm1 ) THEN
         DO 20 m = l1, nm1
            tst = abs( e( m ) )
            IF( tst.EQ.zero )
     $         GO TO 30
            IF( tst.LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
     $          1 ) ) ) )*eps ) THEN
               e( m ) = zero
               GO TO 30
            END IF
   20    CONTINUE
      END IF
      m = n
*
   30 CONTINUE
      l = l1
      lsv = l
      lend = m
      lendsv = lend
      l1 = m + 1
      IF( lend.EQ.l )
     $   GO TO 10
*
*     Scale submatrix in rows and columns L to LEND
*
      anorm = slanst( 'I', lend-l+1, d( l ), e( l ) )
      iscale = 0
      IF( anorm.EQ.zero )
     $   GO TO 10
      IF( anorm.GT.ssfmax ) THEN
         iscale = 1
         CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
     $                info )
         CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
     $                info )
      ELSE IF( anorm.LT.ssfmin ) THEN
         iscale = 2
         CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
     $                info )
         CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
     $                info )
      END IF
*
*     Choose between QL and QR iteration
*
      IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
         lend = lsv
         l = lendsv
      END IF
*
      IF( lend.GT.l ) THEN
*
*        QL Iteration
*
*        Look for small subdiagonal element.
*
   40    CONTINUE
         IF( l.NE.lend ) THEN
            lendm1 = lend - 1
            DO 50 m = l, lendm1
               tst = abs( e( m ) )**2
               IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
     $             safmin )GO TO 60
   50       CONTINUE
         END IF
*
         m = lend
*
   60    CONTINUE
         IF( m.LT.lend )
     $      e( m ) = zero
         p = d( l )
         IF( m.EQ.l )
     $      GO TO 80
*
*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
*        to compute its eigensystem.
*
         IF( m.EQ.l+1 ) THEN
            IF( icompz.GT.0 ) THEN
               CALL slaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
               work( l ) = c
               work( n-1+l ) = s
               CALL clasr( 'R', 'V', 'B', n, 2, work( l ),
     $                     work( n-1+l ), z( 1, l ), ldz )
            ELSE
               CALL slae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
            END IF
            d( l ) = rt1
            d( l+1 ) = rt2
            e( l ) = zero
            l = l + 2
            IF( l.LE.lend )
     $         GO TO 40
            GO TO 140
         END IF
*
         IF( jtot.EQ.nmaxit )
     $      GO TO 140
         jtot = jtot + 1
*
*        Form shift.
*
         g = ( d( l+1 )-p ) / ( two*e( l ) )
         r = slapy2( g, one )
         g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
*
         s = one
         c = one
         p = zero
*
*        Inner loop
*
         mm1 = m - 1
         DO 70 i = mm1, l, -1
            f = s*e( i )
            b = c*e( i )
            CALL slartg( g, f, c, s, r )
            IF( i.NE.m-1 )
     $         e( i+1 ) = r
            g = d( i+1 ) - p
            r = ( d( i )-g )*s + two*c*b
            p = s*r
            d( i+1 ) = g + p
            g = c*r - b
*
*           If eigenvectors are desired, then save rotations.
*
            IF( icompz.GT.0 ) THEN
               work( i ) = c
               work( n-1+i ) = -s
            END IF
*
   70    CONTINUE
*
*        If eigenvectors are desired, then apply saved rotations.
*
         IF( icompz.GT.0 ) THEN
            mm = m - l + 1
            CALL clasr( 'R', 'V', 'B', n, mm, work( l ), work( n-1+l ),
     $                  z( 1, l ), ldz )
         END IF
*
         d( l ) = d( l ) - p
         e( l ) = g
         GO TO 40
*
*        Eigenvalue found.
*
   80    CONTINUE
         d( l ) = p
*
         l = l + 1
         IF( l.LE.lend )
     $      GO TO 40
         GO TO 140
*
      ELSE
*
*        QR Iteration
*
*        Look for small superdiagonal element.
*
   90    CONTINUE
         IF( l.NE.lend ) THEN
            lendp1 = lend + 1
            DO 100 m = l, lendp1, -1
               tst = abs( e( m-1 ) )**2
               IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
     $             safmin )GO TO 110
  100       CONTINUE
         END IF
*
         m = lend
*
  110    CONTINUE
         IF( m.GT.lend )
     $      e( m-1 ) = zero
         p = d( l )
         IF( m.EQ.l )
     $      GO TO 130
*
*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
*        to compute its eigensystem.
*
         IF( m.EQ.l-1 ) THEN
            IF( icompz.GT.0 ) THEN
               CALL slaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
               work( m ) = c
               work( n-1+m ) = s
               CALL clasr( 'R', 'V', 'F', n, 2, work( m ),
     $                     work( n-1+m ), z( 1, l-1 ), ldz )
            ELSE
               CALL slae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
            END IF
            d( l-1 ) = rt1
            d( l ) = rt2
            e( l-1 ) = zero
            l = l - 2
            IF( l.GE.lend )
     $         GO TO 90
            GO TO 140
         END IF
*
         IF( jtot.EQ.nmaxit )
     $      GO TO 140
         jtot = jtot + 1
*
*        Form shift.
*
         g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
         r = slapy2( g, one )
         g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
*
         s = one
         c = one
         p = zero
*
*        Inner loop
*
         lm1 = l - 1
         DO 120 i = m, lm1
            f = s*e( i )
            b = c*e( i )
            CALL slartg( g, f, c, s, r )
            IF( i.NE.m )
     $         e( i-1 ) = r
            g = d( i ) - p
            r = ( d( i+1 )-g )*s + two*c*b
            p = s*r
            d( i ) = g + p
            g = c*r - b
*
*           If eigenvectors are desired, then save rotations.
*
            IF( icompz.GT.0 ) THEN
               work( i ) = c
               work( n-1+i ) = s
            END IF
*
  120    CONTINUE
*
*        If eigenvectors are desired, then apply saved rotations.
*
         IF( icompz.GT.0 ) THEN
            mm = l - m + 1
            CALL clasr( 'R', 'V', 'F', n, mm, work( m ), work( n-1+m ),
     $                  z( 1, m ), ldz )
         END IF
*
         d( l ) = d( l ) - p
         e( lm1 ) = g
         GO TO 90
*
*        Eigenvalue found.
*
  130    CONTINUE
         d( l ) = p
*
         l = l - 1
         IF( l.GE.lend )
     $      GO TO 90
         GO TO 140
*
      END IF
*
*     Undo scaling if necessary
*
  140 CONTINUE
      IF( iscale.EQ.1 ) THEN
         CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
     $                d( lsv ), n, info )
         CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
     $                n, info )
      ELSE IF( iscale.EQ.2 ) THEN
         CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
     $                d( lsv ), n, info )
         CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
     $                n, info )
      END IF
*
*     Check for no convergence to an eigenvalue after a total
*     of N*MAXIT iterations.
*
      IF( jtot.EQ.nmaxit ) THEN
         DO 150 i = 1, n - 1
            IF( e( i ).NE.zero )
     $         info = info + 1
  150    CONTINUE
         RETURN
      END IF
      GO TO 10
*
*     Order eigenvalues and eigenvectors.
*
  160 CONTINUE
      IF( icompz.EQ.0 ) THEN
*
*        Use Quick Sort
*
         CALL slasrt( 'I', n, d, info )
*
      ELSE
*
*        Use Selection Sort to minimize swaps of eigenvectors
*
         DO 180 ii = 2, n
            i = ii - 1
            k = i
            p = d( i )
            DO 170 j = ii, n
               IF( d( j ).LT.p ) THEN
                  k = j
                  p = d( j )
               END IF
  170       CONTINUE
            IF( k.NE.i ) THEN
               d( k ) = d( i )
               d( i ) = p
               CALL cswap( n, z( 1, i ), 1, z( 1, k ), 1 )
            END IF
  180    CONTINUE
      END IF
      RETURN
*
*     End of CSTEQR
*

Here is the call graph for this function:

Here is the caller graph for this function: