ddrges.f

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

      SUBROUTINE DDRGES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
     $                   NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR,
     $                   ALPHAI, BETA, WORK, LWORK, RESULT, 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, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * ), DOTYPE( * )
      INTEGER            ISEED( 4 ), NN( * )
      DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
     $                   B( LDA, * ), BETA( * ), Q( LDQ, * ),
     $                   RESULT( 13 ), S( LDA, * ), T( LDA, * ),
     $                   WORK( * ), Z( LDQ, * )
*     ..
*
*  Purpose
*  =======
*
*  DDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
*  problem driver DGGES.
*
*  DGGES 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(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 DDRGES 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. quasi-triangular form)
*        (no sorting of eigenvalues)
*
*  (6)   if eigenvalues = diagonal blocks of the Schur form (S, T),
*        i.e., test the 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 and j+1-th eigenvalues.
*        (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 blocks of the Schur form (S, T),
*        i.e. test the 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 and j+1-th eigenvalues.
*        (with sorting of eigenvalues).
*
*  (12)  if sorting worked and SDIM is the number of eigenvalues
*        which were SELECTed.
*
*  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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDA, max(NN))
*          The Schur form matrix computed from A by DGGES.  On exit, S
*          contains the Schur form matrix corresponding to the matrix
*          in A.
*
*  T       (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
*          The upper triangular matrix computed from B by DGGES.
*
*  Q       (workspace) DOUBLE PRECISION array, dimension (LDQ, max(NN))
*          The (left) orthogonal matrix computed by DGGES.
*
*  LDQ     (input) INTEGER
*          The leading dimension of Q and Z. It must
*          be at least 1 and at least max( NN ).
*
*  Z       (workspace) DOUBLE PRECISION array, dimension( LDQ, max(NN) )
*          The (right) orthogonal matrix computed by DGGES.
*
*  ALPHAR  (workspace) DOUBLE PRECISION array, dimension (max(NN))
*  ALPHAI  (workspace) DOUBLE PRECISION array, dimension (max(NN))
*  BETA    (workspace) DOUBLE PRECISION array, dimension (max(NN))
*          The generalized eigenvalues of (A,B) computed by DGGES.
*          ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
*          generalized eigenvalue of A and B.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.
*          LWORK >= MAX( 10*(N+1), 3*N*N ), where N is the largest
*          matrix dimension.
*
*  RESULT  (output) DOUBLE PRECISION array, dimension (15)