ddrvsx.f

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

      SUBROUTINE DDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
     $                   NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT,
     $                   WIT, WRTMP, WITMP, VS, LDVS, VS1, 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 ..
      INTEGER            INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
     $                   NTYPES
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * ), DOTYPE( * )
      INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
      DOUBLE PRECISION   A( LDA, * ), H( LDA, * ), HT( LDA, * ),
     $                   RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
     $                   WI( * ), WIT( * ), WITMP( * ), WORK( * ),
     $                   WR( * ), WRT( * ), WRTMP( * )
*     ..
*
*  Purpose
*  =======
*
*     DDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
*     expert driver DGEESX.
*
*     DDRVSX uses both test matrices generated randomly depending on
*     data supplied in the calling sequence, as well as on data
*     read from an input file and including precomputed condition
*     numbers to which it compares the ones it computes.
*
*     When DDRVSX is called, a number of matrix "sizes" ("n's") and a
*     number of matrix "types" are specified.  For each size ("n")
*     and each type of matrix, one matrix will be generated and used
*     to test the nonsymmetric eigenroutines.  For each matrix, 15
*     tests 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
*
*     The "sizes" are specified by an array NN(1:NSIZES); the value of
*     each element NN(j) specifies one size.
*     The "types" are specified by a logical array DOTYPE( 1:NTYPES );
*     if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
*     Currently, the list of possible types is:
*
*     (1)  The zero matrix.
*     (2)  The identity matrix.
*     (3)  A (transposed) Jordan block, with 1's on the diagonal.
*
*     (4)  A diagonal matrix with evenly spaced entries
*          1, ..., ULP  and random signs.
*          (ULP = (first number larger than 1) - 1 )
*     (5)  A diagonal matrix with geometrically spaced entries
*          1, ..., ULP  and random signs.
*     (6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
*          and random signs.
*
*     (7)  Same as (4), but multiplied by a constant near
*          the overflow threshold
*     (8)  Same as (4), but multiplied by a constant near
*          the underflow threshold
*
*     (9)  A matrix of the form  U' T U, where U is orthogonal and
*          T has evenly spaced entries 1, ..., ULP with random signs
*          on the diagonal and random O(1) entries in the upper
*          triangle.
*
*     (10) A matrix of the form  U' T U, where U is orthogonal and
*          T has geometrically spaced entries 1, ..., ULP with random
*          signs on the diagonal and random O(1) entries in the upper
*          triangle.
*
*     (11) A matrix of the form  U' T U, where U is orthogonal and
*          T has "clustered" entries 1, ULP,..., ULP with random
*          signs on the diagonal and random O(1) entries in the upper
*          triangle.
*
*     (12) A matrix of the form  U' T U, where U is orthogonal and
*          T has real or complex conjugate paired eigenvalues randomly
*          chosen from ( ULP, 1 ) and random O(1) entries in the upper
*          triangle.
*
*     (13) A matrix of the form  X' T X, where X has condition
*          SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
*          with random signs on the diagonal and random O(1) entries
*          in the upper triangle.
*
*     (14) A matrix of the form  X' T X, where X has condition
*          SQRT( ULP ) and T has geometrically spaced entries
*          1, ..., ULP with random signs on the diagonal and random
*          O(1) entries in the upper triangle.
*
*     (15) A matrix of the form  X' T X, where X has condition
*          SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
*          with random signs on the diagonal and random O(1) entries
*          in the upper triangle.
*
*     (16) A matrix of the form  X' T X, where X has condition
*          SQRT( ULP ) and T has real or complex conjugate paired
*          eigenvalues randomly chosen from ( ULP, 1 ) and random
*          O(1) entries in the upper triangle.
*
*     (17) Same as (16), but multiplied by a constant
*          near the overflow threshold
*     (18) Same as (16), but multiplied by a constant
*          near the underflow threshold
*
*     (19) Nonsymmetric matrix with random entries chosen from (-1,1).
*          If N is at least 4, all entries in first two rows and last
*          row, and first column and last two columns are zero.
*     (20) Same as (19), but multiplied by a constant
*          near the overflow threshold
*     (21) Same as (19), but multiplied by a constant
*          near the underflow threshold
*
*     In addition, an input file will be read from logical unit number
*     NIUNIT. The file contains matrices along with precomputed
*     eigenvalues and reciprocal condition numbers for the eigenvalue
*     average and right invariant subspace. For these matrices, in
*     addition to tests (1) to (15) we will compute the following two
*     tests:
*
*    (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
*  =========
*
*  NSIZES  (input) INTEGER
*          The number of sizes of matrices to use.  NSIZES must be at
*          least zero. If it is zero, no randomly generated matrices
*          are tested, but any test matrices read from NIUNIT will be
*          tested.
*
*  NN      (input) INTEGER array, dimension (NSIZES)
*          An array containing the sizes to be used for the matrices.
*          Zero values will be skipped.  The values must be at least
*          zero.
*
*  NTYPES  (input) INTEGER
*          The number of elements in DOTYPE. NTYPES must be at least
*          zero. If it is zero, no randomly generated test matrices
*          are tested, but and test matrices read from NIUNIT will be
*          tested. If it is MAXTYP+1 and NSIZES is 1, then an
*          additional type, MAXTYP+1 is defined, which is to use
*          whatever matrix is in A.  This is only useful if
*          DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
*
*  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
*          If DOTYPE(j) is .TRUE., then for each size in NN a
*          matrix of that size and of type j will be generated.
*          If NTYPES is smaller than the maximum number of types
*          defined (PARAMETER MAXTYP), then types NTYPES+1 through
*          MAXTYP will not be generated.  If NTYPES is larger
*          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
*          will be ignored.
*
*  ISEED   (input/output) INTEGER array, dimension (4)
*          On entry ISEED specifies the seed of the random number
*          generator. The array elements should be between 0 and 4095;
*          if not they will be reduced mod 4096.  Also, ISEED(4) must
*          be odd.  The random number generator uses a linear
*          congruential sequence limited to small integers, and so
*          should produce machine independent random numbers. The
*          values of ISEED are changed on exit, and can be used in the
*          next call to DDRVSX to continue the same random number
*          sequence.
*
*  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.
*
*  NIUNIT  (input) INTEGER
*          The FORTRAN unit number for reading in the data file of
*          problems to solve.
*
*  NOUNIT  (input) INTEGER
*          The FORTRAN unit number for printing out error messages
*          (e.g., if a routine returns INFO not equal to 0.)
*
*  A       (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
*          Used to hold the matrix whose eigenvalues are to be
*          computed.  On exit, A contains the last matrix actually used.
*
*  LDA     (input) INTEGER
*          The leading dimension of A, and H. LDA must be at
*          least 1 and at least max( NN ).
*
*  H       (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
*          Another copy of the test matrix A, modified by DGEESX.
*
*  HT      (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
*          Yet another copy of the test matrix A, modified by DGEESX.
*