ddrgsx.f

famous linear algebra library (LAPACK) ports to windows

      SUBROUTINE DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
     $                   BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
     $                   WORK, LWORK, IWORK, LIWORK, 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, LDC, LIWORK, LWORK, NCMAX, NIN,
     $                   NOUT, NSIZE
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * )
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
     $                   ALPHAR( * ), B( LDA, * ), BETA( * ),
     $                   BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
     $                   WORK( * ), Z( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
*  problem expert driver DGGESX.
*
*  DGGESX 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)). Thus, w(j) = alpha(j)/beta(j) is a root of the
*  characteristic equation
*
*      det( A - w(j) B ) = 0
*
*  Optionally it also reorders the eigenvalues so that a selected
*  cluster of eigenvalues appears in the leading diagonal block of the
*  Schur forms; computes a reciprocal condition number for the average
*  of the selected eigenvalues; and computes a reciprocal condition
*  number for the right and left deflating subspaces corresponding to
*  the selected eigenvalues.
*
*  When DDRGSX is called with NSIZE > 0, five (5) types of built-in
*  matrix pairs are used to test the routine DGGESX.
*
*  When DDRGSX is called with NSIZE = 0, it reads in test matrix data
*  to test DGGESX.
*
*  For each matrix pair, the following tests will be performed and
*  compared with the threshhold THRESH except for the tests (7) and (9):
*
*  (1)   | A - Q S Z' | / ( |A| n ulp )
*
*  (2)   | B - Q T Z' | / ( |B| n ulp )
*
*  (3)   | I - QQ' | / ( n ulp )
*
*  (4)   | I - ZZ' | / ( n ulp )
*
*  (5)   if A is in Schur form (i.e. quasi-triangular form)
*
*  (6)   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.
*
*  (7)   if sorting worked and SDIM is the number of eigenvalues
*        which were selected.
*
*  (8)   the estimated value DIF does not differ from the true values of
*        Difu and Difl more than a factor 10*THRESH. If the estimate DIF
*        equals zero the corresponding true values of Difu and Difl
*        should be less than EPS*norm(A, B). If the true value of Difu
*        and Difl equal zero, the estimate DIF should be less than
*        EPS*norm(A, B).
*
*  (9)   If INFO = N+3 is returned by DGGESX, the reordering "failed"
*        and we check that DIF = PL = PR = 0 and that the true value of
*        Difu and Difl is < EPS*norm(A, B). We count the events when
*        INFO=N+3.
*
*  For read-in test matrices, the above tests are run except that the
*  exact value for DIF (and PL) is input data.  Additionally, there is
*  one more test run for read-in test matrices:
*
*  (10)  the estimated value PL does not differ from the true value of
*        PLTRU more than a factor THRESH. If the estimate PL equals
*        zero the corresponding true value of PLTRU should be less than
*        EPS*norm(A, B). If the true value of PLTRU equal zero, the
*        estimate PL should be less than EPS*norm(A, B).
*
*  Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
*  matrix pairs are generated and tested. NSIZE should be kept small.
*
*  SVD (routine DGESVD) is used for computing the true value of DIF_u
*  and DIF_l when testing the built-in test problems.
*
*  Built-in Test Matrices
*  ======================
*
*  All built-in test matrices are the 2 by 2 block of triangular
*  matrices
*
*           A = [ A11 A12 ]    and      B = [ B11 B12 ]
*               [     A22 ]                 [     B22 ]
*
*  where for different type of A11 and A22 are given as the following.
*  A12 and B12 are chosen so that the generalized Sylvester equation
*
*           A11*R - L*A22 = -A12
*           B11*R - L*B22 = -B12
*
*  have prescribed solution R and L.
*
*  Type 1:  A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
*           B11 = I_m, B22 = I_k
*           where J_k(a,b) is the k-by-k Jordan block with ``a'' on
*           diagonal and ``b'' on superdiagonal.
*
*  Type 2:  A11 = (a_ij) = ( 2(.5-sin(i)) ) and
*           B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
*           A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
*           B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
*
*  Type 3:  A11, A22 and B11, B22 are chosen as for Type 2, but each
*           second diagonal block in A_11 and each third diagonal block
*           in A_22 are made as 2 by 2 blocks.
*
*  Type 4:  A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
*              for i=1,...,m,  j=1,...,m and
*           A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
*              for i=m+1,...,k,  j=m+1,...,k
*
*  Type 5:  (A,B) and have potentially close or common eigenvalues and
*           very large departure from block diagonality A_11 is chosen
*           as the m x m leading submatrix of A_1:
*                   |  1  b                            |
*                   | -b  1                            |
*                   |        1+d  b                    |
*                   |         -b 1+d                   |
*            A_1 =  |                  d  1            |
*                   |                 -1  d            |
*                   |                        -d  1     |
*                   |                        -1 -d     |
*                   |                               1  |
*           and A_22 is chosen as the k x k leading submatrix of A_2:
*                   | -1  b                            |
*                   | -b -1                            |
*                   |       1-d  b                     |
*                   |       -b  1-d                    |
*            A_2 =  |                 d 1+b            |
*                   |               -1-b d             |
*                   |                       -d  1+b    |
*                   |                      -1+b  -d    |
*                   |                              1-d |
*           and matrix B are chosen as identity matrices (see DLATM5).
*
*
*  Arguments
*  =========
*
*  NSIZE   (input) INTEGER
*          The maximum size of the matrices to use. NSIZE >= 0.
*          If NSIZE = 0, no built-in tests matrices are used, but
*          read-in test matrices are used to test DGGESX.
*
*  NCMAX   (input) INTEGER
*          Maximum allowable NMAX for generating Kroneker matrix
*          in call to DLAKF2
*
*  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.
*
*  NIN     (input) INTEGER
*          The FORTRAN unit number for reading in the data file of
*          problems to solve.
*
*  NOUT    (input) INTEGER
*          The FORTRAN unit number for printing out error messages
*          (e.g., if a routine returns IINFO not equal to 0.)
*
*  A       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
*          Used to store 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, B, AI, BI, Z and Q,
*          LDA >= max( 1, NSIZE ). For the read-in test,
*          LDA >= max( 1, N ), N is the size of the test matrices.
*
*  B       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
*          Used to store the matrix whose eigenvalues are to be
*          computed.  On exit, B contains the last matrix actually used.
*
*  AI      (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
*          Copy of A, modified by DGGESX.
*
*  BI      (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
*          Copy of B, modified by DGGESX.
*
*  Z       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
*          Z holds the left Schur vectors computed by DGGESX.
*
*  Q       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
*          Q holds the right Schur vectors computed by DGGESX.
*
*  ALPHAR  (workspace) DOUBLE PRECISION array, dimension (NSIZE)
*  ALPHAI  (workspace) DOUBLE PRECISION array, dimension (NSIZE)
*  BETA    (workspace) DOUBLE PRECISION array, dimension (NSIZE)
*          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
*
*  C       (workspace) DOUBLE PRECISION array, dimension (LDC, LDC)
*          Store the matrix generated by subroutine DLAKF2, this is the
*          matrix formed by Kronecker products used for estimating
*          DIF.
*
*  LDC     (input) INTEGER
*          The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
*
*  S       (workspace) DOUBLE PRECISION array, dimension (LDC)
*          Singular values of C
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.
*          LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
*
*  IWORK   (workspace) INTEGER array, dimension (LIWORK)
*
*  LIWORK  (input) INTEGER
*          The dimension of the array IWORK. LIWORK >= NSIZE + 6.
*
*  BWORK   (workspace) LOGICAL array, dimension (LDA)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*          > 0:  A routine returned an error code.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TEN
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ILABAD
      CHARACTER          SENSE
      INTEGER            BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
     $                   MINWRK, MM, MN2, NERRS, NPTKNT, NTEST, NTESTT,
     $                   PRTYPE, QBA, QBB
      DOUBLE PRECISION   ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
     $                   TEMP2, THRSH2, ULP, ULPINV, WEIGHT
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   DIFEST( 2 ), PL( 2 ), RESULT( 10 )
*     ..
*     .. External Functions ..
      LOGICAL            DLCTSX
      INTEGER            ILAENV
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           DLCTSX, ILAENV, DLAMCH, DLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           ALASVM, DGESVD, DGET51, DGET53, DGGESX, DLABAD,
     $                   DLACPY, DLAKF2, DLASET, DLATM5, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, SQRT
*     ..
*     .. Scalars in Common ..
      LOGICAL            FS
      INTEGER            K, M, MPLUSN, N
*     ..
*     .. Common blocks ..
      COMMON             / MN / M, N, MPLUSN, K, FS
*     ..
*     .. Executable Statements ..
*
*     Check for errors
*
      IF( NSIZE.LT.0 ) THEN