zdrvgg.f

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

      SUBROUTINE ZDRVGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
     $                   THRSHN, NOUNIT, A, LDA, B, S, T, S2, T2, Q,
     $                   LDQ, Z, ALPHA1, BETA1, ALPHA2, BETA2, VL, VR,
     $                   WORK, LWORK, RWORK, RESULT, INFO )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
      DOUBLE PRECISION   THRESH, THRSHN
*     ..
*     .. Array Arguments ..
*
*  Purpose
*  =======
*
*  ZDRVGG  checks the nonsymmetric generalized eigenvalue driver
*  routines.
*                                T          T        T
*  ZGEGS factors A and B as Q S Z  and Q T Z , where   means
*  transpose, T is upper triangular, S is in generalized Schur form
*  (upper triangular), and Q and Z are unitary.  It also
*  computes the generalized eigenvalues (alpha(1),beta(1)), ...,
*  (alpha(n),beta(n)), where alpha(j)=S(j,j) and beta(j)=T(j,j) --
*  thus, w(j) = alpha(j)/beta(j) is a root of the generalized
*  eigenvalue problem
*
*      det( A - w(j) B ) = 0
*
*  and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
*  problem
*
*      det( m(j) A - B ) = 0
*
*  ZGEGV computes the generalized eigenvalues (alpha(1),beta(1)), ...,
*  (alpha(n),beta(n)), the matrix L whose columns contain the
*  generalized left eigenvectors l, and the matrix R whose columns
*  contain the generalized right eigenvectors r for the pair (A,B).
*
*  When ZDRVGG 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, 7
*  tests will be performed and compared with the threshhold THRESH:
*
*  Results from ZGEGS:
*
*                   H
*  (1)   | A - Q S Z  | / ( |A| n ulp )
*
*                   H
*  (2)   | B - Q T Z  | / ( |B| n ulp )
*
*                H
*  (3)   | I - QQ  | / ( n ulp )
*
*                H
*  (4)   | I - ZZ  | / ( n ulp )
*
*  (5)   maximum over j of D(j)  where:
*
*                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
*            D(j) = ------------------------ + -----------------------
*                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
*
*  Results from ZGEGV:
*
*  (6)   max over all left eigenvalue/-vector pairs (beta/alpha,l) of
*
*     | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
*
*        where l**H is the conjugate tranpose of l.
*
*  (7)   max over all right eigenvalue/-vector pairs (beta/alpha,r) of
*
*        | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
*
*  Test Matrices
*  ---- --------
*
*  The sizes of the test matrices 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)  ( 0, 0 )         (a pair of zero matrices)
*
*  (2)  ( I, 0 )         (an identity and a zero matrix)
*
*  (3)  ( 0, I )         (an identity and a zero matrix)
*
*  (4)  ( I, I )         (a pair of identity matrices)
*
*          t   t
*  (5)  ( J , J  )       (a pair of transposed Jordan blocks)
*
*                                      t                ( I   0  )
*  (6)  ( X, Y )         where  X = ( J   0  )  and Y = (      t )
*                                   ( 0   I  )          ( 0   J  )
*                        and I is a k x k identity and J a (k+1)x(k+1)
*                        Jordan block; k=(N-1)/2
*
*  (7)  ( D, I )         where D is diag( 0, 1,..., N-1 ) (a diagonal
*                        matrix with those diagonal entries.)
*  (8)  ( I, D )
*
*  (9)  ( big*D, small*I ) where "big" is near overflow and small=1/big
*
*  (10) ( small*D, big*I )
*
*  (11) ( big*I, small*D )
*
*  (12) ( small*I, big*D )
*
*  (13) ( big*D, big*I )
*
*  (14) ( small*D, small*I )
*
*  (15) ( D1, D2 )        where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
*                         D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
*            t   t
*  (16) Q ( J , J ) Z     where Q and Z are random unitary matrices.
*
*  (17) Q ( T1, T2 ) Z    where T1 and T2 are upper triangular matrices
*                         with random O(1) entries above the diagonal
*                         and diagonal entries diag(T1) =
*                         ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
*                         ( 0, N-3, N-4,..., 1, 0, 0 )
*
*  (18) Q ( T1, T2 ) Z    diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
*                         diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
*                         s = machine precision.
*
*  (19) Q ( T1, T2 ) Z    diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
*                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
*
*                                                         N-5
*  (20) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, 1, a, ..., a   =s, 0 )
*                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
*
*  (21) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
*                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
*                         where r1,..., r(N-4) are random.
*
*  (22) Q ( big*T1, small*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
*
*  (23) Q ( small*T1, big*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
*
*  (24) Q ( small*T1, small*T2 ) Z  diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
*
*  (25) Q ( big*T1, big*T2 ) Z      diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
*
*  (26) Q ( T1, T2 ) Z     where T1 and T2 are random upper-triangular
*                          matrices.
*
*  Arguments
*  =========
*
*  NSIZES  (input) INTEGER
*          The number of sizes of matrices to use.  If it is zero,
*          ZDRVGG does nothing.  It must be at least zero.
*
*  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.   If it is zero, ZDRVGG
*          does nothing.  It must be at least zero.  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 ZDRVGG 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.
*
*  THRSHN  (input) DOUBLE PRECISION
*          Threshhold for reporting eigenvector normalization error.
*          If the normalization of any eigenvector differs from 1 by
*          more than THRSHN*ulp, then a special error message will be
*          printed.  (This is handled separately from the other tests,
*          since only a compiler or programming error should cause an
*          error message, at least if THRSHN is at least 5--10.)
*
*  NOUNIT  (input) INTEGER
*          The FORTRAN unit number for printing out error messages
*          (e.g., if a routine returns IINFO not equal to 0.)
*
*  A       (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN))
*          Used to hold the original A matrix.  Used as input only
*          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
*          DOTYPE(MAXTYP+1)=.TRUE.
*
*  LDA     (input) INTEGER
*          The leading dimension of A, B, S, T, S2, and T2.
*          It must be at least 1 and at least max( NN ).
*
*  B       (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN))
*          Used to hold the original B matrix.  Used as input only
*          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
*          DOTYPE(MAXTYP+1)=.TRUE.
*
*  S       (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
*          The upper triangular matrix computed from A by ZGEGS.
*
*  T       (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
*          The upper triangular matrix computed from B by ZGEGS.
*
*  S2      (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
*          The matrix computed from A by ZGEGV.  This will be the
*          Schur (upper triangular) form of some matrix related to A,
*          but will not, in general, be the same as S.
*
*  T2      (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
*          The matrix computed from B by ZGEGV.  This will be the
*          Schur form of some matrix related to B, but will not, in
*          general, be the same as T.
*
*  Q       (workspace) COMPLEX*16 array, dimension (LDQ, max(NN))
*          The (left) unitary matrix computed by ZGEGS.
*
*  LDQ     (input) INTEGER
*          The leading dimension of Q, Z, VL, and VR.  It must
*          be at least 1 and at least max( NN ).
*
*  Z       (workspace) COMPLEX*16 array, dimension (LDQ, max(NN))
*          The (right) unitary matrix computed by ZGEGS.
*
*  ALPHA1  (workspace) COMPLEX*16 array, dimension (max(NN))
*  BETA1   (workspace) COMPLEX*16 array, dimension (max(NN))
*          The generalized eigenvalues of (A,B) computed by ZGEGS.
*          ALPHA1(k) / BETA1(k)  is the k-th generalized eigenvalue of
*          the matrices in A and B.
*
*  ALPHA2  (workspace) COMPLEX*16 array, dimension (max(NN))
*  BETA2   (workspace) COMPLEX*16 array, dimension (max(NN))
*          The generalized eigenvalues of (A,B) computed by ZGEGV.
*          ALPHA2(k) / BETA2(k)  is the k-th generalized eigenvalue of