TESTING/EIG/cdrges.f

   1 *> \brief \b CDRGES
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE CDRGES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
  12 *                          NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHA,
  13 *                          BETA, WORK, LWORK, RWORK, RESULT, BWORK, INFO )
  14 *
  15 *       .. Scalar Arguments ..
  16 *       INTEGER            INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
  17 *       REAL               THRESH
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       LOGICAL            BWORK( * ), DOTYPE( * )
  21 *       INTEGER            ISEED( 4 ), NN( * )
  22 *       REAL               RESULT( 13 ), RWORK( * )
  23 *       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDA, * ),
  24 *      $                   BETA( * ), Q( LDQ, * ), S( LDA, * ),
  25 *      $                   T( LDA, * ), WORK( * ), Z( LDQ, * )
  26 *       ..
  27 *
  28 *
  29 *> \par Purpose:
  30 *  =============
  31 *>
  32 *> \verbatim
  33 *>
  34 *> CDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
  35 *> problem driver CGGES.
  36 *>
  37 *> CGGES factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate
  38 *> transpose, S and T are  upper triangular (i.e., in generalized Schur
  39 *> form), and Q and Z are unitary. It also computes the generalized
  40 *> eigenvalues (alpha(j),beta(j)), j=1,...,n.  Thus,
  41 *> w(j) = alpha(j)/beta(j) is a root of the characteristic equation
  42 *>
  43 *>                 det( A - w(j) B ) = 0
  44 *>
  45 *> Optionally it also reorder the eigenvalues so that a selected
  46 *> cluster of eigenvalues appears in the leading diagonal block of the
  47 *> Schur forms.
  48 *>
  49 *> When CDRGES is called, a number of matrix "sizes" ("N's") and a
  50 *> number of matrix "TYPES" are specified.  For each size ("N")
  51 *> and each TYPE of matrix, a pair of matrices (A, B) will be generated
  52 *> and used for testing. For each matrix pair, the following 13 tests
  53 *> will be performed and compared with the threshold THRESH except
  54 *> the tests (5), (11) and (13).
  55 *>
  56 *>
  57 *> (1)   | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)
  58 *>
  59 *>
  60 *> (2)   | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)
  61 *>
  62 *>
  63 *> (3)   | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)
  64 *>
  65 *>
  66 *> (4)   | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)
  67 *>
  68 *> (5)   if A is in Schur form (i.e. triangular form) (no sorting of
  69 *>       eigenvalues)
  70 *>
  71 *> (6)   if eigenvalues = diagonal elements of the Schur form (S, T),
  72 *>       i.e., test the maximum over j of D(j)  where:
  73 *>
  74 *>                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
  75 *>           D(j) = ------------------------ + -----------------------
  76 *>                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
  77 *>
  78 *>       (no sorting of eigenvalues)
  79 *>
  80 *> (7)   | (A,B) - Q (S,T) Z' | / ( |(A,B)| n ulp )
  81 *>       (with sorting of eigenvalues).
  82 *>
  83 *> (8)   | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).
  84 *>
  85 *> (9)   | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).
  86 *>
  87 *> (10)  if A is in Schur form (i.e. quasi-triangular form)
  88 *>       (with sorting of eigenvalues).
  89 *>
  90 *> (11)  if eigenvalues = diagonal elements of the Schur form (S, T),
  91 *>       i.e. test the maximum over j of D(j)  where:
  92 *>
  93 *>                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
  94 *>           D(j) = ------------------------ + -----------------------
  95 *>                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
  96 *>
  97 *>       (with sorting of eigenvalues).
  98 *>
  99 *> (12)  if sorting worked and SDIM is the number of eigenvalues
 100 *>       which were CELECTed.
 101 *>
 102 *> Test Matrices
 103 *> =============
 104 *>
 105 *> The sizes of the test matrices are specified by an array
 106 *> NN(1:NSIZES); the value of each element NN(j) specifies one size.
 107 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
 108 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
 109 *> Currently, the list of possible types is:
 110 *>
 111 *> (1)  ( 0, 0 )         (a pair of zero matrices)
 112 *>
 113 *> (2)  ( I, 0 )         (an identity and a zero matrix)
 114 *>
 115 *> (3)  ( 0, I )         (an identity and a zero matrix)
 116 *>
 117 *> (4)  ( I, I )         (a pair of identity matrices)
 118 *>
 119 *>         t   t
 120 *> (5)  ( J , J  )       (a pair of transposed Jordan blocks)
 121 *>
 122 *>                                     t                ( I   0  )
 123 *> (6)  ( X, Y )         where  X = ( J   0  )  and Y = (      t )
 124 *>                                  ( 0   I  )          ( 0   J  )
 125 *>                       and I is a k x k identity and J a (k+1)x(k+1)
 126 *>                       Jordan block; k=(N-1)/2
 127 *>
 128 *> (7)  ( D, I )         where D is diag( 0, 1,..., N-1 ) (a diagonal
 129 *>                       matrix with those diagonal entries.)
 130 *> (8)  ( I, D )
 131 *>
 132 *> (9)  ( big*D, small*I ) where "big" is near overflow and small=1/big
 133 *>
 134 *> (10) ( small*D, big*I )
 135 *>
 136 *> (11) ( big*I, small*D )
 137 *>
 138 *> (12) ( small*I, big*D )
 139 *>
 140 *> (13) ( big*D, big*I )
 141 *>
 142 *> (14) ( small*D, small*I )
 143 *>
 144 *> (15) ( D1, D2 )        where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
 145 *>                        D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
 146 *>           t   t
 147 *> (16) Q ( J , J ) Z     where Q and Z are random orthogonal matrices.
 148 *>
 149 *> (17) Q ( T1, T2 ) Z    where T1 and T2 are upper triangular matrices
 150 *>                        with random O(1) entries above the diagonal
 151 *>                        and diagonal entries diag(T1) =
 152 *>                        ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
 153 *>                        ( 0, N-3, N-4,..., 1, 0, 0 )
 154 *>
 155 *> (18) Q ( T1, T2 ) Z    diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
 156 *>                        diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
 157 *>                        s = machine precision.
 158 *>
 159 *> (19) Q ( T1, T2 ) Z    diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
 160 *>                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
 161 *>
 162 *>                                                        N-5
 163 *> (20) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, 1, a, ..., a   =s, 0 )
 164 *>                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
 165 *>
 166 *> (21) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
 167 *>                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
 168 *>                        where r1,..., r(N-4) are random.
 169 *>
 170 *> (22) Q ( big*T1, small*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
 171 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
 172 *>
 173 *> (23) Q ( small*T1, big*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
 174 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
 175 *>
 176 *> (24) Q ( small*T1, small*T2 ) Z  diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
 177 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
 178 *>
 179 *> (25) Q ( big*T1, big*T2 ) Z      diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
 180 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
 181 *>
 182 *> (26) Q ( T1, T2 ) Z     where T1 and T2 are random upper-triangular
 183 *>                         matrices.
 184 *>
 185 *> \endverbatim
 186 *
 187 *  Arguments:
 188 *  ==========
 189 *
 190 *> \param[in] NSIZES
 191 *> \verbatim
 192 *>          NSIZES is INTEGER
 193 *>          The number of sizes of matrices to use.  If it is zero,
 194 *>          SDRGES does nothing.  NSIZES >= 0.
 195 *> \endverbatim
 196 *>
 197 *> \param[in] NN
 198 *> \verbatim
 199 *>          NN is INTEGER array, dimension (NSIZES)
 200 *>          An array containing the sizes to be used for the matrices.
 201 *>          Zero values will be skipped.  NN >= 0.
 202 *> \endverbatim
 203 *>
 204 *> \param[in] NTYPES
 205 *> \verbatim
 206 *>          NTYPES is INTEGER
 207 *>          The number of elements in DOTYPE.   If it is zero, SDRGES
 208 *>          does nothing.  It must be at least zero.  If it is MAXTYP+1
 209 *>          and NSIZES is 1, then an additional type, MAXTYP+1 is
 210 *>          defined, which is to use whatever matrix is in A on input.
 211 *>          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
 212 *>          DOTYPE(MAXTYP+1) is .TRUE. .
 213 *> \endverbatim
 214 *>
 215 *> \param[in] DOTYPE
 216 *> \verbatim
 217 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
 218 *>          If DOTYPE(j) is .TRUE., then for each size in NN a
 219 *>          matrix of that size and of type j will be generated.
 220 *>          If NTYPES is smaller than the maximum number of types
 221 *>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
 222 *>          MAXTYP will not be generated. If NTYPES is larger
 223 *>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
 224 *>          will be ignored.
 225 *> \endverbatim
 226 *>
 227 *> \param[in,out] ISEED
 228 *> \verbatim
 229 *>          ISEED is INTEGER array, dimension (4)
 230 *>          On entry ISEED specifies the seed of the random number
 231 *>          generator. The array elements should be between 0 and 4095;
 232 *>          if not they will be reduced mod 4096. Also, ISEED(4) must
 233 *>          be odd.  The random number generator uses a linear
 234 *>          congruential sequence limited to small integers, and so
 235 *>          should produce machine independent random numbers. The
 236 *>          values of ISEED are changed on exit, and can be used in the
 237 *>          next call to SDRGES to continue the same random number
 238 *>          sequence.
 239 *> \endverbatim
 240 *>
 241 *> \param[in] THRESH
 242 *> \verbatim
 243 *>          THRESH is REAL
 244 *>          A test will count as "failed" if the "error", computed as
 245 *>          described above, exceeds THRESH.  Note that the error is
 246 *>          scaled to be O(1), so THRESH should be a reasonably small
 247 *>          multiple of 1, e.g., 10 or 100.  In particular, it should
 248 *>          not depend on the precision (single vs. double) or the size
 249 *>          of the matrix.  THRESH >= 0.
 250 *> \endverbatim
 251 *>
 252 *> \param[in] NOUNIT
 253 *> \verbatim
 254 *>          NOUNIT is INTEGER
 255 *>          The FORTRAN unit number for printing out error messages
 256 *>          (e.g., if a routine returns IINFO not equal to 0.)
 257 *> \endverbatim
 258 *>
 259 *> \param[in,out] A
 260 *> \verbatim
 261 *>          A is COMPLEX array, dimension(LDA, max(NN))
 262 *>          Used to hold the original A matrix.  Used as input only
 263 *>          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
 264 *>          DOTYPE(MAXTYP+1)=.TRUE.
 265 *> \endverbatim
 266 *>
 267 *> \param[in] LDA
 268 *> \verbatim
 269 *>          LDA is INTEGER
 270 *>          The leading dimension of A, B, S, and T.
 271 *>          It must be at least 1 and at least max( NN ).
 272 *> \endverbatim
 273 *>
 274 *> \param[in,out] B
 275 *> \verbatim
 276 *>          B is COMPLEX array, dimension(LDA, max(NN))
 277 *>          Used to hold the original B matrix.  Used as input only
 278 *>          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
 279 *>          DOTYPE(MAXTYP+1)=.TRUE.
 280 *> \endverbatim
 281 *>
 282 *> \param[out] S
 283 *> \verbatim
 284 *>          S is COMPLEX array, dimension (LDA, max(NN))
 285 *>          The Schur form matrix computed from A by CGGES.  On exit, S
 286 *>          contains the Schur form matrix corresponding to the matrix
 287 *>          in A.
 288 *> \endverbatim
 289 *>
 290 *> \param[out] T
 291 *> \verbatim
 292 *>          T is COMPLEX array, dimension (LDA, max(NN))
 293 *>          The upper triangular matrix computed from B by CGGES.
 294 *> \endverbatim
 295 *>
 296 *> \param[out] Q
 297 *> \verbatim
 298 *>          Q is COMPLEX array, dimension (LDQ, max(NN))
 299 *>          The (left) orthogonal matrix computed by CGGES.
 300 *> \endverbatim
 301 *>
 302 *> \param[in] LDQ
 303 *> \verbatim
 304 *>          LDQ is INTEGER
 305 *>          The leading dimension of Q and Z. It must
 306 *>          be at least 1 and at least max( NN ).
 307 *> \endverbatim
 308 *>
 309 *> \param[out] Z
 310 *> \verbatim
 311 *>          Z is COMPLEX array, dimension( LDQ, max(NN) )
 312 *>          The (right) orthogonal matrix computed by CGGES.
 313 *> \endverbatim
 314 *>
 315 *> \param[out] ALPHA
 316 *> \verbatim
 317 *>          ALPHA is COMPLEX array, dimension (max(NN))
 318 *> \endverbatim
 319 *>
 320 *> \param[out] BETA
 321 *> \verbatim
 322 *>          BETA is COMPLEX array, dimension (max(NN))
 323 *>
 324 *>          The generalized eigenvalues of (A,B) computed by CGGES.
 325 *>          ALPHA(k) / BETA(k) is the k-th generalized eigenvalue of A
 326 *>          and B.
 327 *> \endverbatim
 328 *>
 329 *> \param[out] WORK
 330 *> \verbatim
 331 *>          WORK is COMPLEX array, dimension (LWORK)
 332 *> \endverbatim
 333 *>
 334 *> \param[in] LWORK
 335 *> \verbatim
 336 *>          LWORK is INTEGER
 337 *>          The dimension of the array WORK.  LWORK >= 3*N*N.
 338 *> \endverbatim
 339 *>
 340 *> \param[out] RWORK
 341 *> \verbatim
 342 *>          RWORK is REAL array, dimension ( 8*N )
 343 *>          Real workspace.
 344 *> \endverbatim
 345 *>
 346 *> \param[out] RESULT
 347 *> \verbatim
 348 *>          RESULT is REAL array, dimension (15)
 349 *>          The values computed by the tests described above.
 350 *>          The values are currently limited to 1/ulp, to avoid overflow.
 351 *> \endverbatim
 352 *>
 353 *> \param[out] BWORK
 354 *> \verbatim
 355 *>          BWORK is LOGICAL array, dimension (N)
 356 *> \endverbatim
 357 *>
 358 *> \param[out] INFO
 359 *> \verbatim
 360 *>          INFO is INTEGER
 361 *>          = 0:  successful exit
 362 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 363 *>          > 0:  A routine returned an error code.  INFO is the
 364 *>                absolute value of the INFO value returned.
 365 *> \endverbatim
 366 *
 367 *  Authors:
 368 *  ========
 369 *
 370 *> \author Univ. of Tennessee
 371 *> \author Univ. of California Berkeley
 372 *> \author Univ. of Colorado Denver
 373 *> \author NAG Ltd.
 374 *
 375 *> \date June 2016
 376 *
 377 *> \ingroup complex_eig
 378 *
 379 *  =====================================================================
 380       SUBROUTINE CDRGES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
 381      $                   NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHA,
 382      $                   BETA, WORK, LWORK, RWORK, RESULT, BWORK, INFO )
 383 *
 384 *  -- LAPACK test routine (version 3.6.1) --
 385 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 386 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 387 *     June 2016
 388 *
 389 *     .. Scalar Arguments ..
 390       INTEGER            INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
 391       REAL               THRESH
 392 *     ..
 393 *     .. Array Arguments ..
 394       LOGICAL            BWORK( * ), DOTYPE( * )
 395       INTEGER            ISEED( 4 ), NN( * )
 396       REAL               RESULT( 13 ), RWORK( * )
 397       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDA, * ),
 398      $                   BETA( * ), Q( LDQ, * ), S( LDA, * ),
 399      $                   T( LDA, * ), WORK( * ), Z( LDQ, * )
 400 *     ..
 401 *
 402 *  =====================================================================
 403 *
 404 *     .. Parameters ..
 405       REAL               ZERO, ONE
 406       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 407       COMPLEX            CZERO, CONE
 408       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 409      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
 410       INTEGER            MAXTYP
 411       PARAMETER          ( MAXTYP = 26 )
 412 *     ..
 413 *     .. Local Scalars ..
 414       LOGICAL            BADNN, ILABAD
 415       CHARACTER          SORT
 416       INTEGER            I, IADD, IINFO, IN, ISORT, J, JC, JR, JSIZE,
 417      $                   JTYPE, KNTEIG, MAXWRK, MINWRK, MTYPES, N, N1,
 418      $                   NB, NERRS, NMATS, NMAX, NTEST, NTESTT, RSUB,
 419      $                   SDIM
 420       REAL               SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
 421       COMPLEX            CTEMP, X
 422 *     ..
 423 *     .. Local Arrays ..
 424       LOGICAL            LASIGN( MAXTYP ), LBSIGN( MAXTYP )
 425       INTEGER            IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
 426      $                   KATYPE( MAXTYP ), KAZERO( MAXTYP ),
 427      $                   KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
 428      $                   KBZERO( MAXTYP ), KCLASS( MAXTYP ),
 429      $                   KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
 430       REAL               RMAGN( 0: 3 )
 431 *     ..
 432 *     .. External Functions ..
 433       LOGICAL            CLCTES
 434       INTEGER            ILAENV
 435       REAL               SLAMCH
 436       COMPLEX            CLARND
 437       EXTERNAL           CLCTES, ILAENV, SLAMCH, CLARND
 438 *     ..
 439 *     .. External Subroutines ..
 440       EXTERNAL           ALASVM, CGET51, CGET54, CGGES, CLACPY, CLARFG,
 441      $                   CLASET, CLATM4, CUNM2R, SLABAD, XERBLA
 442 *     ..
 443 *     .. Intrinsic Functions ..
 444       INTRINSIC          ABS, AIMAG, CONJG, MAX, MIN, REAL, SIGN
 445 *     ..
 446 *     .. Statement Functions ..
 447       REAL               ABS1
 448 *     ..
 449 *     .. Statement Function definitions ..
 450       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
 451 *     ..
 452 *     .. Data statements ..
 453       DATA               KCLASS / 15*1, 10*2, 1*3 /
 454       DATA               KZ1 / 0, 1, 2, 1, 3, 3 /
 455       DATA               KZ2 / 0, 0, 1, 2, 1, 1 /
 456       DATA               KADD / 0, 0, 0, 0, 3, 2 /
 457       DATA               KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
 458      $                   4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
 459       DATA               KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
 460      $                   1, 1, -4, 2, -4, 8*8, 0 /
 461       DATA               KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
 462      $                   4*5, 4*3, 1 /
 463       DATA               KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
 464      $                   4*6, 4*4, 1 /
 465       DATA               KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
 466      $                   2, 1 /
 467       DATA               KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
 468      $                   2, 1 /
 469       DATA               KTRIAN / 16*0, 10*1 /
 470       DATA               LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE.,
 471      $                   2*.FALSE., 3*.TRUE., .FALSE., .TRUE.,
 472      $                   3*.FALSE., 5*.TRUE., .FALSE. /
 473       DATA               LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE.,
 474      $                   2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE.,
 475      $                   9*.FALSE. /
 476 *     ..
 477 *     .. Executable Statements ..
 478 *
 479 *     Check for errors
 480 *
 481       INFO = 0
 482 *
 483       BADNN = .FALSE.
 484       NMAX = 1
 485       DO 10 J = 1, NSIZES
 486          NMAX = MAX( NMAX, NN( J ) )
 487          IF( NN( J ).LT.0 )
 488      $      BADNN = .TRUE.
 489    10 CONTINUE
 490 *
 491       IF( NSIZES.LT.0 ) THEN
 492          INFO = -1
 493       ELSE IF( BADNN ) THEN
 494          INFO = -2
 495       ELSE IF( NTYPES.LT.0 ) THEN
 496          INFO = -3
 497       ELSE IF( THRESH.LT.ZERO ) THEN
 498          INFO = -6
 499       ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
 500          INFO = -9
 501       ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN
 502          INFO = -14
 503       END IF
 504 *
 505 *     Compute workspace
 506 *      (Note: Comments in the code beginning "Workspace:" describe the
 507 *       minimal amount of workspace needed at that point in the code,
 508 *       as well as the preferred amount for good performance.
 509 *       NB refers to the optimal block size for the immediately
 510 *       following subroutine, as returned by ILAENV.
 511 *
 512       MINWRK = 1
 513       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
 514          MINWRK = 3*NMAX*NMAX
 515          NB = MAX( 1, ILAENV( 1, 'CGEQRF', ' ', NMAX, NMAX, -1, -1 ),
 516      $        ILAENV( 1, 'CUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ),
 517      $        ILAENV( 1, 'CUNGQR', ' ', NMAX, NMAX, NMAX, -1 ) )
 518          MAXWRK = MAX( NMAX+NMAX*NB, 3*NMAX*NMAX )
 519          WORK( 1 ) = MAXWRK
 520       END IF
 521 *
 522       IF( LWORK.LT.MINWRK )
 523      $   INFO = -19
 524 *
 525       IF( INFO.NE.0 ) THEN
 526          CALL XERBLA( 'CDRGES', -INFO )
 527          RETURN
 528       END IF
 529 *
 530 *     Quick return if possible
 531 *
 532       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
 533      $   RETURN
 534 *
 535       ULP = SLAMCH( 'Precision' )
 536       SAFMIN = SLAMCH( 'Safe minimum' )
 537       SAFMIN = SAFMIN / ULP
 538       SAFMAX = ONE / SAFMIN
 539       CALL SLABAD( SAFMIN, SAFMAX )
 540       ULPINV = ONE / ULP
 541 *
 542 *     The values RMAGN(2:3) depend on N, see below.
 543 *
 544       RMAGN( 0 ) = ZERO
 545       RMAGN( 1 ) = ONE
 546 *
 547 *     Loop over matrix sizes
 548 *
 549       NTESTT = 0
 550       NERRS = 0
 551       NMATS = 0
 552 *
 553       DO 190 JSIZE = 1, NSIZES
 554          N = NN( JSIZE )
 555          N1 = MAX( 1, N )
 556          RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 )
 557          RMAGN( 3 ) = SAFMIN*ULPINV*REAL( N1 )
 558 *
 559          IF( NSIZES.NE.1 ) THEN
 560             MTYPES = MIN( MAXTYP, NTYPES )
 561          ELSE
 562             MTYPES = MIN( MAXTYP+1, NTYPES )
 563          END IF
 564 *
 565 *        Loop over matrix types
 566 *
 567          DO 180 JTYPE = 1, MTYPES
 568             IF( .NOT.DOTYPE( JTYPE ) )
 569      $         GO TO 180
 570             NMATS = NMATS + 1
 571             NTEST = 0
 572 *
 573 *           Save ISEED in case of an error.
 574 *
 575             DO 20 J = 1, 4
 576                IOLDSD( J ) = ISEED( J )
 577    20       CONTINUE
 578 *
 579 *           Initialize RESULT
 580 *
 581             DO 30 J = 1, 13
 582                RESULT( J ) = ZERO
 583    30       CONTINUE
 584 *
 585 *           Generate test matrices A and B
 586 *
 587 *           Description of control parameters:
 588 *
 589 *           KCLASS: =1 means w/o rotation, =2 means w/ rotation,
 590 *                   =3 means random.
 591 *           KATYPE: the "type" to be passed to CLATM4 for computing A.
 592 *           KAZERO: the pattern of zeros on the diagonal for A:
 593 *                   =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
 594 *                   =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
 595 *                   =6: ( 0, 1, 0, xxx, 0 ).  (xxx means a string of
 596 *                   non-zero entries.)
 597 *           KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
 598 *                   =2: large, =3: small.
 599 *           LASIGN: .TRUE. if the diagonal elements of A are to be
 600 *                   multiplied by a random magnitude 1 number.
 601 *           KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
 602 *           KTRIAN: =0: don't fill in the upper triangle, =1: do.
 603 *           KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
 604 *           RMAGN: used to implement KAMAGN and KBMAGN.
 605 *
 606             IF( MTYPES.GT.MAXTYP )
 607      $         GO TO 110
 608             IINFO = 0
 609             IF( KCLASS( JTYPE ).LT.3 ) THEN
 610 *
 611 *              Generate A (w/o rotation)
 612 *
 613                IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
 614                   IN = 2*( ( N-1 ) / 2 ) + 1
 615                   IF( IN.NE.N )
 616      $               CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
 617                ELSE
 618                   IN = N
 619                END IF
 620                CALL CLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
 621      $                      KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
 622      $                      RMAGN( KAMAGN( JTYPE ) ), ULP,
 623      $                      RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2,
 624      $                      ISEED, A, LDA )
 625                IADD = KADD( KAZERO( JTYPE ) )
 626                IF( IADD.GT.0 .AND. IADD.LE.N )
 627      $            A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
 628 *
 629 *              Generate B (w/o rotation)
 630 *
 631                IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
 632                   IN = 2*( ( N-1 ) / 2 ) + 1
 633                   IF( IN.NE.N )
 634      $               CALL CLASET( 'Full', N, N, CZERO, CZERO, B, LDA )
 635                ELSE
 636                   IN = N
 637                END IF
 638                CALL CLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
 639      $                      KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
 640      $                      RMAGN( KBMAGN( JTYPE ) ), ONE,
 641      $                      RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2,
 642      $                      ISEED, B, LDA )
 643                IADD = KADD( KBZERO( JTYPE ) )
 644                IF( IADD.NE.0 .AND. IADD.LE.N )
 645      $            B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
 646 *
 647                IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
 648 *
 649 *                 Include rotations
 650 *
 651 *                 Generate Q, Z as Householder transformations times
 652 *                 a diagonal matrix.
 653 *
 654                   DO 50 JC = 1, N - 1
 655                      DO 40 JR = JC, N
 656                         Q( JR, JC ) = CLARND( 3, ISEED )
 657                         Z( JR, JC ) = CLARND( 3, ISEED )
 658    40                CONTINUE
 659                      CALL CLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1,
 660      $                            WORK( JC ) )
 661                      WORK( 2*N+JC ) = SIGN( ONE, REAL( Q( JC, JC ) ) )
 662                      Q( JC, JC ) = CONE
 663                      CALL CLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1,
 664      $                            WORK( N+JC ) )
 665                      WORK( 3*N+JC ) = SIGN( ONE, REAL( Z( JC, JC ) ) )
 666                      Z( JC, JC ) = CONE
 667    50             CONTINUE
 668                   CTEMP = CLARND( 3, ISEED )
 669                   Q( N, N ) = CONE
 670                   WORK( N ) = CZERO
 671                   WORK( 3*N ) = CTEMP / ABS( CTEMP )
 672                   CTEMP = CLARND( 3, ISEED )
 673                   Z( N, N ) = CONE
 674                   WORK( 2*N ) = CZERO
 675                   WORK( 4*N ) = CTEMP / ABS( CTEMP )
 676 *
 677 *                 Apply the diagonal matrices
 678 *
 679                   DO 70 JC = 1, N
 680                      DO 60 JR = 1, N
 681                         A( JR, JC ) = WORK( 2*N+JR )*
 682      $                                CONJG( WORK( 3*N+JC ) )*
 683      $                                A( JR, JC )
 684                         B( JR, JC ) = WORK( 2*N+JR )*
 685      $                                CONJG( WORK( 3*N+JC ) )*
 686      $                                B( JR, JC )
 687    60                CONTINUE
 688    70             CONTINUE
 689                   CALL CUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
 690      $                         LDA, WORK( 2*N+1 ), IINFO )
 691                   IF( IINFO.NE.0 )
 692      $               GO TO 100
 693                   CALL CUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
 694      $                         A, LDA, WORK( 2*N+1 ), IINFO )
 695                   IF( IINFO.NE.0 )
 696      $               GO TO 100
 697                   CALL CUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
 698      $                         LDA, WORK( 2*N+1 ), IINFO )
 699                   IF( IINFO.NE.0 )
 700      $               GO TO 100
 701                   CALL CUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
 702      $                         B, LDA, WORK( 2*N+1 ), IINFO )
 703                   IF( IINFO.NE.0 )
 704      $               GO TO 100
 705                END IF
 706             ELSE
 707 *
 708 *              Random matrices
 709 *
 710                DO 90 JC = 1, N
 711                   DO 80 JR = 1, N
 712                      A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
 713      $                             CLARND( 4, ISEED )
 714                      B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
 715      $                             CLARND( 4, ISEED )
 716    80             CONTINUE
 717    90          CONTINUE
 718             END IF
 719 *
 720   100       CONTINUE
 721 *
 722             IF( IINFO.NE.0 ) THEN
 723                WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
 724      $            IOLDSD
 725                INFO = ABS( IINFO )
 726                RETURN
 727             END IF
 728 *
 729   110       CONTINUE
 730 *
 731             DO 120 I = 1, 13
 732                RESULT( I ) = -ONE
 733   120       CONTINUE
 734 *
 735 *           Test with and without sorting of eigenvalues
 736 *
 737             DO 150 ISORT = 0, 1
 738                IF( ISORT.EQ.0 ) THEN
 739                   SORT = 'N'
 740                   RSUB = 0
 741                ELSE
 742                   SORT = 'S'
 743                   RSUB = 5
 744                END IF
 745 *
 746 *              Call CGGES to compute H, T, Q, Z, alpha, and beta.
 747 *
 748                CALL CLACPY( 'Full', N, N, A, LDA, S, LDA )
 749                CALL CLACPY( 'Full', N, N, B, LDA, T, LDA )
 750                NTEST = 1 + RSUB + ISORT
 751                RESULT( 1+RSUB+ISORT ) = ULPINV
 752                CALL CGGES( 'V', 'V', SORT, CLCTES, N, S, LDA, T, LDA,
 753      $                     SDIM, ALPHA, BETA, Q, LDQ, Z, LDQ, WORK,
 754      $                     LWORK, RWORK, BWORK, IINFO )
 755                IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
 756                   RESULT( 1+RSUB+ISORT ) = ULPINV
 757                   WRITE( NOUNIT, FMT = 9999 )'CGGES', IINFO, N, JTYPE,
 758      $               IOLDSD
 759                   INFO = ABS( IINFO )
 760                   GO TO 160
 761                END IF
 762 *
 763                NTEST = 4 + RSUB
 764 *
 765 *              Do tests 1--4 (or tests 7--9 when reordering )
 766 *
 767                IF( ISORT.EQ.0 ) THEN
 768                   CALL CGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ,
 769      $                         WORK, RWORK, RESULT( 1 ) )
 770                   CALL CGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ,
 771      $                         WORK, RWORK, RESULT( 2 ) )
 772                ELSE
 773                   CALL CGET54( N, A, LDA, B, LDA, S, LDA, T, LDA, Q,
 774      $                         LDQ, Z, LDQ, WORK, RESULT( 2+RSUB ) )
 775                END IF
 776 *
 777                CALL CGET51( 3, N, B, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK,
 778      $                      RWORK, RESULT( 3+RSUB ) )
 779                CALL CGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK,
 780      $                      RWORK, RESULT( 4+RSUB ) )
 781 *
 782 *              Do test 5 and 6 (or Tests 10 and 11 when reordering):
 783 *              check Schur form of A and compare eigenvalues with
 784 *              diagonals.
 785 *
 786                NTEST = 6 + RSUB
 787                TEMP1 = ZERO
 788 *
 789                DO 130 J = 1, N
 790                   ILABAD = .FALSE.
 791                   TEMP2 = ( ABS1( ALPHA( J )-S( J, J ) ) /
 792      $                    MAX( SAFMIN, ABS1( ALPHA( J ) ), ABS1( S( J,
 793      $                    J ) ) )+ABS1( BETA( J )-T( J, J ) ) /
 794      $                    MAX( SAFMIN, ABS1( BETA( J ) ), ABS1( T( J,
 795      $                    J ) ) ) ) / ULP
 796 *
 797                   IF( J.LT.N ) THEN
 798                      IF( S( J+1, J ).NE.ZERO ) THEN
 799                         ILABAD = .TRUE.
 800                         RESULT( 5+RSUB ) = ULPINV
 801                      END IF
 802                   END IF
 803                   IF( J.GT.1 ) THEN
 804                      IF( S( J, J-1 ).NE.ZERO ) THEN
 805                         ILABAD = .TRUE.
 806                         RESULT( 5+RSUB ) = ULPINV
 807                      END IF
 808                   END IF
 809                   TEMP1 = MAX( TEMP1, TEMP2 )
 810                   IF( ILABAD ) THEN
 811                      WRITE( NOUNIT, FMT = 9998 )J, N, JTYPE, IOLDSD
 812                   END IF
 813   130          CONTINUE
 814                RESULT( 6+RSUB ) = TEMP1
 815 *
 816                IF( ISORT.GE.1 ) THEN
 817 *
 818 *                 Do test 12
 819 *
 820                   NTEST = 12
 821                   RESULT( 12 ) = ZERO
 822                   KNTEIG = 0
 823                   DO 140 I = 1, N
 824                      IF( CLCTES( ALPHA( I ), BETA( I ) ) )
 825      $                  KNTEIG = KNTEIG + 1
 826   140             CONTINUE
 827                   IF( SDIM.NE.KNTEIG )
 828      $               RESULT( 13 ) = ULPINV
 829                END IF
 830 *
 831   150       CONTINUE
 832 *
 833 *           End of Loop -- Check for RESULT(j) > THRESH
 834 *
 835   160       CONTINUE
 836 *
 837             NTESTT = NTESTT + NTEST
 838 *
 839 *           Print out tests which fail.
 840 *
 841             DO 170 JR = 1, NTEST
 842                IF( RESULT( JR ).GE.THRESH ) THEN
 843 *
 844 *                 If this is the first test to fail,
 845 *                 print a header to the data file.
 846 *
 847                   IF( NERRS.EQ.0 ) THEN
 848                      WRITE( NOUNIT, FMT = 9997 )'CGS'
 849 *
 850 *                    Matrix types
 851 *
 852                      WRITE( NOUNIT, FMT = 9996 )
 853                      WRITE( NOUNIT, FMT = 9995 )
 854                      WRITE( NOUNIT, FMT = 9994 )'Unitary'
 855 *
 856 *                    Tests performed
 857 *
 858                      WRITE( NOUNIT, FMT = 9993 )'unitary', '''',
 859      $                  'transpose', ( '''', J = 1, 8 )
 860 *
 861                   END IF
 862                   NERRS = NERRS + 1
 863                   IF( RESULT( JR ).LT.10000.0 ) THEN
 864                      WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
 865      $                  RESULT( JR )
 866                   ELSE
 867                      WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
 868      $                  RESULT( JR )
 869                   END IF
 870                END IF
 871   170       CONTINUE
 872 *
 873   180    CONTINUE
 874   190 CONTINUE
 875 *
 876 *     Summary
 877 *
 878       CALL ALASVM( 'CGS', NOUNIT, NERRS, NTESTT, 0 )
 879 *
 880       WORK( 1 ) = MAXWRK
 881 *
 882       RETURN
 883 *
 884  9999 FORMAT( ' CDRGES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
 885      $      I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' )
 886 *
 887  9998 FORMAT( ' CDRGES: S not in Schur form at eigenvalue ', I6, '.',
 888      $      / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
 889      $      I5, ')' )
 890 *
 891  9997 FORMAT( / 1X, A3, ' -- Complex Generalized Schur from problem ',
 892      $      'driver' )
 893 *
 894  9996 FORMAT( ' Matrix types (see CDRGES for details): ' )
 895 *
 896  9995 FORMAT( ' Special Matrices:', 23X,
 897      $      '(J''=transposed Jordan block)',
 898      $      / '   1=(0,0)  2=(I,0)  3=(0,I)  4=(I,I)  5=(J'',J'')  ',
 899      $      '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices:  ( ',
 900      $      'D=diag(0,1,2,...) )', / '   7=(D,I)   9=(large*D, small*I',
 901      $      ')  11=(large*I, small*D)  13=(large*D, large*I)', /
 902      $      '   8=(I,D)  10=(small*D, large*I)  12=(small*I, large*D) ',
 903      $      ' 14=(small*D, small*I)', / '  15=(D, reversed D)' )
 904  9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
 905      $      / '  16=Transposed Jordan Blocks             19=geometric ',
 906      $      'alpha, beta=0,1', / '  17=arithm. alpha&beta             ',
 907      $      '      20=arithmetic alpha, beta=0,1', / '  18=clustered ',
 908      $      'alpha, beta=0,1            21=random alpha, beta=0,1',
 909      $      / ' Large & Small Matrices:', / '  22=(large, small)   ',
 910      $      '23=(small,large)    24=(small,small)    25=(large,large)',
 911      $      / '  26=random O(1) matrices.' )
 912 *
 913  9993 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
 914      $      'Q and Z are ', A, ',', / 19X,
 915      $      'l and r are the appropriate left and right', / 19X,
 916      $      'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A,
 917      $      ' means ', A, '.)', / ' Without ordering: ',
 918      $      / '  1 = | A - Q S Z', A,
 919      $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
 920      $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
 921      $      ' | / ( n ulp )             4 = | I - ZZ', A,
 922      $      ' | / ( n ulp )', / '  5 = A is in Schur form S',
 923      $      / '  6 = difference between (alpha,beta)',
 924      $      ' and diagonals of (S,T)', / ' With ordering: ',
 925      $      / '  7 = | (A,B) - Q (S,T) Z', A, ' | / ( |(A,B)| n ulp )',
 926      $      / '  8 = | I - QQ', A,
 927      $      ' | / ( n ulp )             9 = | I - ZZ', A,
 928      $      ' | / ( n ulp )', / ' 10 = A is in Schur form S',
 929      $      / ' 11 = difference between (alpha,beta) and diagonals',
 930      $      ' of (S,T)', / ' 12 = SDIM is the correct number of ',
 931      $      'selected eigenvalues', / )
 932  9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
 933      $      4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
 934  9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
 935      $      4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 )
 936 *
 937 *     End of CDRGES
 938 *
 939       END