zget24.f

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

         END IF
*
         RESULT( 5+RSUB ) = ZERO
         DO 60 J = 1, N
            DO 50 I = 1, N
               IF( H( I, J ).NE.HT( I, J ) )
     $            RESULT( 5+RSUB ) = ULPINV
   50       CONTINUE
   60    CONTINUE
*
*        Do Test (6) or Test (12)
*
         RESULT( 6+RSUB ) = ZERO
         DO 70 I = 1, N
            IF( W( I ).NE.WT( I ) )
     $         RESULT( 6+RSUB ) = ULPINV
   70    CONTINUE
*
*        Do Test (13)
*
         IF( ISORT.EQ.1 ) THEN
            RESULT( 13 ) = ZERO
            KNTEIG = 0
            DO 80 I = 1, N
               IF( ZSLECT( W( I ) ) )
     $            KNTEIG = KNTEIG + 1
               IF( I.LT.N ) THEN
                  IF( ZSLECT( W( I+1 ) ) .AND.
     $                ( .NOT.ZSLECT( W( I ) ) ) )RESULT( 13 ) = ULPINV
               END IF
   80       CONTINUE
            IF( SDIM.NE.KNTEIG )
     $         RESULT( 13 ) = ULPINV
         END IF
*
   90 CONTINUE
*
*     If there is enough workspace, perform tests (14) and (15)
*     as well as (10) through (13)
*
      IF( LWORK.GE.( N*( N+1 ) ) / 2 ) THEN
*
*        Compute both RCONDE and RCONDV with VS
*
         SORT = 'S'
         RESULT( 14 ) = ZERO
         RESULT( 15 ) = ZERO
         CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
         CALL ZGEESX( 'V', SORT, ZSLECT, 'B', N, HT, LDA, SDIM1, WT,
     $                VS1, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
     $                BWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 14 ) = ULPINV
            RESULT( 15 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEESX3', IINFO, N, JTYPE,
     $            ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEESX3', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 220
         END IF
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 110 I = 1, N
            IF( W( I ).NE.WT( I ) )
     $         RESULT( 10 ) = ULPINV
            DO 100 J = 1, N
               IF( H( I, J ).NE.HT( I, J ) )
     $            RESULT( 11 ) = ULPINV
               IF( VS( I, J ).NE.VS1( I, J ) )
     $            RESULT( 12 ) = ULPINV
  100       CONTINUE
  110    CONTINUE
         IF( SDIM.NE.SDIM1 )
     $      RESULT( 13 ) = ULPINV
*
*        Compute both RCONDE and RCONDV without VS, and compare
*
         CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
         CALL ZGEESX( 'N', SORT, ZSLECT, 'B', N, HT, LDA, SDIM1, WT,
     $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
     $                BWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 14 ) = ULPINV
            RESULT( 15 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEESX4', IINFO, N, JTYPE,
     $            ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEESX4', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 220
         END IF
*
*        Perform tests (14) and (15)
*
         IF( RCNDE1.NE.RCONDE )
     $      RESULT( 14 ) = ULPINV
         IF( RCNDV1.NE.RCONDV )
     $      RESULT( 15 ) = ULPINV
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 130 I = 1, N
            IF( W( I ).NE.WT( I ) )
     $         RESULT( 10 ) = ULPINV
            DO 120 J = 1, N
               IF( H( I, J ).NE.HT( I, J ) )
     $            RESULT( 11 ) = ULPINV
               IF( VS( I, J ).NE.VS1( I, J ) )
     $            RESULT( 12 ) = ULPINV
  120       CONTINUE
  130    CONTINUE
         IF( SDIM.NE.SDIM1 )
     $      RESULT( 13 ) = ULPINV
*
*        Compute RCONDE with VS, and compare
*
         CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
         CALL ZGEESX( 'V', SORT, ZSLECT, 'E', N, HT, LDA, SDIM1, WT,
     $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
     $                BWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 14 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEESX5', IINFO, N, JTYPE,
     $            ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEESX5', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 220
         END IF
*
*        Perform test (14)
*
         IF( RCNDE1.NE.RCONDE )
     $      RESULT( 14 ) = ULPINV
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 150 I = 1, N
            IF( W( I ).NE.WT( I ) )
     $         RESULT( 10 ) = ULPINV
            DO 140 J = 1, N
               IF( H( I, J ).NE.HT( I, J ) )
     $            RESULT( 11 ) = ULPINV
               IF( VS( I, J ).NE.VS1( I, J ) )
     $            RESULT( 12 ) = ULPINV
  140       CONTINUE
  150    CONTINUE
         IF( SDIM.NE.SDIM1 )
     $      RESULT( 13 ) = ULPINV
*
*        Compute RCONDE without VS, and compare
*
         CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
         CALL ZGEESX( 'N', SORT, ZSLECT, 'E', N, HT, LDA, SDIM1, WT,
     $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
     $                BWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 14 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEESX6', IINFO, N, JTYPE,
     $            ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEESX6', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 220
         END IF
*
*        Perform test (14)
*
         IF( RCNDE1.NE.RCONDE )
     $      RESULT( 14 ) = ULPINV
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 170 I = 1, N
            IF( W( I ).NE.WT( I ) )
     $         RESULT( 10 ) = ULPINV
            DO 160 J = 1, N
               IF( H( I, J ).NE.HT( I, J ) )
     $            RESULT( 11 ) = ULPINV
               IF( VS( I, J ).NE.VS1( I, J ) )
     $            RESULT( 12 ) = ULPINV
  160       CONTINUE
  170    CONTINUE
         IF( SDIM.NE.SDIM1 )
     $      RESULT( 13 ) = ULPINV
*
*        Compute RCONDV with VS, and compare
*
         CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
         CALL ZGEESX( 'V', SORT, ZSLECT, 'V', N, HT, LDA, SDIM1, WT,
     $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
     $                BWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 15 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEESX7', IINFO, N, JTYPE,
     $            ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEESX7', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 220
         END IF
*
*        Perform test (15)
*
         IF( RCNDV1.NE.RCONDV )
     $      RESULT( 15 ) = ULPINV
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 190 I = 1, N
            IF( W( I ).NE.WT( I ) )
     $         RESULT( 10 ) = ULPINV
            DO 180 J = 1, N
               IF( H( I, J ).NE.HT( I, J ) )
     $            RESULT( 11 ) = ULPINV
               IF( VS( I, J ).NE.VS1( I, J ) )
     $            RESULT( 12 ) = ULPINV
  180       CONTINUE
  190    CONTINUE
         IF( SDIM.NE.SDIM1 )
     $      RESULT( 13 ) = ULPINV
*
*        Compute RCONDV without VS, and compare
*
         CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
         CALL ZGEESX( 'N', SORT, ZSLECT, 'V', N, HT, LDA, SDIM1, WT,
     $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
     $                BWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 15 ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'ZGEESX8', IINFO, N, JTYPE,
     $            ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'ZGEESX8', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            GO TO 220
         END IF
*
*        Perform test (15)
*
         IF( RCNDV1.NE.RCONDV )
     $      RESULT( 15 ) = ULPINV
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 210 I = 1, N
            IF( W( I ).NE.WT( I ) )
     $         RESULT( 10 ) = ULPINV
            DO 200 J = 1, N
               IF( H( I, J ).NE.HT( I, J ) )
     $            RESULT( 11 ) = ULPINV
               IF( VS( I, J ).NE.VS1( I, J ) )
     $            RESULT( 12 ) = ULPINV
  200       CONTINUE
  210    CONTINUE
         IF( SDIM.NE.SDIM1 )
     $      RESULT( 13 ) = ULPINV
*
      END IF
*
  220 CONTINUE
*
*     If there are precomputed reciprocal condition numbers, compare
*     computed values with them.
*
      IF( COMP ) THEN
*
*        First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that
*        the logical function ZSLECT selects the eigenvalues specified
*        by NSLCT, ISLCT and ISRT.
*
         SELDIM = N
         SELOPT = 1
         EPS = MAX( ULP, EPSIN )
         DO 230 I = 1, N
            IPNT( I ) = I
            SELVAL( I ) = .FALSE.
            SELWR( I ) = DBLE( WTMP( I ) )
            SELWI( I ) = DIMAG( WTMP( I ) )
  230    CONTINUE
         DO 250 I = 1, N - 1
            KMIN = I
            IF( ISRT.EQ.0 ) THEN
               VRIMIN = DBLE( WTMP( I ) )
            ELSE
               VRIMIN = DIMAG( WTMP( I ) )
            END IF
            DO 240 J = I + 1, N
               IF( ISRT.EQ.0 ) THEN
                  VRICMP = DBLE( WTMP( J ) )
               ELSE
                  VRICMP = DIMAG( WTMP( J ) )
               END IF
               IF( VRICMP.LT.VRIMIN ) THEN
                  KMIN = J
                  VRIMIN = VRICMP
               END IF
  240       CONTINUE
            CTMP = WTMP( KMIN )
            WTMP( KMIN ) = WTMP( I )
            WTMP( I ) = CTMP
            ITMP = IPNT( I )
            IPNT( I ) = IPNT( KMIN )
            IPNT( KMIN ) = ITMP
  250    CONTINUE
         DO 260 I = 1, NSLCT
            SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE.
  260    CONTINUE
*
*        Compute condition numbers
*
         CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
         CALL ZGEESX( 'N', 'S', ZSLECT, 'B', N, HT, LDA, SDIM1, WT, VS1,
     $                LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
     $                IINFO )
         IF( IINFO.NE.0 ) THEN
            RESULT( 16 ) = ULPINV
            RESULT( 17 ) = ULPINV
            WRITE( NOUNIT, FMT = 9999 )'ZGEESX9', IINFO, N, ISEED( 1 )
            INFO = ABS( IINFO )
            GO TO 270
         END IF
*
*        Compare condition number for average of selected eigenvalues
*        taking its condition number into account
*
         ANORM = ZLANGE( '1', N, N, A, LDA, RWORK )
         V = MAX( DBLE( N )*EPS*ANORM, SMLNUM )
         IF( ANORM.EQ.ZERO )
     $      V = ONE
         IF( V.GT.RCONDV ) THEN
            TOL = ONE
         ELSE
            TOL = V / RCONDV
         END IF
         IF( V.GT.RCDVIN ) THEN
            TOLIN = ONE
         ELSE
            TOLIN = V / RCDVIN
         END IF
         TOL = MAX( TOL, SMLNUM / EPS )
         TOLIN = MAX( TOLIN, SMLNUM / EPS )
         IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN
            RESULT( 16 ) = ULPINV
         ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN
            RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
         ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN
            RESULT( 16 ) = ULPINV
         ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN
            RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN )
         ELSE
            RESULT( 16 ) = ONE
         END IF
*
*        Compare condition numbers for right invariant subspace
*        taking its condition number into account
*
         IF( V.GT.RCONDV*RCONDE ) THEN
            TOL = RCONDV
         ELSE
            TOL = V / RCONDE
         END IF
         IF( V.GT.RCDVIN*RCDEIN ) THEN
            TOLIN = RCDVIN
         ELSE
            TOLIN = V / RCDEIN
         END IF
         TOL = MAX( TOL, SMLNUM / EPS )
         TOLIN = MAX( TOLIN, SMLNUM / EPS )
         IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN
            RESULT( 17 ) = ULPINV
         ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN
            RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
         ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN
            RESULT( 17 ) = ULPINV
         ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN
            RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN )
         ELSE
            RESULT( 17 ) = ONE
         END IF
*
  270    CONTINUE
*
      END IF
*
 9999 FORMAT( ' ZGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', INPUT EXAMPLE NUMBER = ', I4 )
 9998 FORMAT( ' ZGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
*
      RETURN
*
*     End of ZGET24
*
      END