1*> \brief \b CCHKGG 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 CCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 12* TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1, 13* S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1, 14* ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK, 15* RWORK, LLWORK, RESULT, INFO ) 16* 17* .. Scalar Arguments .. 18* LOGICAL TSTDIF 19* INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES 20* REAL THRESH, THRSHN 21* .. 22* .. Array Arguments .. 23* LOGICAL DOTYPE( * ), LLWORK( * ) 24* INTEGER ISEED( 4 ), NN( * ) 25* REAL RESULT( 15 ), RWORK( * ) 26* COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ), 27* $ B( LDA, * ), BETA1( * ), BETA3( * ), 28* $ EVECTL( LDU, * ), EVECTR( LDU, * ), 29* $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ), 30* $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ), 31* $ T( LDA, * ), U( LDU, * ), V( LDU, * ), 32* $ WORK( * ), Z( LDU, * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> CCHKGG checks the nonsymmetric generalized eigenvalue problem 42*> routines. 43*> H H H 44*> CGGHRD factors A and B as U H V and U T V , where means conjugate 45*> transpose, H is hessenberg, T is triangular and U and V are unitary. 46*> 47*> H H 48*> CHGEQZ factors H and T as Q S Z and Q P Z , where P and S are upper 49*> triangular and Q and Z are unitary. It also computes the generalized 50*> eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)), where 51*> alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus, w(j) = alpha(j)/beta(j) 52*> is a root of the generalized eigenvalue problem 53*> 54*> det( A - w(j) B ) = 0 55*> 56*> and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent 57*> problem 58*> 59*> det( m(j) A - B ) = 0 60*> 61*> CTGEVC computes the matrix L of left eigenvectors and the matrix R 62*> of right eigenvectors for the matrix pair ( S, P ). In the 63*> description below, l and r are left and right eigenvectors 64*> corresponding to the generalized eigenvalues (alpha,beta). 65*> 66*> When CCHKGG is called, a number of matrix "sizes" ("n's") and a 67*> number of matrix "types" are specified. For each size ("n") 68*> and each type of matrix, one matrix will be generated and used 69*> to test the nonsymmetric eigenroutines. For each matrix, 13 70*> tests will be performed. The first twelve "test ratios" should be 71*> small -- O(1). They will be compared with the threshold THRESH: 72*> 73*> H 74*> (1) | A - U H V | / ( |A| n ulp ) 75*> 76*> H 77*> (2) | B - U T V | / ( |B| n ulp ) 78*> 79*> H 80*> (3) | I - UU | / ( n ulp ) 81*> 82*> H 83*> (4) | I - VV | / ( n ulp ) 84*> 85*> H 86*> (5) | H - Q S Z | / ( |H| n ulp ) 87*> 88*> H 89*> (6) | T - Q P Z | / ( |T| n ulp ) 90*> 91*> H 92*> (7) | I - QQ | / ( n ulp ) 93*> 94*> H 95*> (8) | I - ZZ | / ( n ulp ) 96*> 97*> (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of 98*> H 99*> | (beta A - alpha B) l | / ( ulp max( |beta A|, |alpha B| ) ) 100*> 101*> (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of 102*> H 103*> | (beta H - alpha T) l' | / ( ulp max( |beta H|, |alpha T| ) ) 104*> 105*> where the eigenvectors l' are the result of passing Q to 106*> STGEVC and back transforming (JOB='B'). 107*> 108*> (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of 109*> 110*> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) ) 111*> 112*> (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of 113*> 114*> | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) ) 115*> 116*> where the eigenvectors r' are the result of passing Z to 117*> STGEVC and back transforming (JOB='B'). 118*> 119*> The last three test ratios will usually be small, but there is no 120*> mathematical requirement that they be so. They are therefore 121*> compared with THRESH only if TSTDIF is .TRUE. 122*> 123*> (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp ) 124*> 125*> (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp ) 126*> 127*> (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| , 128*> |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp 129*> 130*> In addition, the normalization of L and R are checked, and compared 131*> with the threshold THRSHN. 132*> 133*> Test Matrices 134*> ---- -------- 135*> 136*> The sizes of the test matrices are specified by an array 137*> NN(1:NSIZES); the value of each element NN(j) specifies one size. 138*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if 139*> DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 140*> Currently, the list of possible types is: 141*> 142*> (1) ( 0, 0 ) (a pair of zero matrices) 143*> 144*> (2) ( I, 0 ) (an identity and a zero matrix) 145*> 146*> (3) ( 0, I ) (an identity and a zero matrix) 147*> 148*> (4) ( I, I ) (a pair of identity matrices) 149*> 150*> t t 151*> (5) ( J , J ) (a pair of transposed Jordan blocks) 152*> 153*> t ( I 0 ) 154*> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t ) 155*> ( 0 I ) ( 0 J ) 156*> and I is a k x k identity and J a (k+1)x(k+1) 157*> Jordan block; k=(N-1)/2 158*> 159*> (7) ( D, I ) where D is P*D1, P is a random unitary diagonal 160*> matrix (i.e., with random magnitude 1 entries 161*> on the diagonal), and D1=diag( 0, 1,..., N-1 ) 162*> (i.e., a diagonal matrix with D1(1,1)=0, 163*> D1(2,2)=1, ..., D1(N,N)=N-1.) 164*> (8) ( I, D ) 165*> 166*> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big 167*> 168*> (10) ( small*D, big*I ) 169*> 170*> (11) ( big*I, small*D ) 171*> 172*> (12) ( small*I, big*D ) 173*> 174*> (13) ( big*D, big*I ) 175*> 176*> (14) ( small*D, small*I ) 177*> 178*> (15) ( D1, D2 ) where D1=P*diag( 0, 0, 1, ..., N-3, 0 ) and 179*> D2=Q*diag( 0, N-3, N-4,..., 1, 0, 0 ), and 180*> P and Q are random unitary diagonal matrices. 181*> t t 182*> (16) U ( J , J ) V where U and V are random unitary matrices. 183*> 184*> (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices 185*> with random O(1) entries above the diagonal 186*> and diagonal entries diag(T1) = 187*> P*( 0, 0, 1, ..., N-3, 0 ) and diag(T2) = 188*> Q*( 0, N-3, N-4,..., 1, 0, 0 ) 189*> 190*> (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 ) 191*> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 ) 192*> s = machine precision. 193*> 194*> (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 ) 195*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 ) 196*> 197*> N-5 198*> (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 ) 199*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 200*> 201*> (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 ) 202*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 203*> where r1,..., r(N-4) are random. 204*> 205*> (22) U ( big*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 206*> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 207*> 208*> (23) U ( small*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 209*> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 210*> 211*> (24) U ( small*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 212*> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 213*> 214*> (25) U ( big*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 215*> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 216*> 217*> (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular 218*> matrices. 219*> \endverbatim 220* 221* Arguments: 222* ========== 223* 224*> \param[in] NSIZES 225*> \verbatim 226*> NSIZES is INTEGER 227*> The number of sizes of matrices to use. If it is zero, 228*> CCHKGG does nothing. It must be at least zero. 229*> \endverbatim 230*> 231*> \param[in] NN 232*> \verbatim 233*> NN is INTEGER array, dimension (NSIZES) 234*> An array containing the sizes to be used for the matrices. 235*> Zero values will be skipped. The values must be at least 236*> zero. 237*> \endverbatim 238*> 239*> \param[in] NTYPES 240*> \verbatim 241*> NTYPES is INTEGER 242*> The number of elements in DOTYPE. If it is zero, CCHKGG 243*> does nothing. It must be at least zero. If it is MAXTYP+1 244*> and NSIZES is 1, then an additional type, MAXTYP+1 is 245*> defined, which is to use whatever matrix is in A. This 246*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 247*> DOTYPE(MAXTYP+1) is .TRUE. . 248*> \endverbatim 249*> 250*> \param[in] DOTYPE 251*> \verbatim 252*> DOTYPE is LOGICAL array, dimension (NTYPES) 253*> If DOTYPE(j) is .TRUE., then for each size in NN a 254*> matrix of that size and of type j will be generated. 255*> If NTYPES is smaller than the maximum number of types 256*> defined (PARAMETER MAXTYP), then types NTYPES+1 through 257*> MAXTYP will not be generated. If NTYPES is larger 258*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 259*> will be ignored. 260*> \endverbatim 261*> 262*> \param[in,out] ISEED 263*> \verbatim 264*> ISEED is INTEGER array, dimension (4) 265*> On entry ISEED specifies the seed of the random number 266*> generator. The array elements should be between 0 and 4095; 267*> if not they will be reduced mod 4096. Also, ISEED(4) must 268*> be odd. The random number generator uses a linear 269*> congruential sequence limited to small integers, and so 270*> should produce machine independent random numbers. The 271*> values of ISEED are changed on exit, and can be used in the 272*> next call to CCHKGG to continue the same random number 273*> sequence. 274*> \endverbatim 275*> 276*> \param[in] THRESH 277*> \verbatim 278*> THRESH is REAL 279*> A test will count as "failed" if the "error", computed as 280*> described above, exceeds THRESH. Note that the error 281*> is scaled to be O(1), so THRESH should be a reasonably 282*> small multiple of 1, e.g., 10 or 100. In particular, 283*> it should not depend on the precision (single vs. double) 284*> or the size of the matrix. It must be at least zero. 285*> \endverbatim 286*> 287*> \param[in] TSTDIF 288*> \verbatim 289*> TSTDIF is LOGICAL 290*> Specifies whether test ratios 13-15 will be computed and 291*> compared with THRESH. 292*> = .FALSE.: Only test ratios 1-12 will be computed and tested. 293*> Ratios 13-15 will be set to zero. 294*> = .TRUE.: All the test ratios 1-15 will be computed and 295*> tested. 296*> \endverbatim 297*> 298*> \param[in] THRSHN 299*> \verbatim 300*> THRSHN is REAL 301*> Threshold for reporting eigenvector normalization error. 302*> If the normalization of any eigenvector differs from 1 by 303*> more than THRSHN*ulp, then a special error message will be 304*> printed. (This is handled separately from the other tests, 305*> since only a compiler or programming error should cause an 306*> error message, at least if THRSHN is at least 5--10.) 307*> \endverbatim 308*> 309*> \param[in] NOUNIT 310*> \verbatim 311*> NOUNIT is INTEGER 312*> The FORTRAN unit number for printing out error messages 313*> (e.g., if a routine returns IINFO not equal to 0.) 314*> \endverbatim 315*> 316*> \param[in,out] A 317*> \verbatim 318*> A is COMPLEX array, dimension (LDA, max(NN)) 319*> Used to hold the original A matrix. Used as input only 320*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 321*> DOTYPE(MAXTYP+1)=.TRUE. 322*> \endverbatim 323*> 324*> \param[in] LDA 325*> \verbatim 326*> LDA is INTEGER 327*> The leading dimension of A, B, H, T, S1, P1, S2, and P2. 328*> It must be at least 1 and at least max( NN ). 329*> \endverbatim 330*> 331*> \param[in,out] B 332*> \verbatim 333*> B is COMPLEX array, dimension (LDA, max(NN)) 334*> Used to hold the original B matrix. Used as input only 335*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 336*> DOTYPE(MAXTYP+1)=.TRUE. 337*> \endverbatim 338*> 339*> \param[out] H 340*> \verbatim 341*> H is COMPLEX array, dimension (LDA, max(NN)) 342*> The upper Hessenberg matrix computed from A by CGGHRD. 343*> \endverbatim 344*> 345*> \param[out] T 346*> \verbatim 347*> T is COMPLEX array, dimension (LDA, max(NN)) 348*> The upper triangular matrix computed from B by CGGHRD. 349*> \endverbatim 350*> 351*> \param[out] S1 352*> \verbatim 353*> S1 is COMPLEX array, dimension (LDA, max(NN)) 354*> The Schur (upper triangular) matrix computed from H by CHGEQZ 355*> when Q and Z are also computed. 356*> \endverbatim 357*> 358*> \param[out] S2 359*> \verbatim 360*> S2 is COMPLEX array, dimension (LDA, max(NN)) 361*> The Schur (upper triangular) matrix computed from H by CHGEQZ 362*> when Q and Z are not computed. 363*> \endverbatim 364*> 365*> \param[out] P1 366*> \verbatim 367*> P1 is COMPLEX array, dimension (LDA, max(NN)) 368*> The upper triangular matrix computed from T by CHGEQZ 369*> when Q and Z are also computed. 370*> \endverbatim 371*> 372*> \param[out] P2 373*> \verbatim 374*> P2 is COMPLEX array, dimension (LDA, max(NN)) 375*> The upper triangular matrix computed from T by CHGEQZ 376*> when Q and Z are not computed. 377*> \endverbatim 378*> 379*> \param[out] U 380*> \verbatim 381*> U is COMPLEX array, dimension (LDU, max(NN)) 382*> The (left) unitary matrix computed by CGGHRD. 383*> \endverbatim 384*> 385*> \param[in] LDU 386*> \verbatim 387*> LDU is INTEGER 388*> The leading dimension of U, V, Q, Z, EVECTL, and EVECTR. It 389*> must be at least 1 and at least max( NN ). 390*> \endverbatim 391*> 392*> \param[out] V 393*> \verbatim 394*> V is COMPLEX array, dimension (LDU, max(NN)) 395*> The (right) unitary matrix computed by CGGHRD. 396*> \endverbatim 397*> 398*> \param[out] Q 399*> \verbatim 400*> Q is COMPLEX array, dimension (LDU, max(NN)) 401*> The (left) unitary matrix computed by CHGEQZ. 402*> \endverbatim 403*> 404*> \param[out] Z 405*> \verbatim 406*> Z is COMPLEX array, dimension (LDU, max(NN)) 407*> The (left) unitary matrix computed by CHGEQZ. 408*> \endverbatim 409*> 410*> \param[out] ALPHA1 411*> \verbatim 412*> ALPHA1 is COMPLEX array, dimension (max(NN)) 413*> \endverbatim 414*> 415*> \param[out] BETA1 416*> \verbatim 417*> BETA1 is COMPLEX array, dimension (max(NN)) 418*> The generalized eigenvalues of (A,B) computed by CHGEQZ 419*> when Q, Z, and the full Schur matrices are computed. 420*> \endverbatim 421*> 422*> \param[out] ALPHA3 423*> \verbatim 424*> ALPHA3 is COMPLEX array, dimension (max(NN)) 425*> \endverbatim 426*> 427*> \param[out] BETA3 428*> \verbatim 429*> BETA3 is COMPLEX array, dimension (max(NN)) 430*> The generalized eigenvalues of (A,B) computed by CHGEQZ 431*> when neither Q, Z, nor the Schur matrices are computed. 432*> \endverbatim 433*> 434*> \param[out] EVECTL 435*> \verbatim 436*> EVECTL is COMPLEX array, dimension (LDU, max(NN)) 437*> The (lower triangular) left eigenvector matrix for the 438*> matrices in S1 and P1. 439*> \endverbatim 440*> 441*> \param[out] EVECTR 442*> \verbatim 443*> EVECTR is COMPLEX array, dimension (LDU, max(NN)) 444*> The (upper triangular) right eigenvector matrix for the 445*> matrices in S1 and P1. 446*> \endverbatim 447*> 448*> \param[out] WORK 449*> \verbatim 450*> WORK is COMPLEX array, dimension (LWORK) 451*> \endverbatim 452*> 453*> \param[in] LWORK 454*> \verbatim 455*> LWORK is INTEGER 456*> The number of entries in WORK. This must be at least 457*> max( 4*N, 2 * N**2, 1 ), for all N=NN(j). 458*> \endverbatim 459*> 460*> \param[out] RWORK 461*> \verbatim 462*> RWORK is REAL array, dimension (2*max(NN)) 463*> \endverbatim 464*> 465*> \param[out] LLWORK 466*> \verbatim 467*> LLWORK is LOGICAL array, dimension (max(NN)) 468*> \endverbatim 469*> 470*> \param[out] RESULT 471*> \verbatim 472*> RESULT is REAL array, dimension (15) 473*> The values computed by the tests described above. 474*> The values are currently limited to 1/ulp, to avoid 475*> overflow. 476*> \endverbatim 477*> 478*> \param[out] INFO 479*> \verbatim 480*> INFO is INTEGER 481*> = 0: successful exit. 482*> < 0: if INFO = -i, the i-th argument had an illegal value. 483*> > 0: A routine returned an error code. INFO is the 484*> absolute value of the INFO value returned. 485*> \endverbatim 486* 487* Authors: 488* ======== 489* 490*> \author Univ. of Tennessee 491*> \author Univ. of California Berkeley 492*> \author Univ. of Colorado Denver 493*> \author NAG Ltd. 494* 495*> \ingroup complex_eig 496* 497* ===================================================================== 498 SUBROUTINE CCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 499 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1, 500 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1, 501 $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK, 502 $ RWORK, LLWORK, RESULT, INFO ) 503* 504* -- LAPACK test routine -- 505* -- LAPACK is a software package provided by Univ. of Tennessee, -- 506* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 507* 508* .. Scalar Arguments .. 509 LOGICAL TSTDIF 510 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES 511 REAL THRESH, THRSHN 512* .. 513* .. Array Arguments .. 514 LOGICAL DOTYPE( * ), LLWORK( * ) 515 INTEGER ISEED( 4 ), NN( * ) 516 REAL RESULT( 15 ), RWORK( * ) 517 COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ), 518 $ B( LDA, * ), BETA1( * ), BETA3( * ), 519 $ EVECTL( LDU, * ), EVECTR( LDU, * ), 520 $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ), 521 $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ), 522 $ T( LDA, * ), U( LDU, * ), V( LDU, * ), 523 $ WORK( * ), Z( LDU, * ) 524* .. 525* 526* ===================================================================== 527* 528* .. Parameters .. 529 REAL ZERO, ONE 530 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 531 COMPLEX CZERO, CONE 532 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 533 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 534 INTEGER MAXTYP 535 PARAMETER ( MAXTYP = 26 ) 536* .. 537* .. Local Scalars .. 538 LOGICAL BADNN 539 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE, 540 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX, 541 $ NTEST, NTESTT 542 REAL ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2, 543 $ ULP, ULPINV 544 COMPLEX CTEMP 545* .. 546* .. Local Arrays .. 547 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP ) 548 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 549 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 550 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 551 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 552 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 553 REAL DUMMA( 4 ), RMAGN( 0: 3 ) 554 COMPLEX CDUMMA( 4 ) 555* .. 556* .. External Functions .. 557 REAL CLANGE, SLAMCH 558 COMPLEX CLARND 559 EXTERNAL CLANGE, SLAMCH, CLARND 560* .. 561* .. External Subroutines .. 562 EXTERNAL CGEQR2, CGET51, CGET52, CGGHRD, CHGEQZ, CLACPY, 563 $ CLARFG, CLASET, CLATM4, CTGEVC, CUNM2R, SLABAD, 564 $ SLASUM, XERBLA 565* .. 566* .. Intrinsic Functions .. 567 INTRINSIC ABS, CONJG, MAX, MIN, REAL, SIGN 568* .. 569* .. Data statements .. 570 DATA KCLASS / 15*1, 10*2, 1*3 / 571 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 572 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 573 DATA KADD / 0, 0, 0, 0, 3, 2 / 574 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 575 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 576 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 577 $ 1, 1, -4, 2, -4, 8*8, 0 / 578 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 579 $ 4*5, 4*3, 1 / 580 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 581 $ 4*6, 4*4, 1 / 582 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 583 $ 2, 1 / 584 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 585 $ 2, 1 / 586 DATA KTRIAN / 16*0, 10*1 / 587 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE., 588 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE., 589 $ 3*.FALSE., 5*.TRUE., .FALSE. / 590 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE., 591 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE., 592 $ 9*.FALSE. / 593* .. 594* .. Executable Statements .. 595* 596* Check for errors 597* 598 INFO = 0 599* 600 BADNN = .FALSE. 601 NMAX = 1 602 DO 10 J = 1, NSIZES 603 NMAX = MAX( NMAX, NN( J ) ) 604 IF( NN( J ).LT.0 ) 605 $ BADNN = .TRUE. 606 10 CONTINUE 607* 608 LWKOPT = MAX( 2*NMAX*NMAX, 4*NMAX, 1 ) 609* 610* Check for errors 611* 612 IF( NSIZES.LT.0 ) THEN 613 INFO = -1 614 ELSE IF( BADNN ) THEN 615 INFO = -2 616 ELSE IF( NTYPES.LT.0 ) THEN 617 INFO = -3 618 ELSE IF( THRESH.LT.ZERO ) THEN 619 INFO = -6 620 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 621 INFO = -10 622 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN 623 INFO = -19 624 ELSE IF( LWKOPT.GT.LWORK ) THEN 625 INFO = -30 626 END IF 627* 628 IF( INFO.NE.0 ) THEN 629 CALL XERBLA( 'CCHKGG', -INFO ) 630 RETURN 631 END IF 632* 633* Quick return if possible 634* 635 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 636 $ RETURN 637* 638 SAFMIN = SLAMCH( 'Safe minimum' ) 639 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 640 SAFMIN = SAFMIN / ULP 641 SAFMAX = ONE / SAFMIN 642 CALL SLABAD( SAFMIN, SAFMAX ) 643 ULPINV = ONE / ULP 644* 645* The values RMAGN(2:3) depend on N, see below. 646* 647 RMAGN( 0 ) = ZERO 648 RMAGN( 1 ) = ONE 649* 650* Loop over sizes, types 651* 652 NTESTT = 0 653 NERRS = 0 654 NMATS = 0 655* 656 DO 240 JSIZE = 1, NSIZES 657 N = NN( JSIZE ) 658 N1 = MAX( 1, N ) 659 RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 ) 660 RMAGN( 3 ) = SAFMIN*ULPINV*N1 661* 662 IF( NSIZES.NE.1 ) THEN 663 MTYPES = MIN( MAXTYP, NTYPES ) 664 ELSE 665 MTYPES = MIN( MAXTYP+1, NTYPES ) 666 END IF 667* 668 DO 230 JTYPE = 1, MTYPES 669 IF( .NOT.DOTYPE( JTYPE ) ) 670 $ GO TO 230 671 NMATS = NMATS + 1 672 NTEST = 0 673* 674* Save ISEED in case of an error. 675* 676 DO 20 J = 1, 4 677 IOLDSD( J ) = ISEED( J ) 678 20 CONTINUE 679* 680* Initialize RESULT 681* 682 DO 30 J = 1, 15 683 RESULT( J ) = ZERO 684 30 CONTINUE 685* 686* Compute A and B 687* 688* Description of control parameters: 689* 690* KCLASS: =1 means w/o rotation, =2 means w/ rotation, 691* =3 means random. 692* KATYPE: the "type" to be passed to CLATM4 for computing A. 693* KAZERO: the pattern of zeros on the diagonal for A: 694* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 695* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 696* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 697* non-zero entries.) 698* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 699* =2: large, =3: small. 700* LASIGN: .TRUE. if the diagonal elements of A are to be 701* multiplied by a random magnitude 1 number. 702* KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. 703* KTRIAN: =0: don't fill in the upper triangle, =1: do. 704* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 705* RMAGN: used to implement KAMAGN and KBMAGN. 706* 707 IF( MTYPES.GT.MAXTYP ) 708 $ GO TO 110 709 IINFO = 0 710 IF( KCLASS( JTYPE ).LT.3 ) THEN 711* 712* Generate A (w/o rotation) 713* 714 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 715 IN = 2*( ( N-1 ) / 2 ) + 1 716 IF( IN.NE.N ) 717 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 718 ELSE 719 IN = N 720 END IF 721 CALL CLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 722 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 723 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 724 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4, 725 $ ISEED, A, LDA ) 726 IADD = KADD( KAZERO( JTYPE ) ) 727 IF( IADD.GT.0 .AND. IADD.LE.N ) 728 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 729* 730* Generate B (w/o rotation) 731* 732 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 733 IN = 2*( ( N-1 ) / 2 ) + 1 734 IF( IN.NE.N ) 735 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 736 ELSE 737 IN = N 738 END IF 739 CALL CLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 740 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 741 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 742 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4, 743 $ ISEED, B, LDA ) 744 IADD = KADD( KBZERO( JTYPE ) ) 745 IF( IADD.NE.0 ) 746 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 747* 748 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 749* 750* Include rotations 751* 752* Generate U, V as Householder transformations times a 753* diagonal matrix. (Note that CLARFG makes U(j,j) and 754* V(j,j) real.) 755* 756 DO 50 JC = 1, N - 1 757 DO 40 JR = JC, N 758 U( JR, JC ) = CLARND( 3, ISEED ) 759 V( JR, JC ) = CLARND( 3, ISEED ) 760 40 CONTINUE 761 CALL CLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1, 762 $ WORK( JC ) ) 763 WORK( 2*N+JC ) = SIGN( ONE, REAL( U( JC, JC ) ) ) 764 U( JC, JC ) = CONE 765 CALL CLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1, 766 $ WORK( N+JC ) ) 767 WORK( 3*N+JC ) = SIGN( ONE, REAL( V( JC, JC ) ) ) 768 V( JC, JC ) = CONE 769 50 CONTINUE 770 CTEMP = CLARND( 3, ISEED ) 771 U( N, N ) = CONE 772 WORK( N ) = CZERO 773 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 774 CTEMP = CLARND( 3, ISEED ) 775 V( N, N ) = CONE 776 WORK( 2*N ) = CZERO 777 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 778* 779* Apply the diagonal matrices 780* 781 DO 70 JC = 1, N 782 DO 60 JR = 1, N 783 A( JR, JC ) = WORK( 2*N+JR )* 784 $ CONJG( WORK( 3*N+JC ) )* 785 $ A( JR, JC ) 786 B( JR, JC ) = WORK( 2*N+JR )* 787 $ CONJG( WORK( 3*N+JC ) )* 788 $ B( JR, JC ) 789 60 CONTINUE 790 70 CONTINUE 791 CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A, 792 $ LDA, WORK( 2*N+1 ), IINFO ) 793 IF( IINFO.NE.0 ) 794 $ GO TO 100 795 CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ), 796 $ A, LDA, WORK( 2*N+1 ), IINFO ) 797 IF( IINFO.NE.0 ) 798 $ GO TO 100 799 CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B, 800 $ LDA, WORK( 2*N+1 ), IINFO ) 801 IF( IINFO.NE.0 ) 802 $ GO TO 100 803 CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ), 804 $ B, LDA, WORK( 2*N+1 ), IINFO ) 805 IF( IINFO.NE.0 ) 806 $ GO TO 100 807 END IF 808 ELSE 809* 810* Random matrices 811* 812 DO 90 JC = 1, N 813 DO 80 JR = 1, N 814 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 815 $ CLARND( 4, ISEED ) 816 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 817 $ CLARND( 4, ISEED ) 818 80 CONTINUE 819 90 CONTINUE 820 END IF 821* 822 ANORM = CLANGE( '1', N, N, A, LDA, RWORK ) 823 BNORM = CLANGE( '1', N, N, B, LDA, RWORK ) 824* 825 100 CONTINUE 826* 827 IF( IINFO.NE.0 ) THEN 828 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 829 $ IOLDSD 830 INFO = ABS( IINFO ) 831 RETURN 832 END IF 833* 834 110 CONTINUE 835* 836* Call CGEQR2, CUNM2R, and CGGHRD to compute H, T, U, and V 837* 838 CALL CLACPY( ' ', N, N, A, LDA, H, LDA ) 839 CALL CLACPY( ' ', N, N, B, LDA, T, LDA ) 840 NTEST = 1 841 RESULT( 1 ) = ULPINV 842* 843 CALL CGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO ) 844 IF( IINFO.NE.0 ) THEN 845 WRITE( NOUNIT, FMT = 9999 )'CGEQR2', IINFO, N, JTYPE, 846 $ IOLDSD 847 INFO = ABS( IINFO ) 848 GO TO 210 849 END IF 850* 851 CALL CUNM2R( 'L', 'C', N, N, N, T, LDA, WORK, H, LDA, 852 $ WORK( N+1 ), IINFO ) 853 IF( IINFO.NE.0 ) THEN 854 WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE, 855 $ IOLDSD 856 INFO = ABS( IINFO ) 857 GO TO 210 858 END IF 859* 860 CALL CLASET( 'Full', N, N, CZERO, CONE, U, LDU ) 861 CALL CUNM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU, 862 $ WORK( N+1 ), IINFO ) 863 IF( IINFO.NE.0 ) THEN 864 WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE, 865 $ IOLDSD 866 INFO = ABS( IINFO ) 867 GO TO 210 868 END IF 869* 870 CALL CGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V, 871 $ LDU, IINFO ) 872 IF( IINFO.NE.0 ) THEN 873 WRITE( NOUNIT, FMT = 9999 )'CGGHRD', IINFO, N, JTYPE, 874 $ IOLDSD 875 INFO = ABS( IINFO ) 876 GO TO 210 877 END IF 878 NTEST = 4 879* 880* Do tests 1--4 881* 882 CALL CGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK, 883 $ RWORK, RESULT( 1 ) ) 884 CALL CGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK, 885 $ RWORK, RESULT( 2 ) ) 886 CALL CGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK, 887 $ RWORK, RESULT( 3 ) ) 888 CALL CGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK, 889 $ RWORK, RESULT( 4 ) ) 890* 891* Call CHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests. 892* 893* Compute T1 and UZ 894* 895* Eigenvalues only 896* 897 CALL CLACPY( ' ', N, N, H, LDA, S2, LDA ) 898 CALL CLACPY( ' ', N, N, T, LDA, P2, LDA ) 899 NTEST = 5 900 RESULT( 5 ) = ULPINV 901* 902 CALL CHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA, 903 $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK, 904 $ RWORK, IINFO ) 905 IF( IINFO.NE.0 ) THEN 906 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(E)', IINFO, N, JTYPE, 907 $ IOLDSD 908 INFO = ABS( IINFO ) 909 GO TO 210 910 END IF 911* 912* Eigenvalues and Full Schur Form 913* 914 CALL CLACPY( ' ', N, N, H, LDA, S2, LDA ) 915 CALL CLACPY( ' ', N, N, T, LDA, P2, LDA ) 916* 917 CALL CHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA, 918 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK, 919 $ RWORK, IINFO ) 920 IF( IINFO.NE.0 ) THEN 921 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(S)', IINFO, N, JTYPE, 922 $ IOLDSD 923 INFO = ABS( IINFO ) 924 GO TO 210 925 END IF 926* 927* Eigenvalues, Schur Form, and Schur Vectors 928* 929 CALL CLACPY( ' ', N, N, H, LDA, S1, LDA ) 930 CALL CLACPY( ' ', N, N, T, LDA, P1, LDA ) 931* 932 CALL CHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA, 933 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK, 934 $ RWORK, IINFO ) 935 IF( IINFO.NE.0 ) THEN 936 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(V)', IINFO, N, JTYPE, 937 $ IOLDSD 938 INFO = ABS( IINFO ) 939 GO TO 210 940 END IF 941* 942 NTEST = 8 943* 944* Do Tests 5--8 945* 946 CALL CGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK, 947 $ RWORK, RESULT( 5 ) ) 948 CALL CGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK, 949 $ RWORK, RESULT( 6 ) ) 950 CALL CGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK, 951 $ RWORK, RESULT( 7 ) ) 952 CALL CGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK, 953 $ RWORK, RESULT( 8 ) ) 954* 955* Compute the Left and Right Eigenvectors of (S1,P1) 956* 957* 9: Compute the left eigenvector Matrix without 958* back transforming: 959* 960 NTEST = 9 961 RESULT( 9 ) = ULPINV 962* 963* To test "SELECT" option, compute half of the eigenvectors 964* in one call, and half in another 965* 966 I1 = N / 2 967 DO 120 J = 1, I1 968 LLWORK( J ) = .TRUE. 969 120 CONTINUE 970 DO 130 J = I1 + 1, N 971 LLWORK( J ) = .FALSE. 972 130 CONTINUE 973* 974 CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL, 975 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 976 IF( IINFO.NE.0 ) THEN 977 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S1)', IINFO, N, 978 $ JTYPE, IOLDSD 979 INFO = ABS( IINFO ) 980 GO TO 210 981 END IF 982* 983 I1 = IN 984 DO 140 J = 1, I1 985 LLWORK( J ) = .FALSE. 986 140 CONTINUE 987 DO 150 J = I1 + 1, N 988 LLWORK( J ) = .TRUE. 989 150 CONTINUE 990* 991 CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, 992 $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN, 993 $ WORK, RWORK, IINFO ) 994 IF( IINFO.NE.0 ) THEN 995 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S2)', IINFO, N, 996 $ JTYPE, IOLDSD 997 INFO = ABS( IINFO ) 998 GO TO 210 999 END IF 1000* 1001 CALL CGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU, 1002 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 1003 RESULT( 9 ) = DUMMA( 1 ) 1004 IF( DUMMA( 2 ).GT.THRSHN ) THEN 1005 WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=S)', 1006 $ DUMMA( 2 ), N, JTYPE, IOLDSD 1007 END IF 1008* 1009* 10: Compute the left eigenvector Matrix with 1010* back transforming: 1011* 1012 NTEST = 10 1013 RESULT( 10 ) = ULPINV 1014 CALL CLACPY( 'F', N, N, Q, LDU, EVECTL, LDU ) 1015 CALL CTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL, 1016 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 1017 IF( IINFO.NE.0 ) THEN 1018 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,B)', IINFO, N, 1019 $ JTYPE, IOLDSD 1020 INFO = ABS( IINFO ) 1021 GO TO 210 1022 END IF 1023* 1024 CALL CGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1, 1025 $ BETA1, WORK, RWORK, DUMMA( 1 ) ) 1026 RESULT( 10 ) = DUMMA( 1 ) 1027 IF( DUMMA( 2 ).GT.THRSHN ) THEN 1028 WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=B)', 1029 $ DUMMA( 2 ), N, JTYPE, IOLDSD 1030 END IF 1031* 1032* 11: Compute the right eigenvector Matrix without 1033* back transforming: 1034* 1035 NTEST = 11 1036 RESULT( 11 ) = ULPINV 1037* 1038* To test "SELECT" option, compute half of the eigenvectors 1039* in one call, and half in another 1040* 1041 I1 = N / 2 1042 DO 160 J = 1, I1 1043 LLWORK( J ) = .TRUE. 1044 160 CONTINUE 1045 DO 170 J = I1 + 1, N 1046 LLWORK( J ) = .FALSE. 1047 170 CONTINUE 1048* 1049 CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 1050 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 1051 IF( IINFO.NE.0 ) THEN 1052 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S1)', IINFO, N, 1053 $ JTYPE, IOLDSD 1054 INFO = ABS( IINFO ) 1055 GO TO 210 1056 END IF 1057* 1058 I1 = IN 1059 DO 180 J = 1, I1 1060 LLWORK( J ) = .FALSE. 1061 180 CONTINUE 1062 DO 190 J = I1 + 1, N 1063 LLWORK( J ) = .TRUE. 1064 190 CONTINUE 1065* 1066 CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 1067 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK, 1068 $ RWORK, IINFO ) 1069 IF( IINFO.NE.0 ) THEN 1070 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S2)', IINFO, N, 1071 $ JTYPE, IOLDSD 1072 INFO = ABS( IINFO ) 1073 GO TO 210 1074 END IF 1075* 1076 CALL CGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU, 1077 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 1078 RESULT( 11 ) = DUMMA( 1 ) 1079 IF( DUMMA( 2 ).GT.THRESH ) THEN 1080 WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=S)', 1081 $ DUMMA( 2 ), N, JTYPE, IOLDSD 1082 END IF 1083* 1084* 12: Compute the right eigenvector Matrix with 1085* back transforming: 1086* 1087 NTEST = 12 1088 RESULT( 12 ) = ULPINV 1089 CALL CLACPY( 'F', N, N, Z, LDU, EVECTR, LDU ) 1090 CALL CTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 1091 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 1092 IF( IINFO.NE.0 ) THEN 1093 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,B)', IINFO, N, 1094 $ JTYPE, IOLDSD 1095 INFO = ABS( IINFO ) 1096 GO TO 210 1097 END IF 1098* 1099 CALL CGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU, 1100 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 1101 RESULT( 12 ) = DUMMA( 1 ) 1102 IF( DUMMA( 2 ).GT.THRESH ) THEN 1103 WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=B)', 1104 $ DUMMA( 2 ), N, JTYPE, IOLDSD 1105 END IF 1106* 1107* Tests 13--15 are done only on request 1108* 1109 IF( TSTDIF ) THEN 1110* 1111* Do Tests 13--14 1112* 1113 CALL CGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU, 1114 $ WORK, RWORK, RESULT( 13 ) ) 1115 CALL CGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU, 1116 $ WORK, RWORK, RESULT( 14 ) ) 1117* 1118* Do Test 15 1119* 1120 TEMP1 = ZERO 1121 TEMP2 = ZERO 1122 DO 200 J = 1, N 1123 TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) ) 1124 TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) ) 1125 200 CONTINUE 1126* 1127 TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) ) 1128 TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) ) 1129 RESULT( 15 ) = MAX( TEMP1, TEMP2 ) 1130 NTEST = 15 1131 ELSE 1132 RESULT( 13 ) = ZERO 1133 RESULT( 14 ) = ZERO 1134 RESULT( 15 ) = ZERO 1135 NTEST = 12 1136 END IF 1137* 1138* End of Loop -- Check for RESULT(j) > THRESH 1139* 1140 210 CONTINUE 1141* 1142 NTESTT = NTESTT + NTEST 1143* 1144* Print out tests which fail. 1145* 1146 DO 220 JR = 1, NTEST 1147 IF( RESULT( JR ).GE.THRESH ) THEN 1148* 1149* If this is the first test to fail, 1150* print a header to the data file. 1151* 1152 IF( NERRS.EQ.0 ) THEN 1153 WRITE( NOUNIT, FMT = 9997 )'CGG' 1154* 1155* Matrix types 1156* 1157 WRITE( NOUNIT, FMT = 9996 ) 1158 WRITE( NOUNIT, FMT = 9995 ) 1159 WRITE( NOUNIT, FMT = 9994 )'Unitary' 1160* 1161* Tests performed 1162* 1163 WRITE( NOUNIT, FMT = 9993 )'unitary', '*', 1164 $ 'conjugate transpose', ( '*', J = 1, 10 ) 1165* 1166 END IF 1167 NERRS = NERRS + 1 1168 IF( RESULT( JR ).LT.10000.0 ) THEN 1169 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 1170 $ RESULT( JR ) 1171 ELSE 1172 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 1173 $ RESULT( JR ) 1174 END IF 1175 END IF 1176 220 CONTINUE 1177* 1178 230 CONTINUE 1179 240 CONTINUE 1180* 1181* Summary 1182* 1183 CALL SLASUM( 'CGG', NOUNIT, NERRS, NTESTT ) 1184 RETURN 1185* 1186 9999 FORMAT( ' CCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 1187 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 1188* 1189 9998 FORMAT( ' CCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ', 1190 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 1191 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 1192 $ ')' ) 1193* 1194 9997 FORMAT( 1X, A3, ' -- Complex Generalized eigenvalue problem' ) 1195* 1196 9996 FORMAT( ' Matrix types (see CCHKGG for details): ' ) 1197* 1198 9995 FORMAT( ' Special Matrices:', 23X, 1199 $ '(J''=transposed Jordan block)', 1200 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 1201 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 1202 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 1203 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 1204 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 1205 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 1206 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 1207 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 1208 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 1209 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 1210 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 1211 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 1212 $ '23=(small,large) 24=(small,small) 25=(large,large)', 1213 $ / ' 26=random O(1) matrices.' ) 1214* 1215 9993 FORMAT( / ' Tests performed: (H is Hessenberg, S is Schur, B, ', 1216 $ 'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A, 1217 $ ', l and r are the', / 20X, 1218 $ 'appropriate left and right eigenvectors, resp., a is', 1219 $ / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)', 1220 $ / ' 1 = | A - U H V', A, 1221 $ ' | / ( |A| n ulp ) 2 = | B - U T V', A, 1222 $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', A, 1223 $ ' | / ( n ulp ) 4 = | I - VV', A, 1224 $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', A, 1225 $ ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A, 1226 $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A, 1227 $ ' | / ( n ulp ) 8 = | I - ZZ', A, 1228 $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A, 1229 $ ' l | / const. 10 = max | ( b H - a T )', A, 1230 $ ' l | / const.', / 1231 $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H', 1232 $ ' - a T ) r | / const.', / 1X ) 1233* 1234 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 1235 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 1236 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 1237 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 ) 1238* 1239* End of CCHKGG 1240* 1241 END 1242