dget24.f

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

      SUBROUTINE DGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
     $                   H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
     $                   LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
     $                   RESULT, WORK, LWORK, IWORK, BWORK, INFO )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      LOGICAL            COMP
      INTEGER            INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
      DOUBLE PRECISION   RCDEIN, RCDVIN, THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * )
      INTEGER            ISEED( 4 ), ISLCT( * ), IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), H( LDA, * ), HT( LDA, * ),
     $                   RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
     $                   WI( * ), WIT( * ), WITMP( * ), WORK( * ),
     $                   WR( * ), WRT( * ), WRTMP( * )
*     ..
*
*  Purpose
*  =======
*
*     DGET24 checks the nonsymmetric eigenvalue (Schur form) problem
*     expert driver DGEESX.
*
*     If COMP = .FALSE., the first 13 of the following tests will be
*     be performed on the input matrix A, and also tests 14 and 15
*     if LWORK is sufficiently large.
*     If COMP = .TRUE., all 17 test will be performed.
*
*     (1)     0 if T is in Schur form, 1/ulp otherwise
*            (no sorting of eigenvalues)
*
*     (2)     | A - VS T VS' | / ( n |A| ulp )
*
*       Here VS is the matrix of Schur eigenvectors, and T is in Schur
*       form  (no sorting of eigenvalues).
*
*     (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
*
*     (4)     0     if WR+sqrt(-1)*WI are eigenvalues of T
*             1/ulp otherwise
*             (no sorting of eigenvalues)
*
*     (5)     0     if T(with VS) = T(without VS),
*             1/ulp otherwise
*             (no sorting of eigenvalues)
*
*     (6)     0     if eigenvalues(with VS) = eigenvalues(without VS),
*             1/ulp otherwise
*             (no sorting of eigenvalues)
*
*     (7)     0 if T is in Schur form, 1/ulp otherwise
*             (with sorting of eigenvalues)
*
*     (8)     | A - VS T VS' | / ( n |A| ulp )
*
*       Here VS is the matrix of Schur eigenvectors, and T is in Schur
*       form  (with sorting of eigenvalues).
*
*     (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
*
*     (10)    0     if WR+sqrt(-1)*WI are eigenvalues of T
*             1/ulp otherwise
*             If workspace sufficient, also compare WR, WI with and
*             without reciprocal condition numbers
*             (with sorting of eigenvalues)
*
*     (11)    0     if T(with VS) = T(without VS),
*             1/ulp otherwise
*             If workspace sufficient, also compare T with and without
*             reciprocal condition numbers
*             (with sorting of eigenvalues)
*
*     (12)    0     if eigenvalues(with VS) = eigenvalues(without VS),
*             1/ulp otherwise
*             If workspace sufficient, also compare VS with and without
*             reciprocal condition numbers
*             (with sorting of eigenvalues)
*
*     (13)    if sorting worked and SDIM is the number of
*             eigenvalues which were SELECTed
*             If workspace sufficient, also compare SDIM with and
*             without reciprocal condition numbers
*
*     (14)    if RCONDE the same no matter if VS and/or RCONDV computed
*
*     (15)    if RCONDV the same no matter if VS and/or RCONDE computed
*
*     (16)  |RCONDE - RCDEIN| / cond(RCONDE)
*
*        RCONDE is the reciprocal average eigenvalue condition number
*        computed by DGEESX and RCDEIN (the precomputed true value)
*        is supplied as input.  cond(RCONDE) is the condition number
*        of RCONDE, and takes errors in computing RCONDE into account,
*        so that the resulting quantity should be O(ULP). cond(RCONDE)
*        is essentially given by norm(A)/RCONDV.
*
*     (17)  |RCONDV - RCDVIN| / cond(RCONDV)
*
*        RCONDV is the reciprocal right invariant subspace condition
*        number computed by DGEESX and RCDVIN (the precomputed true
*        value) is supplied as input. cond(RCONDV) is the condition
*        number of RCONDV, and takes errors in computing RCONDV into
*        account, so that the resulting quantity should be O(ULP).
*        cond(RCONDV) is essentially given by norm(A)/RCONDE.
*
*  Arguments
*  =========
*
*  COMP    (input) LOGICAL
*          COMP describes which input tests to perform:
*            = .FALSE. if the computed condition numbers are not to
*                      be tested against RCDVIN and RCDEIN
*            = .TRUE.  if they are to be compared
*
*  JTYPE   (input) INTEGER
*          Type of input matrix. Used to label output if error occurs.
*
*  ISEED   (input) INTEGER array, dimension (4)
*          If COMP = .FALSE., the random number generator seed
*          used to produce matrix.
*          If COMP = .TRUE., ISEED(1) = the number of the example.
*          Used to label output if error occurs.
*
*  THRESH  (input) DOUBLE PRECISION
*          A test will count as "failed" if the "error", computed as
*          described above, exceeds THRESH.  Note that the error
*          is scaled to be O(1), so THRESH should be a reasonably
*          small multiple of 1, e.g., 10 or 100.  In particular,
*          it should not depend on the precision (single vs. double)
*          or the size of the matrix.  It must be at least zero.
*
*  NOUNIT  (input) INTEGER
*          The FORTRAN unit number for printing out error messages
*          (e.g., if a routine returns INFO not equal to 0.)
*
*  N       (input) INTEGER
*          The dimension of A. N must be at least 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
*          Used to hold the matrix whose eigenvalues are to be
*          computed.
*
*  LDA     (input) INTEGER
*          The leading dimension of A, and H. LDA must be at
*          least 1 and at least N.
*
*  H       (workspace) DOUBLE PRECISION array, dimension (LDA, N)
*          Another copy of the test matrix A, modified by DGEESX.
*
*  HT      (workspace) DOUBLE PRECISION array, dimension (LDA, N)
*          Yet another copy of the test matrix A, modified by DGEESX.
*
*  WR      (workspace) DOUBLE PRECISION array, dimension (N)
*  WI      (workspace) DOUBLE PRECISION array, dimension (N)
*          The real and imaginary parts of the eigenvalues of A.
*          On exit, WR + WI*i are the eigenvalues of the matrix in A.
*
*  WRT     (workspace) DOUBLE PRECISION array, dimension (N)
*  WIT     (workspace) DOUBLE PRECISION array, dimension (N)
*          Like WR, WI, these arrays contain the eigenvalues of A,
*          but those computed when DGEESX only computes a partial
*          eigendecomposition, i.e. not Schur vectors
*
*  WRTMP   (workspace) DOUBLE PRECISION array, dimension (N)
*  WITMP   (workspace) DOUBLE PRECISION array, dimension (N)
*          Like WR, WI, these arrays contain the eigenvalues of A,
*          but sorted by increasing real part.
*
*  VS      (workspace) DOUBLE PRECISION array, dimension (LDVS, N)
*          VS holds the computed Schur vectors.
*
*  LDVS    (input) INTEGER
*          Leading dimension of VS. Must be at least max(1, N).
*
*  VS1     (workspace) DOUBLE PRECISION array, dimension (LDVS, N)
*          VS1 holds another copy of the computed Schur vectors.
*
*  RCDEIN  (input) DOUBLE PRECISION
*          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
*          condition number for the average of selected eigenvalues.
*
*  RCDVIN  (input) DOUBLE PRECISION
*          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
*          condition number for the selected right invariant subspace.
*
*  NSLCT   (input) INTEGER
*          When COMP = .TRUE. the number of selected eigenvalues
*          corresponding to the precomputed values RCDEIN and RCDVIN.
*
*  ISLCT   (input) INTEGER array, dimension (NSLCT)
*          When COMP = .TRUE. ISLCT selects the eigenvalues of the
*          input matrix corresponding to the precomputed values RCDEIN
*          and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
*          eigenvalue with the J-th largest real part is selected.
*          Not referenced if COMP = .FALSE.
*
*  RESULT  (output) DOUBLE PRECISION array, dimension (17)
*          The values computed by the 17 tests described above.
*          The values are currently limited to 1/ulp, to avoid
*          overflow.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The number of entries in WORK to be passed to DGEESX. This
*          must be at least 3*N, and N+N**2 if tests 14--16 are to
*          be performed.
*
*  IWORK   (workspace) INTEGER array, dimension (N*N)
*
*  BWORK   (workspace) LOGICAL array, dimension (N)
*
*  INFO    (output) INTEGER
*          If 0,  successful exit.
*          If <0, input parameter -INFO had an incorrect value.
*          If >0, DGEESX returned an error code, the absolute
*                 value of which is returned.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
      DOUBLE PRECISION   EPSIN
      PARAMETER          ( EPSIN = 5.9605D-8 )
*     ..
*     .. Local Scalars ..
      CHARACTER          SORT
      INTEGER            I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
     $                   RSUB, SDIM, SDIM1
      DOUBLE PRECISION   ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
     $                   SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN,
     $                   VRMIN, WNORM
*     ..
*     .. Local Arrays ..
      INTEGER            IPNT( 20 )
*     ..
*     .. Arrays in Common ..
      LOGICAL            SELVAL( 20 )
      DOUBLE PRECISION   SELWI( 20 ), SELWR( 20 )
*     ..
*     .. Scalars in Common ..
      INTEGER            SELDIM, SELOPT
*     ..
*     .. Common blocks ..
      COMMON             / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
*     ..
*     .. External Functions ..
      LOGICAL            DSLECT
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           DSLECT, DLAMCH, DLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           DCOPY, DGEESX, DGEMM, DLACPY, DORT01, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, MAX, MIN, SIGN, SQRT
*     ..
*     .. Executable Statements ..
*
*     Check for errors
*
      INFO = 0
      IF( THRESH.LT.ZERO ) THEN
         INFO = -3
      ELSE IF( NOUNIT.LE.0 ) THEN
         INFO = -5
      ELSE IF( N.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
         INFO = -8
      ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN
         INFO = -18
      ELSE IF( LWORK.LT.3*N ) THEN
         INFO = -26
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGET24', -INFO )
         RETURN
      END IF
*
*     Quick return if nothing to do
*
      DO 10 I = 1, 17
         RESULT( I ) = -ONE
   10 CONTINUE
*
      IF( N.EQ.0 )
     $   RETURN
*
*     Important constants
*
      SMLNUM = DLAMCH( 'Safe minimum' )
      ULP = DLAMCH( 'Precision' )
      ULPINV = ONE / ULP
*
*     Perform tests (1)-(13)
*
      SELOPT = 0
      LIWORK = N*N
      DO 120 ISORT = 0, 1
         IF( ISORT.EQ.0 ) THEN
            SORT = 'N'
            RSUB = 0
         ELSE
            SORT = 'S'
            RSUB = 6
         END IF
*
*        Compute Schur form and Schur vectors, and test them
*
         CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
         CALL DGEESX( 'V', SORT, DSLECT, 'N', N, H, LDA, SDIM, WR, WI,
     $                VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK,
     $                LIWORK, BWORK, IINFO )
         IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
            RESULT( 1+RSUB ) = ULPINV
            IF( JTYPE.NE.22 ) THEN
               WRITE( NOUNIT, FMT = 9998 )'DGEESX1', IINFO, N, JTYPE,
     $            ISEED
            ELSE
               WRITE( NOUNIT, FMT = 9999 )'DGEESX1', IINFO, N,
     $            ISEED( 1 )
            END IF
            INFO = ABS( IINFO )
            RETURN
         END IF
         IF( ISORT.EQ.0 ) THEN
            CALL DCOPY( N, WR, 1, WRTMP, 1 )
            CALL DCOPY( N, WI, 1, WITMP, 1 )
         END IF
*
*        Do Test (1) or Test (7)
*
         RESULT( 1+RSUB ) = ZERO
         DO 30 J = 1, N - 2
            DO 20 I = J + 2, N
               IF( H( I, J ).NE.ZERO )
     $            RESULT( 1+RSUB ) = ULPINV
   20       CONTINUE
   30    CONTINUE
         DO 40 I = 1, N - 2
            IF( H( I+1, I ).NE.ZERO .AND. H( I+2, I+1 ).NE.ZERO )
     $         RESULT( 1+RSUB ) = ULPINV
   40    CONTINUE
         DO 50 I = 1, N - 1
            IF( H( I+1, I ).NE.ZERO ) THEN
               IF( H( I, I ).NE.H( I+1, I+1 ) .OR. H( I, I+1 ).EQ.
     $             ZERO .OR. SIGN( ONE, H( I+1, I ) ).EQ.
     $             SIGN( ONE, H( I, I+1 ) ) )RESULT( 1+RSUB ) = ULPINV
            END IF
   50    CONTINUE
*
*        Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
*
*        Copy A to VS1, used as workspace
*
         CALL DLACPY( ' ', N, N, A, LDA, VS1, LDVS )
*
*        Compute Q*H and store in HT.
*
         CALL DGEMM( 'No transpose', 'No transpose', N, N, N, ONE, VS,
     $               LDVS, H, LDA, ZERO, HT, LDA )
*
*        Compute A - Q*H*Q'
*
         CALL DGEMM( 'No transpose', 'Transpose', N, N, N, -ONE, HT,
     $               LDA, VS, LDVS, ONE, VS1, LDVS )
*
         ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), SMLNUM )
         WNORM = DLANGE( '1', N, N, VS1, LDVS, WORK )
*
         IF( ANORM.GT.WNORM ) THEN
            RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP )
         ELSE
            IF( ANORM.LT.ONE ) THEN
               RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) /
     $                            ( N*ULP )
            ELSE
               RESULT( 2+RSUB ) = MIN( WNORM / ANORM, DBLE( N ) ) /
     $                            ( N*ULP )
            END IF
         END IF
*
*        Test (3) or (9):  Compute norm( I - Q'*Q ) / ( N * ULP )
*
         CALL DORT01( 'Columns', N, N, VS, LDVS, WORK, LWORK,
     $                RESULT( 3+RSUB ) )
*
*        Do Test (4) or Test (10)
*
         RESULT( 4+RSUB ) = ZERO
         DO 60 I = 1, N
            IF( H( I, I ).NE.WR( I ) )
     $         RESULT( 4+RSUB ) = ULPINV
   60    CONTINUE
         IF( N.GT.1 ) THEN
            IF( H( 2, 1 ).EQ.ZERO .AND. WI( 1 ).NE.ZERO )
     $         RESULT( 4+RSUB ) = ULPINV
            IF( H( N, N-1 ).EQ.ZERO .AND. WI( N ).NE.ZERO )
     $         RESULT( 4+RSUB ) = ULPINV
         END IF
         DO 70 I = 1, N - 1
            IF( H( I+1, I ).NE.ZERO ) THEN
               TMP = SQRT( ABS( H( I+1, I ) ) )*
     $               SQRT( ABS( H( I, I+1 ) ) )
               RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
     $                            ABS( WI( I )-TMP ) /
     $                            MAX( ULP*TMP, SMLNUM ) )
               RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
     $                            ABS( WI( I+1 )+TMP ) /
     $                            MAX( ULP*TMP, SMLNUM ) )
            ELSE IF( I.GT.1 ) THEN
               IF( H( I+1, I ).EQ.ZERO .AND. H( I, I-1 ).EQ.ZERO .AND.
     $             WI( I ).NE.ZERO )RESULT( 4+RSUB ) = ULPINV
            END IF
   70    CONTINUE
*
*        Do Test (5) or Test (11)
*
         CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )