1*> \brief \b ZDRGES3 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 ZDRGES3( 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*> ZDRGES3 checks the nonsymmetric generalized eigenvalue (Schur form) 35*> problem driver ZGGES3. 36*> 37*> ZGGES3 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 ZDRGES3 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*> DDRGES3 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, DDRGES3 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 DDRGES3 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 ZGGES3. 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 ZGGES3. 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 ZGGES3. 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 ZGGES3. 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 ZGGES3. 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 February 2015 376* 377*> \ingroup complex16_eig 378* 379* ===================================================================== 380 SUBROUTINE ZDRGES3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 381 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHA, 382 $ BETA, WORK, LWORK, RWORK, RESULT, BWORK, 383 $ INFO ) 384* 385* -- LAPACK test routine (version 3.6.0) -- 386* -- LAPACK is a software package provided by Univ. of Tennessee, -- 387* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 388* February 2015 389* 390* .. Scalar Arguments .. 391 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES 392 DOUBLE PRECISION THRESH 393* .. 394* .. Array Arguments .. 395 LOGICAL BWORK( * ), DOTYPE( * ) 396 INTEGER ISEED( 4 ), NN( * ) 397 DOUBLE PRECISION RESULT( 13 ), RWORK( * ) 398 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDA, * ), 399 $ BETA( * ), Q( LDQ, * ), S( LDA, * ), 400 $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 401* .. 402* 403* ===================================================================== 404* 405* .. Parameters .. 406 DOUBLE PRECISION ZERO, ONE 407 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 408 COMPLEX*16 CZERO, CONE 409 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 410 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 411 INTEGER MAXTYP 412 PARAMETER ( MAXTYP = 26 ) 413* .. 414* .. Local Scalars .. 415 LOGICAL BADNN, ILABAD 416 CHARACTER SORT 417 INTEGER I, IADD, IINFO, IN, ISORT, J, JC, JR, JSIZE, 418 $ JTYPE, KNTEIG, MAXWRK, MINWRK, MTYPES, N, N1, 419 $ NB, NERRS, NMATS, NMAX, NTEST, NTESTT, RSUB, 420 $ SDIM 421 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV 422 COMPLEX*16 CTEMP, X 423* .. 424* .. Local Arrays .. 425 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP ) 426 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 427 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 428 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 429 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 430 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 431 DOUBLE PRECISION RMAGN( 0: 3 ) 432* .. 433* .. External Functions .. 434 LOGICAL ZLCTES 435 INTEGER ILAENV 436 DOUBLE PRECISION DLAMCH 437 COMPLEX*16 ZLARND 438 EXTERNAL ZLCTES, ILAENV, DLAMCH, ZLARND 439* .. 440* .. External Subroutines .. 441 EXTERNAL ALASVM, DLABAD, XERBLA, ZGET51, ZGET54, ZGGES3, 442 $ ZLACPY, ZLARFG, ZLASET, ZLATM4, ZUNM2R 443* .. 444* .. Intrinsic Functions .. 445 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SIGN 446* .. 447* .. Statement Functions .. 448 DOUBLE PRECISION ABS1 449* .. 450* .. Statement Function definitions .. 451 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) 452* .. 453* .. Data statements .. 454 DATA KCLASS / 15*1, 10*2, 1*3 / 455 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 456 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 457 DATA KADD / 0, 0, 0, 0, 3, 2 / 458 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 459 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 460 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 461 $ 1, 1, -4, 2, -4, 8*8, 0 / 462 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 463 $ 4*5, 4*3, 1 / 464 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 465 $ 4*6, 4*4, 1 / 466 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 467 $ 2, 1 / 468 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 469 $ 2, 1 / 470 DATA KTRIAN / 16*0, 10*1 / 471 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE., 472 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE., 473 $ 3*.FALSE., 5*.TRUE., .FALSE. / 474 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE., 475 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE., 476 $ 9*.FALSE. / 477* .. 478* .. Executable Statements .. 479* 480* Check for errors 481* 482 INFO = 0 483* 484 BADNN = .FALSE. 485 NMAX = 1 486 DO 10 J = 1, NSIZES 487 NMAX = MAX( NMAX, NN( J ) ) 488 IF( NN( J ).LT.0 ) 489 $ BADNN = .TRUE. 490 10 CONTINUE 491* 492 IF( NSIZES.LT.0 ) THEN 493 INFO = -1 494 ELSE IF( BADNN ) THEN 495 INFO = -2 496 ELSE IF( NTYPES.LT.0 ) THEN 497 INFO = -3 498 ELSE IF( THRESH.LT.ZERO ) THEN 499 INFO = -6 500 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 501 INFO = -9 502 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN 503 INFO = -14 504 END IF 505* 506* Compute workspace 507* (Note: Comments in the code beginning "Workspace:" describe the 508* minimal amount of workspace needed at that point in the code, 509* as well as the preferred amount for good performance. 510* NB refers to the optimal block size for the immediately 511* following subroutine, as returned by ILAENV. 512* 513 MINWRK = 1 514 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 515 MINWRK = 3*NMAX*NMAX 516 NB = MAX( 1, ILAENV( 1, 'ZGEQRF', ' ', NMAX, NMAX, -1, -1 ), 517 $ ILAENV( 1, 'ZUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ), 518 $ ILAENV( 1, 'ZUNGQR', ' ', NMAX, NMAX, NMAX, -1 ) ) 519 MAXWRK = MAX( NMAX+NMAX*NB, 3*NMAX*NMAX ) 520 WORK( 1 ) = MAXWRK 521 END IF 522* 523 IF( LWORK.LT.MINWRK ) 524 $ INFO = -19 525* 526 IF( INFO.NE.0 ) THEN 527 CALL XERBLA( 'ZDRGES3', -INFO ) 528 RETURN 529 END IF 530* 531* Quick return if possible 532* 533 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 534 $ RETURN 535* 536 ULP = DLAMCH( 'Precision' ) 537 SAFMIN = DLAMCH( 'Safe minimum' ) 538 SAFMIN = SAFMIN / ULP 539 SAFMAX = ONE / SAFMIN 540 CALL DLABAD( SAFMIN, SAFMAX ) 541 ULPINV = ONE / ULP 542* 543* The values RMAGN(2:3) depend on N, see below. 544* 545 RMAGN( 0 ) = ZERO 546 RMAGN( 1 ) = ONE 547* 548* Loop over matrix sizes 549* 550 NTESTT = 0 551 NERRS = 0 552 NMATS = 0 553* 554 DO 190 JSIZE = 1, NSIZES 555 N = NN( JSIZE ) 556 N1 = MAX( 1, N ) 557 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 ) 558 RMAGN( 3 ) = SAFMIN*ULPINV*DBLE( N1 ) 559* 560 IF( NSIZES.NE.1 ) THEN 561 MTYPES = MIN( MAXTYP, NTYPES ) 562 ELSE 563 MTYPES = MIN( MAXTYP+1, NTYPES ) 564 END IF 565* 566* Loop over matrix types 567* 568 DO 180 JTYPE = 1, MTYPES 569 IF( .NOT.DOTYPE( JTYPE ) ) 570 $ GO TO 180 571 NMATS = NMATS + 1 572 NTEST = 0 573* 574* Save ISEED in case of an error. 575* 576 DO 20 J = 1, 4 577 IOLDSD( J ) = ISEED( J ) 578 20 CONTINUE 579* 580* Initialize RESULT 581* 582 DO 30 J = 1, 13 583 RESULT( J ) = ZERO 584 30 CONTINUE 585* 586* Generate test matrices A and B 587* 588* Description of control parameters: 589* 590* KZLASS: =1 means w/o rotation, =2 means w/ rotation, 591* =3 means random. 592* KATYPE: the "type" to be passed to ZLATM4 for computing A. 593* KAZERO: the pattern of zeros on the diagonal for A: 594* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 595* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 596* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 597* non-zero entries.) 598* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 599* =2: large, =3: small. 600* LASIGN: .TRUE. if the diagonal elements of A are to be 601* multiplied by a random magnitude 1 number. 602* KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. 603* KTRIAN: =0: don't fill in the upper triangle, =1: do. 604* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 605* RMAGN: used to implement KAMAGN and KBMAGN. 606* 607 IF( MTYPES.GT.MAXTYP ) 608 $ GO TO 110 609 IINFO = 0 610 IF( KCLASS( JTYPE ).LT.3 ) THEN 611* 612* Generate A (w/o rotation) 613* 614 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 615 IN = 2*( ( N-1 ) / 2 ) + 1 616 IF( IN.NE.N ) 617 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 618 ELSE 619 IN = N 620 END IF 621 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 622 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 623 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 624 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 625 $ ISEED, A, LDA ) 626 IADD = KADD( KAZERO( JTYPE ) ) 627 IF( IADD.GT.0 .AND. IADD.LE.N ) 628 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 629* 630* Generate B (w/o rotation) 631* 632 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 633 IN = 2*( ( N-1 ) / 2 ) + 1 634 IF( IN.NE.N ) 635 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 636 ELSE 637 IN = N 638 END IF 639 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 640 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 641 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 642 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 643 $ ISEED, B, LDA ) 644 IADD = KADD( KBZERO( JTYPE ) ) 645 IF( IADD.NE.0 .AND. IADD.LE.N ) 646 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 647* 648 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 649* 650* Include rotations 651* 652* Generate Q, Z as Householder transformations times 653* a diagonal matrix. 654* 655 DO 50 JC = 1, N - 1 656 DO 40 JR = JC, N 657 Q( JR, JC ) = ZLARND( 3, ISEED ) 658 Z( JR, JC ) = ZLARND( 3, ISEED ) 659 40 CONTINUE 660 CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 661 $ WORK( JC ) ) 662 WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) ) 663 Q( JC, JC ) = CONE 664 CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 665 $ WORK( N+JC ) ) 666 WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) ) 667 Z( JC, JC ) = CONE 668 50 CONTINUE 669 CTEMP = ZLARND( 3, ISEED ) 670 Q( N, N ) = CONE 671 WORK( N ) = CZERO 672 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 673 CTEMP = ZLARND( 3, ISEED ) 674 Z( N, N ) = CONE 675 WORK( 2*N ) = CZERO 676 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 677* 678* Apply the diagonal matrices 679* 680 DO 70 JC = 1, N 681 DO 60 JR = 1, N 682 A( JR, JC ) = WORK( 2*N+JR )* 683 $ DCONJG( WORK( 3*N+JC ) )* 684 $ A( JR, JC ) 685 B( JR, JC ) = WORK( 2*N+JR )* 686 $ DCONJG( WORK( 3*N+JC ) )* 687 $ B( JR, JC ) 688 60 CONTINUE 689 70 CONTINUE 690 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A, 691 $ LDA, WORK( 2*N+1 ), IINFO ) 692 IF( IINFO.NE.0 ) 693 $ GO TO 100 694 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 695 $ A, LDA, WORK( 2*N+1 ), IINFO ) 696 IF( IINFO.NE.0 ) 697 $ GO TO 100 698 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B, 699 $ LDA, WORK( 2*N+1 ), IINFO ) 700 IF( IINFO.NE.0 ) 701 $ GO TO 100 702 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 703 $ B, LDA, WORK( 2*N+1 ), IINFO ) 704 IF( IINFO.NE.0 ) 705 $ GO TO 100 706 END IF 707 ELSE 708* 709* Random matrices 710* 711 DO 90 JC = 1, N 712 DO 80 JR = 1, N 713 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 714 $ ZLARND( 4, ISEED ) 715 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 716 $ ZLARND( 4, ISEED ) 717 80 CONTINUE 718 90 CONTINUE 719 END IF 720* 721 100 CONTINUE 722* 723 IF( IINFO.NE.0 ) THEN 724 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 725 $ IOLDSD 726 INFO = ABS( IINFO ) 727 RETURN 728 END IF 729* 730 110 CONTINUE 731* 732 DO 120 I = 1, 13 733 RESULT( I ) = -ONE 734 120 CONTINUE 735* 736* Test with and without sorting of eigenvalues 737* 738 DO 150 ISORT = 0, 1 739 IF( ISORT.EQ.0 ) THEN 740 SORT = 'N' 741 RSUB = 0 742 ELSE 743 SORT = 'S' 744 RSUB = 5 745 END IF 746* 747* Call ZGGES3 to compute H, T, Q, Z, alpha, and beta. 748* 749 CALL ZLACPY( 'Full', N, N, A, LDA, S, LDA ) 750 CALL ZLACPY( 'Full', N, N, B, LDA, T, LDA ) 751 NTEST = 1 + RSUB + ISORT 752 RESULT( 1+RSUB+ISORT ) = ULPINV 753 CALL ZGGES3( 'V', 'V', SORT, ZLCTES, N, S, LDA, T, LDA, 754 $ SDIM, ALPHA, BETA, Q, LDQ, Z, LDQ, WORK, 755 $ LWORK, RWORK, BWORK, IINFO ) 756 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 757 RESULT( 1+RSUB+ISORT ) = ULPINV 758 WRITE( NOUNIT, FMT = 9999 )'ZGGES3', IINFO, N, JTYPE, 759 $ IOLDSD 760 INFO = ABS( IINFO ) 761 GO TO 160 762 END IF 763* 764 NTEST = 4 + RSUB 765* 766* Do tests 1--4 (or tests 7--9 when reordering ) 767* 768 IF( ISORT.EQ.0 ) THEN 769 CALL ZGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ, 770 $ WORK, RWORK, RESULT( 1 ) ) 771 CALL ZGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ, 772 $ WORK, RWORK, RESULT( 2 ) ) 773 ELSE 774 CALL ZGET54( N, A, LDA, B, LDA, S, LDA, T, LDA, Q, 775 $ LDQ, Z, LDQ, WORK, RESULT( 2+RSUB ) ) 776 END IF 777* 778 CALL ZGET51( 3, N, B, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK, 779 $ RWORK, RESULT( 3+RSUB ) ) 780 CALL ZGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK, 781 $ RWORK, RESULT( 4+RSUB ) ) 782* 783* Do test 5 and 6 (or Tests 10 and 11 when reordering): 784* check Schur form of A and compare eigenvalues with 785* diagonals. 786* 787 NTEST = 6 + RSUB 788 TEMP1 = ZERO 789* 790 DO 130 J = 1, N 791 ILABAD = .FALSE. 792 TEMP2 = ( ABS1( ALPHA( J )-S( J, J ) ) / 793 $ MAX( SAFMIN, ABS1( ALPHA( J ) ), ABS1( S( J, 794 $ J ) ) )+ABS1( BETA( J )-T( J, J ) ) / 795 $ MAX( SAFMIN, ABS1( BETA( J ) ), ABS1( T( J, 796 $ J ) ) ) ) / ULP 797* 798 IF( J.LT.N ) THEN 799 IF( S( J+1, J ).NE.ZERO ) THEN 800 ILABAD = .TRUE. 801 RESULT( 5+RSUB ) = ULPINV 802 END IF 803 END IF 804 IF( J.GT.1 ) THEN 805 IF( S( J, J-1 ).NE.ZERO ) THEN 806 ILABAD = .TRUE. 807 RESULT( 5+RSUB ) = ULPINV 808 END IF 809 END IF 810 TEMP1 = MAX( TEMP1, TEMP2 ) 811 IF( ILABAD ) THEN 812 WRITE( NOUNIT, FMT = 9998 )J, N, JTYPE, IOLDSD 813 END IF 814 130 CONTINUE 815 RESULT( 6+RSUB ) = TEMP1 816* 817 IF( ISORT.GE.1 ) THEN 818* 819* Do test 12 820* 821 NTEST = 12 822 RESULT( 12 ) = ZERO 823 KNTEIG = 0 824 DO 140 I = 1, N 825 IF( ZLCTES( ALPHA( I ), BETA( I ) ) ) 826 $ KNTEIG = KNTEIG + 1 827 140 CONTINUE 828 IF( SDIM.NE.KNTEIG ) 829 $ RESULT( 13 ) = ULPINV 830 END IF 831* 832 150 CONTINUE 833* 834* End of Loop -- Check for RESULT(j) > THRESH 835* 836 160 CONTINUE 837* 838 NTESTT = NTESTT + NTEST 839* 840* Print out tests which fail. 841* 842 DO 170 JR = 1, NTEST 843 IF( RESULT( JR ).GE.THRESH ) THEN 844* 845* If this is the first test to fail, 846* print a header to the data file. 847* 848 IF( NERRS.EQ.0 ) THEN 849 WRITE( NOUNIT, FMT = 9997 )'ZGS' 850* 851* Matrix types 852* 853 WRITE( NOUNIT, FMT = 9996 ) 854 WRITE( NOUNIT, FMT = 9995 ) 855 WRITE( NOUNIT, FMT = 9994 )'Unitary' 856* 857* Tests performed 858* 859 WRITE( NOUNIT, FMT = 9993 )'unitary', '''', 860 $ 'transpose', ( '''', J = 1, 8 ) 861* 862 END IF 863 NERRS = NERRS + 1 864 IF( RESULT( JR ).LT.10000.0D0 ) THEN 865 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 866 $ RESULT( JR ) 867 ELSE 868 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 869 $ RESULT( JR ) 870 END IF 871 END IF 872 170 CONTINUE 873* 874 180 CONTINUE 875 190 CONTINUE 876* 877* Summary 878* 879 CALL ALASVM( 'ZGS', NOUNIT, NERRS, NTESTT, 0 ) 880* 881 WORK( 1 ) = MAXWRK 882* 883 RETURN 884* 885 9999 FORMAT( ' ZDRGES3: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 886 $ I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' ) 887* 888 9998 FORMAT( ' ZDRGES3: S not in Schur form at eigenvalue ', I6, '.', 889 $ / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), 890 $ I5, ')' ) 891* 892 9997 FORMAT( / 1X, A3, ' -- Complex Generalized Schur from problem ', 893 $ 'driver' ) 894* 895 9996 FORMAT( ' Matrix types (see ZDRGES3 for details): ' ) 896* 897 9995 FORMAT( ' Special Matrices:', 23X, 898 $ '(J''=transposed Jordan block)', 899 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 900 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 901 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 902 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 903 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 904 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 905 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 906 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 907 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 908 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 909 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 910 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 911 $ '23=(small,large) 24=(small,small) 25=(large,large)', 912 $ / ' 26=random O(1) matrices.' ) 913* 914 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ', 915 $ 'Q and Z are ', A, ',', / 19X, 916 $ 'l and r are the appropriate left and right', / 19X, 917 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A, 918 $ ' means ', A, '.)', / ' Without ordering: ', 919 $ / ' 1 = | A - Q S Z', A, 920 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A, 921 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A, 922 $ ' | / ( n ulp ) 4 = | I - ZZ', A, 923 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S', 924 $ / ' 6 = difference between (alpha,beta)', 925 $ ' and diagonals of (S,T)', / ' With ordering: ', 926 $ / ' 7 = | (A,B) - Q (S,T) Z', A, ' | / ( |(A,B)| n ulp )', 927 $ / ' 8 = | I - QQ', A, 928 $ ' | / ( n ulp ) 9 = | I - ZZ', A, 929 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S', 930 $ / ' 11 = difference between (alpha,beta) and diagonals', 931 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ', 932 $ 'selected eigenvalues', / ) 933 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 934 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 935 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 936 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 ) 937* 938* End of ZDRGES3 939* 940 END 941