cget37.f

Go to the documentation of this file.

*> \brief \b CGET37
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE CGET37( RMAX, LMAX, NINFO, KNT, NIN )
* 
*       .. Scalar Arguments ..
*       INTEGER            KNT, NIN
*       ..
*       .. Array Arguments ..
*       INTEGER            LMAX( 3 ), NINFO( 3 )
*       REAL               RMAX( 3 )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CGET37 tests CTRSNA, a routine for estimating condition numbers of
*> eigenvalues and/or right eigenvectors of a matrix.
*>
*> The test matrices are read from a file with logical unit number NIN.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[out] RMAX
*> \verbatim
*>          RMAX is REAL array, dimension (3)
*>          Value of the largest test ratio.
*>          RMAX(1) = largest ratio comparing different calls to CTRSNA
*>          RMAX(2) = largest error in reciprocal condition
*>                    numbers taking their conditioning into account
*>          RMAX(3) = largest error in reciprocal condition
*>                    numbers not taking their conditioning into
*>                    account (may be larger than RMAX(2))
*> \endverbatim
*>
*> \param[out] LMAX
*> \verbatim
*>          LMAX is INTEGER array, dimension (3)
*>          LMAX(i) is example number where largest test ratio
*>          RMAX(i) is achieved. Also:
*>          If CGEHRD returns INFO nonzero on example i, LMAX(1)=i
*>          If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
*>          If CTRSNA returns INFO nonzero on example i, LMAX(3)=i
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*>          NINFO is INTEGER array, dimension (3)
*>          NINFO(1) = No. of times CGEHRD returned INFO nonzero
*>          NINFO(2) = No. of times CHSEQR returned INFO nonzero
*>          NINFO(3) = No. of times CTRSNA returned INFO nonzero
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*>          KNT is INTEGER
*>          Total number of examples tested.
*> \endverbatim
*>
*> \param[in] NIN
*> \verbatim
*>          NIN is INTEGER
*>          Input logical unit number
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup complex_eig
*
*  =====================================================================
      SUBROUTINE cget37( RMAX, LMAX, NINFO, KNT, NIN )
*
*  -- LAPACK test 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 ..
      INTEGER            KNT, NIN
*     ..
*     .. Array Arguments ..
      INTEGER            LMAX( 3 ), NINFO( 3 )
      REAL               RMAX( 3 )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO
      parameter                ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
      REAL               EPSIN
      parameter                ( epsin = 5.9605e-8 )
      INTEGER            LDT, LWORK
      parameter                ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ICMP, INFO, ISCL, ISRT, J, KMIN, M, N
      REAL               BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
     $                   vcmin, vmax, vmin, vmul
*     ..
*     .. Local Arrays ..
      LOGICAL            SELECT( ldt )
      INTEGER            LCMP( 3 )
      REAL               DUM( 1 ), RWORK( 2*ldt ), S( ldt ), SEP( ldt ),
     $                   sepin( ldt ), septmp( ldt ), sin( ldt ),
     $                   stmp( ldt ), val( 3 ), wiin( ldt ),
     $                   wrin( ldt ), wsrt( ldt )
      COMPLEX            CDUM( 1 ), LE( ldt, ldt ), RE( ldt, ldt ),
     $                   t( ldt, ldt ), tmp( ldt, ldt ), w( ldt ),
     $                   work( lwork ), wtmp( ldt )
*     ..
*     .. External Functions ..
      REAL               CLANGE, SLAMCH
      EXTERNAL           clange, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           ccopy, cgehrd, chseqr, clacpy, csscal, ctrevc,
     $                   ctrsna, scopy, slabad, sscal
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          aimag, max, REAL, SQRT
*     ..
*     .. Executable Statements ..
*
      eps = slamch( 'P' )
      smlnum = slamch( 'S' ) / eps
      bignum = one / smlnum
      CALL slabad( smlnum, bignum )
*
*     EPSIN = 2**(-24) = precision to which input data computed
*
      eps = max( eps, epsin )
      rmax( 1 ) = zero
      rmax( 2 ) = zero
      rmax( 3 ) = zero
      lmax( 1 ) = 0
      lmax( 2 ) = 0
      lmax( 3 ) = 0
      knt = 0
      ninfo( 1 ) = 0
      ninfo( 2 ) = 0
      ninfo( 3 ) = 0
      val( 1 ) = sqrt( smlnum )
      val( 2 ) = one
      val( 3 ) = sqrt( bignum )
*
*     Read input data until N=0.  Assume input eigenvalues are sorted
*     lexicographically (increasing by real part if ISRT = 0,
*     increasing by imaginary part if ISRT = 1)
*
   10 CONTINUE
      READ( nin, fmt = * )n, isrt
      IF( n.EQ.0 )
     $   RETURN
      DO 20 i = 1, n
         READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
   20 CONTINUE
      DO 30 i = 1, n
         READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
   30 CONTINUE
      tnrm = clange( 'M', n, n, tmp, ldt, rwork )
      DO 260 iscl = 1, 3
*
*        Scale input matrix
*
         knt = knt + 1
         CALL clacpy( 'F', n, n, tmp, ldt, t, ldt )
         vmul = val( iscl )
         DO 40 i = 1, n
            CALL csscal( n, vmul, t( 1, i ), 1 )
   40    CONTINUE
         IF( tnrm.EQ.zero )
     $      vmul = one
*
*        Compute eigenvalues and eigenvectors
*
         CALL cgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
     $                info )
         IF( info.NE.0 ) THEN
            lmax( 1 ) = knt
            ninfo( 1 ) = ninfo( 1 ) + 1
            GO TO 260
         END IF
         DO 60 j = 1, n - 2
            DO 50 i = j + 2, n
               t( i, j ) = zero
   50       CONTINUE
   60    CONTINUE
*
*        Compute Schur form
*
         CALL chseqr( 'S', 'N', n, 1, n, t, ldt, w, cdum, 1, work,
     $                lwork, info )
         IF( info.NE.0 ) THEN
            lmax( 2 ) = knt
            ninfo( 2 ) = ninfo( 2 ) + 1
            GO TO 260
         END IF
*
*        Compute eigenvectors
*
         DO 70 i = 1, n
            SELECT( i ) = .true.
   70    CONTINUE
         CALL ctrevc( 'B', 'A', SELECT, n, t, ldt, le, ldt, re, ldt, n,
     $                m, work, rwork, info )
*
*        Compute condition numbers
*
         CALL ctrsna( 'B', 'A', SELECT, n, t, ldt, le, ldt, re, ldt, s,
     $                sep, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
*
*        Sort eigenvalues and condition numbers lexicographically
*        to compare with inputs
*
         CALL ccopy( n, w, 1, wtmp, 1 )
         IF( isrt.EQ.0 ) THEN
*
*           Sort by increasing real part
*
            DO 80 i = 1, n
               wsrt( i ) = REAL( W( I ) )
   80       CONTINUE
         ELSE
*
*           Sort by increasing imaginary part
*
            DO 90 i = 1, n
               wsrt( i ) = aimag( w( i ) )
   90       CONTINUE
         END IF
         CALL scopy( n, s, 1, stmp, 1 )
         CALL scopy( n, sep, 1, septmp, 1 )
         CALL sscal( n, one / vmul, septmp, 1 )
         DO 110 i = 1, n - 1
            kmin = i
            vmin = wsrt( i )
            DO 100 j = i + 1, n
               IF( wsrt( j ).LT.vmin ) THEN
                  kmin = j
                  vmin = wsrt( j )
               END IF
  100       CONTINUE
            wsrt( kmin ) = wsrt( i )
            wsrt( i ) = vmin
            vcmin = wtmp( i )
            wtmp( i ) = w( kmin )
            wtmp( kmin ) = vcmin
            vmin = stmp( kmin )
            stmp( kmin ) = stmp( i )
            stmp( i ) = vmin
            vmin = septmp( kmin )
            septmp( kmin ) = septmp( i )
            septmp( i ) = vmin
  110    CONTINUE
*
*        Compare condition numbers for eigenvalues
*        taking their condition numbers into account
*
         v = max( two*REAL( n )*EPS*TNRM, SMLNUM )
         IF( tnrm.EQ.zero )
     $      v = one
         DO 120 i = 1, n
            IF( v.GT.septmp( i ) ) THEN
               tol = one
            ELSE
               tol = v / septmp( i )
            END IF
            IF( v.GT.sepin( i ) ) THEN
               tolin = one
            ELSE
               tolin = v / sepin( i )
            END IF
            tol = max( tol, smlnum / eps )
            tolin = max( tolin, smlnum / eps )
            IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol ) THEN
               vmax = one / eps
            ELSE IF( sin( i )-tolin.GT.stmp( i )+tol ) THEN
               vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
            ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) ) THEN
               vmax = one / eps
            ELSE IF( sin( i )+tolin.LT.stmp( i )-tol ) THEN
               vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
            ELSE
               vmax = one
            END IF
            IF( vmax.GT.rmax( 2 ) ) THEN
               rmax( 2 ) = vmax
               IF( ninfo( 2 ).EQ.0 )
     $            lmax( 2 ) = knt
            END IF
  120    CONTINUE
*
*        Compare condition numbers for eigenvectors
*        taking their condition numbers into account
*
         DO 130 i = 1, n
            IF( v.GT.septmp( i )*stmp( i ) ) THEN
               tol = septmp( i )
            ELSE
               tol = v / stmp( i )
            END IF
            IF( v.GT.sepin( i )*sin( i ) ) THEN
               tolin = sepin( i )
            ELSE
               tolin = v / sin( i )
            END IF
            tol = max( tol, smlnum / eps )
            tolin = max( tolin, smlnum / eps )
            IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol ) THEN
               vmax = one / eps
            ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol ) THEN
               vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
            ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) ) THEN
               vmax = one / eps
            ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol ) THEN
               vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
            ELSE
               vmax = one
            END IF
            IF( vmax.GT.rmax( 2 ) ) THEN
               rmax( 2 ) = vmax
               IF( ninfo( 2 ).EQ.0 )
     $            lmax( 2 ) = knt
            END IF
  130    CONTINUE
*
*        Compare condition numbers for eigenvalues
*        without taking their condition numbers into account
*
         DO 140 i = 1, n
            IF( sin( i ).LE.REAL( 2*n )*EPS .AND. STMP( i ).LE.
     $          REAL( 2*n )*EPS ) then
               vmax = one
            ELSE IF( eps*sin( i ).GT.stmp( i ) ) THEN
               vmax = one / eps
            ELSE IF( sin( i ).GT.stmp( i ) ) THEN
               vmax = sin( i ) / stmp( i )
            ELSE IF( sin( i ).LT.eps*stmp( i ) ) THEN
               vmax = one / eps
            ELSE IF( sin( i ).LT.stmp( i ) ) THEN
               vmax = stmp( i ) / sin( i )
            ELSE
               vmax = one
            END IF
            IF( vmax.GT.rmax( 3 ) ) THEN
               rmax( 3 ) = vmax
               IF( ninfo( 3 ).EQ.0 )
     $            lmax( 3 ) = knt
            END IF
  140    CONTINUE
*
*        Compare condition numbers for eigenvectors
*        without taking their condition numbers into account
*
         DO 150 i = 1, n
            IF( sepin( i ).LE.v .AND. septmp( i ).LE.v ) THEN
               vmax = one
            ELSE IF( eps*sepin( i ).GT.septmp( i ) ) THEN
               vmax = one / eps
            ELSE IF( sepin( i ).GT.septmp( i ) ) THEN
               vmax = sepin( i ) / septmp( i )
            ELSE IF( sepin( i ).LT.eps*septmp( i ) ) THEN
               vmax = one / eps
            ELSE IF( sepin( i ).LT.septmp( i ) ) THEN
               vmax = septmp( i ) / sepin( i )
            ELSE
               vmax = one
            END IF
            IF( vmax.GT.rmax( 3 ) ) THEN
               rmax( 3 ) = vmax
               IF( ninfo( 3 ).EQ.0 )
     $            lmax( 3 ) = knt
            END IF
  150    CONTINUE
*
*        Compute eigenvalue condition numbers only and compare
*
         vmax = zero
         dum( 1 ) = -one
         CALL scopy( n, dum, 0, stmp, 1 )
         CALL scopy( n, dum, 0, septmp, 1 )
         CALL ctrsna( 'E', 'A', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 160 i = 1, n
            IF( stmp( i ).NE.s( i ) )
     $         vmax = one / eps
            IF( septmp( i ).NE.dum( 1 ) )
     $         vmax = one / eps
  160    CONTINUE
*
*        Compute eigenvector condition numbers only and compare
*
         CALL scopy( n, dum, 0, stmp, 1 )
         CALL scopy( n, dum, 0, septmp, 1 )
         CALL ctrsna( 'V', 'A', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 170 i = 1, n
            IF( stmp( i ).NE.dum( 1 ) )
     $         vmax = one / eps
            IF( septmp( i ).NE.sep( i ) )
     $         vmax = one / eps
  170    CONTINUE
*
*        Compute all condition numbers using SELECT and compare
*
         DO 180 i = 1, n
            SELECT( i ) = .true.
  180    CONTINUE
         CALL scopy( n, dum, 0, stmp, 1 )
         CALL scopy( n, dum, 0, septmp, 1 )
         CALL ctrsna( 'B', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 190 i = 1, n
            IF( septmp( i ).NE.sep( i ) )
     $         vmax = one / eps
            IF( stmp( i ).NE.s( i ) )
     $         vmax = one / eps
  190    CONTINUE
*
*        Compute eigenvalue condition numbers using SELECT and compare
*
         CALL scopy( n, dum, 0, stmp, 1 )
         CALL scopy( n, dum, 0, septmp, 1 )
         CALL ctrsna( 'E', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 200 i = 1, n
            IF( stmp( i ).NE.s( i ) )
     $         vmax = one / eps
            IF( septmp( i ).NE.dum( 1 ) )
     $         vmax = one / eps
  200    CONTINUE
*
*        Compute eigenvector condition numbers using SELECT and compare
*
         CALL scopy( n, dum, 0, stmp, 1 )
         CALL scopy( n, dum, 0, septmp, 1 )
         CALL ctrsna( 'V', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 210 i = 1, n
            IF( stmp( i ).NE.dum( 1 ) )
     $         vmax = one / eps
            IF( septmp( i ).NE.sep( i ) )
     $         vmax = one / eps
  210    CONTINUE
         IF( vmax.GT.rmax( 1 ) ) THEN
            rmax( 1 ) = vmax
            IF( ninfo( 1 ).EQ.0 )
     $         lmax( 1 ) = knt
         END IF
*
*        Select second and next to last eigenvalues
*
         DO 220 i = 1, n
            SELECT( i ) = .false.
  220    CONTINUE
         icmp = 0
         IF( n.GT.1 ) THEN
            icmp = 1
            lcmp( 1 ) = 2
            SELECT( 2 ) = .true.
            CALL ccopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
            CALL ccopy( n, le( 1, 2 ), 1, le( 1, 1 ), 1 )
         END IF
         IF( n.GT.3 ) THEN
            icmp = 2
            lcmp( 2 ) = n - 1
            SELECT( n-1 ) = .true.
            CALL ccopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
            CALL ccopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
         END IF
*
*        Compute all selected condition numbers
*
         CALL scopy( icmp, dum, 0, stmp, 1 )
         CALL scopy( icmp, dum, 0, septmp, 1 )
         CALL ctrsna( 'B', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 230 i = 1, icmp
            j = lcmp( i )
            IF( septmp( i ).NE.sep( j ) )
     $         vmax = one / eps
            IF( stmp( i ).NE.s( j ) )
     $         vmax = one / eps
  230    CONTINUE
*
*        Compute selected eigenvalue condition numbers
*
         CALL scopy( icmp, dum, 0, stmp, 1 )
         CALL scopy( icmp, dum, 0, septmp, 1 )
         CALL ctrsna( 'E', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 240 i = 1, icmp
            j = lcmp( i )
            IF( stmp( i ).NE.s( j ) )
     $         vmax = one / eps
            IF( septmp( i ).NE.dum( 1 ) )
     $         vmax = one / eps
  240    CONTINUE
*
*        Compute selected eigenvector condition numbers
*
         CALL scopy( icmp, dum, 0, stmp, 1 )
         CALL scopy( icmp, dum, 0, septmp, 1 )
         CALL ctrsna( 'V', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
     $                stmp, septmp, n, m, work, n, rwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 260
         END IF
         DO 250 i = 1, icmp
            j = lcmp( i )
            IF( stmp( i ).NE.dum( 1 ) )
     $         vmax = one / eps
            IF( septmp( i ).NE.sep( j ) )
     $         vmax = one / eps
  250    CONTINUE
         IF( vmax.GT.rmax( 1 ) ) THEN
            rmax( 1 ) = vmax
            IF( ninfo( 1 ).EQ.0 )
     $         lmax( 1 ) = knt
         END IF
  260 CONTINUE
      GO TO 10
*
*     End of CGET37
*
      END