cget37.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 537 行 · 第 1/2 页

      SUBROUTINE CGET37( RMAX, LMAX, NINFO, KNT, NIN )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            KNT, NIN
*     ..
*     .. Array Arguments ..
      INTEGER            LMAX( 3 ), NINFO( 3 )
      REAL               RMAX( 3 )
*     ..
*
*  Purpose
*  =======
*
*  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.
*
*  Arguments
*  ==========
*
*  RMAX    (output) 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))
*
*  LMAX    (output) 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
*
*  NINFO   (output) 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
*
*  KNT     (output) INTEGER
*          Total number of examples tested.
*
*  NIN     (input) INTEGER
*          Input logical unit number
*
*  =====================================================================
*
*     .. 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