sdrvgg.f

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

      SUBROUTINE SDRVGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
     $                   THRSHN, NOUNIT, A, LDA, B, S, T, S2, T2, Q,
     $                   LDQ, Z, ALPHR1, ALPHI1, BETA1, ALPHR2, ALPHI2,
     $                   BETA2, VL, VR, WORK, LWORK, 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
      REAL               THRESH, THRSHN
*     ..
*     .. Array Arguments ..
      LOGICAL            DOTYPE( * )
      INTEGER            ISEED( 4 ), NN( * )
      REAL               A( LDA, * ), ALPHI1( * ), ALPHI2( * ),
     $                   ALPHR1( * ), ALPHR2( * ), B( LDA, * ),
     $                   BETA1( * ), BETA2( * ), Q( LDQ, * ),
     $                   RESULT( * ), S( LDA, * ), S2( LDA, * ),
     $                   T( LDA, * ), T2( LDA, * ), VL( LDQ, * ),
     $                   VR( LDQ, * ), WORK( * ), Z( LDQ, * )
*     ..
*
*  Purpose
*  =======
*
*  SDRVGG  checks the nonsymmetric generalized eigenvalue driver
*  routines.
*                                T          T        T
*  SGEGS 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
*  (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
*  the 2x2 blocks corresponding to complex conjugate pairs of
*  generalized eigenvalues), and Q and Z are orthogonal.  It also
*  computes the generalized eigenvalues (alpha(1),beta(1)), ...,
*  (alpha(n),beta(n)), where alpha(j)=S(j,j) and beta(j)=P(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
*
*  SGEGV 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 SDRVGG 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 SGEGS:
*
*                   T
*  (1)   | A - Q S Z  | / ( |A| n ulp )
*
*                   T
*  (2)   | B - Q T Z  | / ( |B| n ulp )
*
*                T
*  (3)   | I - QQ  | / ( n ulp )
*
*                T
*  (4)   | I - ZZ  | / ( n ulp )
*
*  (5)   maximum over j of D(j)  where:
*
*  if alpha(j) is real:
*                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
*            D(j) = ------------------------ + -----------------------
*                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
*
*  if alpha(j) is complex:
*                                  | det( s S - w T ) |
*            D(j) = ---------------------------------------------------
*                   ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
*
*            and S and T are here the 2 x 2 diagonal blocks of S and T
*            corresponding to the j-th eigenvalue.
*
*  Results from SGEGV:
*
*  (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 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,
*          SDRVGG 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, SDRVGG
*          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 SDRVGG to continue the same random number
*          sequence.
*
*  THRESH  (input) REAL
*          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) REAL
*          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) REAL 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) REAL 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) REAL array, dimension (LDA, max(NN))
*          The Schur form matrix computed from A by SGEGS.  On exit, S
*          contains the Schur form matrix corresponding to the matrix
*          in A.
*
*  T       (workspace) REAL array, dimension (LDA, max(NN))
*          The upper triangular matrix computed from B by SGEGS.
*
*  S2      (workspace) REAL array, dimension (LDA, max(NN))
*          The matrix computed from A by SGEGV.  This will be the
*          Schur form of some matrix related to A, but will not, in
*          general, be the same as S.
*
*  T2      (workspace) REAL array, dimension (LDA, max(NN))
*          The matrix computed from B by SGEGV.  This will be the
*          Schur form of some matrix related to B, but will not, in
*          general, be the same as T.
*
*  Q       (workspace) REAL array, dimension (LDQ, max(NN))
*          The (left) orthogonal matrix computed by SGEGS.
*
*  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) REAL array of
*                             dimension( LDQ, max(NN) )
*          The (right) orthogonal matrix computed by SGEGS.
*
*  ALPHR1  (workspace) REAL array, dimension (max(NN))
*  ALPHI1  (workspace) REAL array, dimension (max(NN))
*  BETA1   (workspace) REAL array, dimension (max(NN))
*
*          The generalized eigenvalues of (A,B) computed by SGEGS.
*          ( ALPHR1(k)+ALPHI1(k)*i ) / BETA1(k) is the k-th
*          generalized eigenvalue of the matrices in A and B.
*
*  ALPHR2  (workspace) REAL array, dimension (max(NN))
*  ALPHI2  (workspace) REAL array, dimension (max(NN))