zdrges.f

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

      SUBROUTINE ZDRGES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
     $                   NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHA,
     $                   BETA, WORK, LWORK, RWORK, RESULT, BWORK, INFO )
*
*  -- LAPACK test routine (version 3.1.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     February 2007
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * ), DOTYPE( * )
      INTEGER            ISEED( 4 ), NN( * )
      DOUBLE PRECISION   RESULT( 13 ), RWORK( * )
      COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDA, * ),
     $                   BETA( * ), Q( LDQ, * ), S( LDA, * ),
     $                   T( LDA, * ), WORK( * ), Z( LDQ, * )
*     ..
*
*  Purpose
*  =======
*
*  ZDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
*  problem driver ZGGES.
*
*  ZGGES factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate
*  transpose, S and T are  upper triangular (i.e., in generalized Schur
*  form), and Q and Z are unitary. It also computes the generalized
*  eigenvalues (alpha(j),beta(j)), j=1,...,n.  Thus,
*  w(j) = alpha(j)/beta(j) is a root of the characteristic equation
*
*                  det( A - w(j) B ) = 0
*
*  Optionally it also reorder the eigenvalues so that a selected
*  cluster of eigenvalues appears in the leading diagonal block of the
*  Schur forms.
*
*  When ZDRGES 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, a pair of matrices (A, B) will be generated
*  and used for testing. For each matrix pair, the following 13 tests
*  will be performed and compared with the threshhold THRESH except
*  the tests (5), (11) and (13).
*
*
*  (1)   | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)
*
*
*  (2)   | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)
*
*
*  (3)   | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)
*
*
*  (4)   | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)
*
*  (5)   if A is in Schur form (i.e. triangular form) (no sorting of
*        eigenvalues)
*
*  (6)   if eigenvalues = diagonal elements of the Schur form (S, T),
*        i.e., test the 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)|)
*
*        (no sorting of eigenvalues)
*
*  (7)   | (A,B) - Q (S,T) Z' | / ( |(A,B)| n ulp )
*        (with sorting of eigenvalues).
*
*  (8)   | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).
*
*  (9)   | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).
*
*  (10)  if A is in Schur form (i.e. quasi-triangular form)
*        (with sorting of eigenvalues).
*
*  (11)  if eigenvalues = diagonal elements of the Schur form (S, T),
*        i.e. test the 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)|)
*
*        (with sorting of eigenvalues).
*
*  (12)  if sorting worked and SDIM is the number of eigenvalues
*        which were CELECTed.
*
*  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 orthogonal 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,
*          DDRGES does nothing.  NSIZES >= 0.
*
*  NN      (input) INTEGER array, dimension (NSIZES)
*          An array containing the sizes to be used for the matrices.
*          Zero values will be skipped.  NN >= 0.
*
*  NTYPES  (input) INTEGER
*          The number of elements in DOTYPE.   If it is zero, DDRGES
*          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 on input.
*          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 DDRGES 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.  THRESH >= 0.
*
*  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, and T.
*          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 Schur form matrix computed from A by ZGGES.  On exit, S
*          contains the Schur form matrix corresponding to the matrix
*          in A.
*
*  T       (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
*          The upper triangular matrix computed from B by ZGGES.
*
*  Q       (workspace) COMPLEX*16 array, dimension (LDQ, max(NN))
*          The (left) orthogonal matrix computed by ZGGES.
*
*  LDQ     (input) INTEGER
*          The leading dimension of Q and Z. It must
*          be at least 1 and at least max( NN ).
*
*  Z       (workspace) COMPLEX*16 array, dimension( LDQ, max(NN) )
*          The (right) orthogonal matrix computed by ZGGES.
*
*  ALPHA   (workspace) COMPLEX*16 array, dimension (max(NN))
*  BETA    (workspace) COMPLEX*16 array, dimension (max(NN))
*          The generalized eigenvalues of (A,B) computed by ZGGES.
*          ALPHA(k) / BETA(k) is the k-th generalized eigenvalue of A
*          and B.
*
*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.  LWORK >= 3*N*N.
*
*  RWORK   (workspace) DOUBLE PRECISION array, dimension ( 8*N )
*          Real workspace.
*
*  RESULT  (output) DOUBLE PRECISION array, dimension (15)
*          The values computed by the tests described above.