zget23.f

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

*     Do Test (1)
*
      CALL ZGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, W, WORK, RWORK,
     $             RES )
      RESULT( 1 ) = RES( 1 )
*
*     Do Test (2)
*
      CALL ZGET22( 'C', 'N', 'C', N, A, LDA, VL, LDVL, W, WORK, RWORK,
     $             RES )
      RESULT( 2 ) = RES( 1 )
*
*     Do Test (3)
*
      DO 30 J = 1, N
         TNRM = DZNRM2( N, VR( 1, J ), 1 )
         RESULT( 3 ) = MAX( RESULT( 3 ),
     $                 MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
         VMX = ZERO
         VRMX = ZERO
         DO 20 JJ = 1, N
            VTST = ABS( VR( JJ, J ) )
            IF( VTST.GT.VMX )
     $         VMX = VTST
            IF( DIMAG( VR( JJ, J ) ).EQ.ZERO .AND.
     $          ABS( DBLE( VR( JJ, J ) ) ).GT.VRMX )
     $          VRMX = ABS( DBLE( VR( JJ, J ) ) )
   20    CONTINUE
         IF( VRMX / VMX.LT.ONE-TWO*ULP )
     $      RESULT( 3 ) = ULPINV
   30 CONTINUE
*
*     Do Test (4)
*
      DO 50 J = 1, N
         TNRM = DZNRM2( N, VL( 1, J ), 1 )
         RESULT( 4 ) = MAX( RESULT( 4 ),
     $                 MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
         VMX = ZERO
         VRMX = ZERO
         DO 40 JJ = 1, N
            VTST = ABS( VL( JJ, J ) )
            IF( VTST.GT.VMX )
     $         VMX = VTST
            IF( DIMAG( VL( JJ, J ) ).EQ.ZERO .AND.
     $          ABS( DBLE( VL( JJ, J ) ) ).GT.VRMX )
     $          VRMX = ABS( DBLE( VL( JJ, J ) ) )
   40    CONTINUE
         IF( VRMX / VMX.LT.ONE-TWO*ULP )
     $      RESULT( 4 ) = ULPINV
   50 CONTINUE
*
*     Test for all options of computing condition numbers
*
      DO 200 ISENS = 1, ISENSM
*
         SENSE = SENS( ISENS )
*
*        Compute eigenvalues only, and test them
*
         CALL ZLACPY( 'F', N, N, A, LDA, H, LDA )
         CALL ZGEEVX( BALANC, 'N', 'N', SENSE, N, H, LDA, W1, CDUM, 1,
     $                CDUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
     $                RCNDV1, WORK, LWORK, RWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 1 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEEVX2', IINFO, N, JTYPE,
     $            BALANC, ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEEVX2', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 190
         END IF
*
*        Do Test (5)
*
         DO 60 J = 1, N
            IF( W( J ).NE.W1( J ) )
     $         RESULT( 5 ) = ULPINV
   60    CONTINUE
*
*        Do Test (8)
*
         IF( .NOT.NOBAL ) THEN
            DO 70 J = 1, N
               IF( SCALE( J ).NE.SCALE1( J ) )
     $            RESULT( 8 ) = ULPINV
   70       CONTINUE
            IF( ILO.NE.ILO1 )
     $         RESULT( 8 ) = ULPINV
            IF( IHI.NE.IHI1 )
     $         RESULT( 8 ) = ULPINV
            IF( ABNRM.NE.ABNRM1 )
     $         RESULT( 8 ) = ULPINV
         END IF
*
*        Do Test (9)
*
         IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
            DO 80 J = 1, N
               IF( RCONDV( J ).NE.RCNDV1( J ) )
     $            RESULT( 9 ) = ULPINV
   80       CONTINUE
         END IF
*
*        Compute eigenvalues and right eigenvectors, and test them
*
         CALL ZLACPY( 'F', N, N, A, LDA, H, LDA )
         CALL ZGEEVX( BALANC, 'N', 'V', SENSE, N, H, LDA, W1, CDUM, 1,
     $                LRE, LDLRE, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
     $                RCNDV1, WORK, LWORK, RWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 1 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEEVX3', IINFO, N, JTYPE,
     $            BALANC, ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEEVX3', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 190
         END IF
*
*        Do Test (5) again
*
         DO 90 J = 1, N
            IF( W( J ).NE.W1( J ) )
     $         RESULT( 5 ) = ULPINV
   90    CONTINUE
*
*        Do Test (6)
*
         DO 110 J = 1, N
            DO 100 JJ = 1, N
               IF( VR( J, JJ ).NE.LRE( J, JJ ) )
     $            RESULT( 6 ) = ULPINV
  100       CONTINUE
  110    CONTINUE
*
*        Do Test (8) again
*
         IF( .NOT.NOBAL ) THEN
            DO 120 J = 1, N
               IF( SCALE( J ).NE.SCALE1( J ) )
     $            RESULT( 8 ) = ULPINV
  120       CONTINUE
            IF( ILO.NE.ILO1 )
     $         RESULT( 8 ) = ULPINV
            IF( IHI.NE.IHI1 )
     $         RESULT( 8 ) = ULPINV
            IF( ABNRM.NE.ABNRM1 )
     $         RESULT( 8 ) = ULPINV
         END IF
*
*        Do Test (9) again
*
         IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
            DO 130 J = 1, N
               IF( RCONDV( J ).NE.RCNDV1( J ) )
     $            RESULT( 9 ) = ULPINV
  130       CONTINUE
         END IF
*
*        Compute eigenvalues and left eigenvectors, and test them
*
         CALL ZLACPY( 'F', N, N, A, LDA, H, LDA )
         CALL ZGEEVX( BALANC, 'V', 'N', SENSE, N, H, LDA, W1, LRE,
     $                LDLRE, CDUM, 1, ILO1, IHI1, SCALE1, ABNRM1,
     $                RCNDE1, RCNDV1, WORK, LWORK, RWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 1 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEEVX4', IINFO, N, JTYPE,
     $            BALANC, ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEEVX4', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 190
         END IF
*
*        Do Test (5) again
*
         DO 140 J = 1, N
            IF( W( J ).NE.W1( J ) )
     $         RESULT( 5 ) = ULPINV
  140    CONTINUE
*
*        Do Test (7)
*
         DO 160 J = 1, N
            DO 150 JJ = 1, N
               IF( VL( J, JJ ).NE.LRE( J, JJ ) )
     $            RESULT( 7 ) = ULPINV
  150       CONTINUE
  160    CONTINUE
*
*        Do Test (8) again
*
         IF( .NOT.NOBAL ) THEN
            DO 170 J = 1, N
               IF( SCALE( J ).NE.SCALE1( J ) )
     $            RESULT( 8 ) = ULPINV
  170       CONTINUE
            IF( ILO.NE.ILO1 )
     $         RESULT( 8 ) = ULPINV
            IF( IHI.NE.IHI1 )
     $         RESULT( 8 ) = ULPINV
            IF( ABNRM.NE.ABNRM1 )
     $         RESULT( 8 ) = ULPINV
         END IF
*
*        Do Test (9) again
*
         IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
            DO 180 J = 1, N
               IF( RCONDV( J ).NE.RCNDV1( J ) )
     $            RESULT( 9 ) = ULPINV
  180       CONTINUE
         END IF
*
  190    CONTINUE
*
  200 CONTINUE
*
*     If COMP, compare condition numbers to precomputed ones
*
      IF( COMP ) THEN
         CALL ZLACPY( 'F', N, N, A, LDA, H, LDA )
         CALL ZGEEVX( 'N', 'V', 'V', 'B', N, H, LDA, W, VL, LDVL, VR,
     $                LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV,
     $                WORK, LWORK, RWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 1 ) = ULPINV
            WRITE( NOUNIT, FMT = 9999 )'ZGEEVX5', IINFO, N, ISEED( 1 )
            INFO = ABS( IINFO )
            GO TO 250
         END IF
*
*        Sort eigenvalues and condition numbers lexicographically
*        to compare with inputs
*
         DO 220 I = 1, N - 1
            KMIN = I
            IF( ISRT.EQ.0 ) THEN
               VRIMIN = DBLE( W( I ) )
            ELSE
               VRIMIN = DIMAG( W( I ) )
            END IF
            DO 210 J = I + 1, N
               IF( ISRT.EQ.0 ) THEN
                  VRICMP = DBLE( W( J ) )
               ELSE
                  VRICMP = DIMAG( W( J ) )
               END IF
               IF( VRICMP.LT.VRIMIN ) THEN
                  KMIN = J
                  VRIMIN = VRICMP
               END IF
  210       CONTINUE
            CTMP = W( KMIN )
            W( KMIN ) = W( I )
            W( I ) = CTMP
            VRIMIN = RCONDE( KMIN )
            RCONDE( KMIN ) = RCONDE( I )
            RCONDE( I ) = VRIMIN
            VRIMIN = RCONDV( KMIN )
            RCONDV( KMIN ) = RCONDV( I )
            RCONDV( I ) = VRIMIN
  220    CONTINUE
*
*        Compare condition numbers for eigenvectors
*        taking their condition numbers into account
*
         RESULT( 10 ) = ZERO
         EPS = MAX( EPSIN, ULP )
         V = MAX( DBLE( N )*EPS*ABNRM, SMLNUM )
         IF( ABNRM.EQ.ZERO )
     $      V = ONE
         DO 230 I = 1, N
            IF( V.GT.RCONDV( I )*RCONDE( I ) ) THEN
               TOL = RCONDV( I )
            ELSE
               TOL = V / RCONDE( I )
            END IF
            IF( V.GT.RCDVIN( I )*RCDEIN( I ) ) THEN
               TOLIN = RCDVIN( I )
            ELSE
               TOLIN = V / RCDEIN( I )
            END IF
            TOL = MAX( TOL, SMLNUM / EPS )
            TOLIN = MAX( TOLIN, SMLNUM / EPS )
            IF( EPS*( RCDVIN( I )-TOLIN ).GT.RCONDV( I )+TOL ) THEN
               VMAX = ONE / EPS
            ELSE IF( RCDVIN( I )-TOLIN.GT.RCONDV( I )+TOL ) THEN
               VMAX = ( RCDVIN( I )-TOLIN ) / ( RCONDV( I )+TOL )
            ELSE IF( RCDVIN( I )+TOLIN.LT.EPS*( RCONDV( I )-TOL ) ) THEN
               VMAX = ONE / EPS
            ELSE IF( RCDVIN( I )+TOLIN.LT.RCONDV( I )-TOL ) THEN
               VMAX = ( RCONDV( I )-TOL ) / ( RCDVIN( I )+TOLIN )
            ELSE
               VMAX = ONE
            END IF
            RESULT( 10 ) = MAX( RESULT( 10 ), VMAX )
  230    CONTINUE
*
*        Compare condition numbers for eigenvalues
*        taking their condition numbers into account
*
         RESULT( 11 ) = ZERO
         DO 240 I = 1, N
            IF( V.GT.RCONDV( I ) ) THEN
               TOL = ONE
            ELSE
               TOL = V / RCONDV( I )
            END IF
            IF( V.GT.RCDVIN( I ) ) THEN
               TOLIN = ONE
            ELSE
               TOLIN = V / RCDVIN( I )
            END IF
            TOL = MAX( TOL, SMLNUM / EPS )
            TOLIN = MAX( TOLIN, SMLNUM / EPS )
            IF( EPS*( RCDEIN( I )-TOLIN ).GT.RCONDE( I )+TOL ) THEN
               VMAX = ONE / EPS
            ELSE IF( RCDEIN( I )-TOLIN.GT.RCONDE( I )+TOL ) THEN
               VMAX = ( RCDEIN( I )-TOLIN ) / ( RCONDE( I )+TOL )
            ELSE IF( RCDEIN( I )+TOLIN.LT.EPS*( RCONDE( I )-TOL ) ) THEN
               VMAX = ONE / EPS
            ELSE IF( RCDEIN( I )+TOLIN.LT.RCONDE( I )-TOL ) THEN
               VMAX = ( RCONDE( I )-TOL ) / ( RCDEIN( I )+TOLIN )
            ELSE
               VMAX = ONE
            END IF
            RESULT( 11 ) = MAX( RESULT( 11 ), VMAX )
  240    CONTINUE
  250    CONTINUE
*
      END IF
*
 9999 FORMAT( ' ZGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', INPUT EXAMPLE NUMBER = ', I4 )
 9998 FORMAT( ' ZGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', BALANC = ', A, ', ISEED=(',
     $      3( I5, ',' ), I5, ')' )
*
      RETURN
*
*     End of ZGET23
*
      END