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