1*> \brief \b ZDRGES 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 ZDRGES( 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* DOUBLE PRECISION THRESH 18* .. 19* .. Array Arguments .. 20* LOGICAL BWORK( * ), DOTYPE( * ) 21* INTEGER ISEED( 4 ), NN( * ) 22* DOUBLE PRECISION RESULT( 13 ), RWORK( * ) 23* COMPLEX*16 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*> ZDRGES checks the nonsymmetric generalized eigenvalue (Schur form) 35*> problem driver ZGGES. 36*> 37*> ZGGES 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 ZDRGES 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 threshhold 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*> DDRGES 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, DDRGES 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 DDRGES to continue the same random number 238*> sequence. 239*> \endverbatim 240*> 241*> \param[in] THRESH 242*> \verbatim 243*> THRESH is DOUBLE PRECISION 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*16 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*16 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*16 array, dimension (LDA, max(NN)) 285*> The Schur form matrix computed from A by ZGGES. 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*16 array, dimension (LDA, max(NN)) 293*> The upper triangular matrix computed from B by ZGGES. 294*> \endverbatim 295*> 296*> \param[out] Q 297*> \verbatim 298*> Q is COMPLEX*16 array, dimension (LDQ, max(NN)) 299*> The (left) orthogonal matrix computed by ZGGES. 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*16 array, dimension( LDQ, max(NN) ) 312*> The (right) orthogonal matrix computed by ZGGES. 313*> \endverbatim 314*> 315*> \param[out] ALPHA 316*> \verbatim 317*> ALPHA is COMPLEX*16 array, dimension (max(NN)) 318*> \endverbatim 319*> 320*> \param[out] BETA 321*> \verbatim 322*> BETA is COMPLEX*16 array, dimension (max(NN)) 323*> 324*> The generalized eigenvalues of (A,B) computed by ZGGES. 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*16 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 DOUBLE PRECISION array, dimension ( 8*N ) 343*> Real workspace. 344*> \endverbatim 345*> 346*> \param[out] RESULT 347*> \verbatim 348*> RESULT is DOUBLE PRECISION 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 November 2011 376* 377*> \ingroup complex16_eig 378* 379* ===================================================================== 380 SUBROUTINE ZDRGES( 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.4.0) -- 385* -- LAPACK is a software package provided by Univ. of Tennessee, -- 386* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 387* November 2011 388* 389* .. Scalar Arguments .. 390 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES 391 DOUBLE PRECISION THRESH 392* .. 393* .. Array Arguments .. 394 LOGICAL BWORK( * ), DOTYPE( * ) 395 INTEGER ISEED( 4 ), NN( * ) 396 DOUBLE PRECISION RESULT( 13 ), RWORK( * ) 397 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDA, * ), 398 $ BETA( * ), Q( LDQ, * ), S( LDA, * ), 399 $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 400* .. 401* 402* ===================================================================== 403* 404* .. Parameters .. 405 DOUBLE PRECISION ZERO, ONE 406 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 407 COMPLEX*16 CZERO, CONE 408 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 409 $ CONE = ( 1.0D+0, 0.0D+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 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV 421 COMPLEX*16 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 DOUBLE PRECISION RMAGN( 0: 3 ) 431* .. 432* .. External Functions .. 433 LOGICAL ZLCTES 434 INTEGER ILAENV 435 DOUBLE PRECISION DLAMCH 436 COMPLEX*16 ZLARND 437 EXTERNAL ZLCTES, ILAENV, DLAMCH, ZLARND 438* .. 439* .. External Subroutines .. 440 EXTERNAL ALASVM, DLABAD, XERBLA, ZGET51, ZGET54, ZGGES, 441 $ ZLACPY, ZLARFG, ZLASET, ZLATM4, ZUNM2R 442* .. 443* .. Intrinsic Functions .. 444 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SIGN 445* .. 446* .. Statement Functions .. 447 DOUBLE PRECISION ABS1 448* .. 449* .. Statement Function definitions .. 450 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( 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, 'ZGEQRF', ' ', NMAX, NMAX, -1, -1 ), 516 $ ILAENV( 1, 'ZUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ), 517 $ ILAENV( 1, 'ZUNGQR', ' ', 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( 'ZDRGES', -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 = DLAMCH( 'Precision' ) 536 SAFMIN = DLAMCH( 'Safe minimum' ) 537 SAFMIN = SAFMIN / ULP 538 SAFMAX = ONE / SAFMIN 539 CALL DLABAD( 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 / DBLE( N1 ) 557 RMAGN( 3 ) = SAFMIN*ULPINV*DBLE( 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* KZLASS: =1 means w/o rotation, =2 means w/ rotation, 590* =3 means random. 591* KATYPE: the "type" to be passed to ZLATM4 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 ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 617 ELSE 618 IN = N 619 END IF 620 CALL ZLATM4( 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 ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 635 ELSE 636 IN = N 637 END IF 638 CALL ZLATM4( 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 ) = ZLARND( 3, ISEED ) 657 Z( JR, JC ) = ZLARND( 3, ISEED ) 658 40 CONTINUE 659 CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 660 $ WORK( JC ) ) 661 WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) ) 662 Q( JC, JC ) = CONE 663 CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 664 $ WORK( N+JC ) ) 665 WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) ) 666 Z( JC, JC ) = CONE 667 50 CONTINUE 668 CTEMP = ZLARND( 3, ISEED ) 669 Q( N, N ) = CONE 670 WORK( N ) = CZERO 671 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 672 CTEMP = ZLARND( 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 $ DCONJG( WORK( 3*N+JC ) )* 683 $ A( JR, JC ) 684 B( JR, JC ) = WORK( 2*N+JR )* 685 $ DCONJG( WORK( 3*N+JC ) )* 686 $ B( JR, JC ) 687 60 CONTINUE 688 70 CONTINUE 689 CALL ZUNM2R( '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 ZUNM2R( '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 ZUNM2R( '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 ZUNM2R( '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 $ ZLARND( 4, ISEED ) 714 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 715 $ ZLARND( 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 ZGGES to compute H, T, Q, Z, alpha, and beta. 747* 748 CALL ZLACPY( 'Full', N, N, A, LDA, S, LDA ) 749 CALL ZLACPY( 'Full', N, N, B, LDA, T, LDA ) 750 NTEST = 1 + RSUB + ISORT 751 RESULT( 1+RSUB+ISORT ) = ULPINV 752 CALL ZGGES( 'V', 'V', SORT, ZLCTES, 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 )'ZGGES', 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 ZGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ, 769 $ WORK, RWORK, RESULT( 1 ) ) 770 CALL ZGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ, 771 $ WORK, RWORK, RESULT( 2 ) ) 772 ELSE 773 CALL ZGET54( N, A, LDA, B, LDA, S, LDA, T, LDA, Q, 774 $ LDQ, Z, LDQ, WORK, RESULT( 2+RSUB ) ) 775 END IF 776* 777 CALL ZGET51( 3, N, B, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK, 778 $ RWORK, RESULT( 3+RSUB ) ) 779 CALL ZGET51( 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( ZLCTES( 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 )'ZGS' 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.0D0 ) 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( 'ZGS', NOUNIT, NERRS, NTESTT, 0 ) 879* 880 WORK( 1 ) = MAXWRK 881* 882 RETURN 883* 884 9999 FORMAT( ' ZDRGES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 885 $ I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' ) 886* 887 9998 FORMAT( ' ZDRGES: 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 ZDRGES 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, D10.3 ) 936* 937* End of ZDRGES 938* 939 END 940