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