1*> \brief \b ZDRGEV3 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 ZDRGEV3( 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* DOUBLE PRECISION THRESH 20* .. 21* .. Array Arguments .. 22* LOGICAL DOTYPE( * ) 23* INTEGER ISEED( 4 ), NN( * ) 24* DOUBLE PRECISION RESULT( * ), RWORK( * ) 25* COMPLEX*16 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*> ZDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver 38*> routine ZGGEV3. 39*> 40*> ZGGEV3 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 ZDRGEV3 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 ZGGEV3: 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*> ZDRGEV3 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, ZDRGEV3 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 ZDRGES to continue the same random number 233*> sequence. 234*> \endverbatim 235*> 236*> \param[in] THRESH 237*> \verbatim 238*> THRESH is DOUBLE PRECISION 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*16 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*16 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*16 array, dimension (LDA, max(NN)) 280*> The Schur form matrix computed from A by ZGGEV3. 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*16 array, dimension (LDA, max(NN)) 288*> The upper triangular matrix computed from B by ZGGEV3. 289*> \endverbatim 290*> 291*> \param[out] Q 292*> \verbatim 293*> Q is COMPLEX*16 array, dimension (LDQ, max(NN)) 294*> The (left) eigenvectors matrix computed by ZGGEV3. 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*16 array, dimension( LDQ, max(NN) ) 307*> The (right) orthogonal matrix computed by ZGGEV3. 308*> \endverbatim 309*> 310*> \param[out] QE 311*> \verbatim 312*> QE is COMPLEX*16 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*16 array, dimension (max(NN)) 325*> \endverbatim 326*> 327*> \param[out] BETA 328*> \verbatim 329*> BETA is COMPLEX*16 array, dimension (max(NN)) 330*> 331*> The generalized eigenvalues of (A,B) computed by ZGGEV3. 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*16 array, dimension (max(NN)) 339*> \endverbatim 340*> 341*> \param[out] BETA1 342*> \verbatim 343*> BETA1 is COMPLEX*16 array, dimension (max(NN)) 344*> 345*> Like ALPHAR, ALPHAI, BETA, these arrays contain the 346*> eigenvalues of A and B, but those computed when ZGGEV3 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*16 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 DOUBLE PRECISION array, dimension (8*N) 365*> Real workspace. 366*> \endverbatim 367*> 368*> \param[out] RESULT 369*> \verbatim 370*> RESULT is DOUBLE PRECISION 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 complex16_eig 393* 394* ===================================================================== 395 SUBROUTINE ZDRGEV3( 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 DOUBLE PRECISION THRESH 408* .. 409* .. Array Arguments .. 410 LOGICAL DOTYPE( * ) 411 INTEGER ISEED( 4 ), NN( * ) 412 DOUBLE PRECISION RESULT( * ), RWORK( * ) 413 COMPLEX*16 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 DOUBLE PRECISION ZERO, ONE 423 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 424 COMPLEX*16 CZERO, CONE 425 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 426 $ CONE = ( 1.0D+0, 0.0D+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 DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV 436 COMPLEX*16 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 DOUBLE PRECISION RMAGN( 0: 3 ) 446* .. 447* .. External Functions .. 448 INTEGER ILAENV 449 DOUBLE PRECISION DLAMCH 450 COMPLEX*16 ZLARND 451 EXTERNAL ILAENV, DLAMCH, ZLARND 452* .. 453* .. External Subroutines .. 454 EXTERNAL ALASVM, DLABAD, XERBLA, ZGET52, ZGGEV3, ZLACPY, 455 $ ZLARFG, ZLASET, ZLATM4, ZUNM2R 456* .. 457* .. Intrinsic Functions .. 458 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, 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, 'ZGEQRF', ' ', NMAX, NMAX, -1, -1 ), 526 $ ILAENV( 1, 'ZUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ), 527 $ ILAENV( 1, 'ZUNGQR', ' ', 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( 'ZDRGEV3', -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 = DLAMCH( 'Precision' ) 546 SAFMIN = DLAMCH( 'Safe minimum' ) 547 SAFMIN = SAFMIN / ULP 548 SAFMAX = ONE / SAFMIN 549 CALL DLABAD( 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 / DBLE( 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* KZLASS: =1 means w/o rotation, =2 means w/ rotation, 591* =3 means random. 592* KATYPE: the "type" to be passed to ZLATM4 for computing A. 593* KAZERO: the pattern of zeros on the diagonal for A: 594* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 595* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 596* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 597* non-zero entries.) 598* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 599* =2: large, =3: small. 600* LASIGN: .TRUE. if the diagonal elements of A are to be 601* multiplied by a random magnitude 1 number. 602* KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. 603* KTRIAN: =0: don't fill in the upper triangle, =1: do. 604* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 605* RMAGN: used to implement KAMAGN and KBMAGN. 606* 607 IF( MTYPES.GT.MAXTYP ) 608 $ GO TO 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 ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 618 ELSE 619 IN = N 620 END IF 621 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 622 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 623 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 624 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 625 $ ISEED, A, LDA ) 626 IADD = KADD( KAZERO( JTYPE ) ) 627 IF( IADD.GT.0 .AND. IADD.LE.N ) 628 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 629* 630* Generate B (w/o rotation) 631* 632 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 633 IN = 2*( ( N-1 ) / 2 ) + 1 634 IF( IN.NE.N ) 635 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 636 ELSE 637 IN = N 638 END IF 639 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 640 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 641 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 642 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 643 $ ISEED, B, LDA ) 644 IADD = KADD( KBZERO( JTYPE ) ) 645 IF( IADD.NE.0 .AND. IADD.LE.N ) 646 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 647* 648 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 649* 650* Include rotations 651* 652* Generate Q, Z as Householder transformations times 653* a diagonal matrix. 654* 655 DO 40 JC = 1, N - 1 656 DO 30 JR = JC, N 657 Q( JR, JC ) = ZLARND( 3, ISEED ) 658 Z( JR, JC ) = ZLARND( 3, ISEED ) 659 30 CONTINUE 660 CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 661 $ WORK( JC ) ) 662 WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) ) 663 Q( JC, JC ) = CONE 664 CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 665 $ WORK( N+JC ) ) 666 WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) ) 667 Z( JC, JC ) = CONE 668 40 CONTINUE 669 CTEMP = ZLARND( 3, ISEED ) 670 Q( N, N ) = CONE 671 WORK( N ) = CZERO 672 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 673 CTEMP = ZLARND( 3, ISEED ) 674 Z( N, N ) = CONE 675 WORK( 2*N ) = CZERO 676 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 677* 678* Apply the diagonal matrices 679* 680 DO 60 JC = 1, N 681 DO 50 JR = 1, N 682 A( JR, JC ) = WORK( 2*N+JR )* 683 $ DCONJG( WORK( 3*N+JC ) )* 684 $ A( JR, JC ) 685 B( JR, JC ) = WORK( 2*N+JR )* 686 $ DCONJG( WORK( 3*N+JC ) )* 687 $ B( JR, JC ) 688 50 CONTINUE 689 60 CONTINUE 690 CALL ZUNM2R( '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 ZUNM2R( '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 ZUNM2R( '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 ZUNM2R( '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 $ ZLARND( 4, ISEED ) 715 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 716 $ ZLARND( 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 ZLAQZ0 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 ZGGEV3 to compute eigenvalues and eigenvectors. 745* 746 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 747 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 748 CALL ZGGEV3( '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 )'ZGGEV31', 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 ZGET52( .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', 'ZGGEV31', 764 $ RESULT( 2 ), N, JTYPE, IOLDSD 765 END IF 766* 767* Do the tests (3) and (4) 768* 769 CALL ZGET52( .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', 'ZGGEV31', 773 $ RESULT( 4 ), N, JTYPE, IOLDSD 774 END IF 775* 776* Do test (5) 777* 778 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 779 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 780 CALL ZGGEV3( '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 )'ZGGEV32', 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 test (6): Compute eigenvalues and left eigenvectors, 796* and test them 797* 798 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 799 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 800 CALL ZGGEV3( '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 )'ZGGEV33', IERR, N, JTYPE, 805 $ IOLDSD 806 INFO = ABS( IERR ) 807 GO TO 190 808 END IF 809* 810 DO 130 J = 1, N 811 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 812 $ BETA1( J ) )RESULT( 6 ) = ULPINV 813 130 CONTINUE 814* 815 DO 150 J = 1, N 816 DO 140 JC = 1, N 817 IF( Q( J, JC ).NE.QE( J, JC ) ) 818 $ RESULT( 6 ) = ULPINV 819 140 CONTINUE 820 150 CONTINUE 821* 822* Do test (7): Compute eigenvalues and right eigenvectors, 823* and test them 824* 825 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 826 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 827 CALL ZGGEV3( 'N', 'V', N, S, LDA, T, LDA, ALPHA1, BETA1, Q, 828 $ LDQ, QE, LDQE, WORK, LWORK, RWORK, IERR ) 829 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 830 RESULT( 1 ) = ULPINV 831 WRITE( NOUNIT, FMT = 9999 )'ZGGEV34', IERR, N, JTYPE, 832 $ IOLDSD 833 INFO = ABS( IERR ) 834 GO TO 190 835 END IF 836* 837 DO 160 J = 1, N 838 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 839 $ BETA1( J ) )RESULT( 7 ) = ULPINV 840 160 CONTINUE 841* 842 DO 180 J = 1, N 843 DO 170 JC = 1, N 844 IF( Z( J, JC ).NE.QE( J, JC ) ) 845 $ RESULT( 7 ) = ULPINV 846 170 CONTINUE 847 180 CONTINUE 848* 849* End of Loop -- Check for RESULT(j) > THRESH 850* 851 190 CONTINUE 852* 853 NTESTT = NTESTT + 7 854* 855* Print out tests which fail. 856* 857 DO 200 JR = 1, 7 858 IF( RESULT( JR ).GE.THRESH ) THEN 859* 860* If this is the first test to fail, 861* print a header to the data file. 862* 863 IF( NERRS.EQ.0 ) THEN 864 WRITE( NOUNIT, FMT = 9997 )'ZGV' 865* 866* Matrix types 867* 868 WRITE( NOUNIT, FMT = 9996 ) 869 WRITE( NOUNIT, FMT = 9995 ) 870 WRITE( NOUNIT, FMT = 9994 )'Orthogonal' 871* 872* Tests performed 873* 874 WRITE( NOUNIT, FMT = 9993 ) 875* 876 END IF 877 NERRS = NERRS + 1 878 IF( RESULT( JR ).LT.10000.0D0 ) THEN 879 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 880 $ RESULT( JR ) 881 ELSE 882 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 883 $ RESULT( JR ) 884 END IF 885 END IF 886 200 CONTINUE 887* 888 210 CONTINUE 889 220 CONTINUE 890* 891* Summary 892* 893 CALL ALASVM( 'ZGV3', NOUNIT, NERRS, NTESTT, 0 ) 894* 895 WORK( 1 ) = MAXWRK 896* 897 RETURN 898* 899 9999 FORMAT( ' ZDRGEV3: ', A, ' returned INFO=', I6, '.', / 3X, 'N=', 900 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 901* 902 9998 FORMAT( ' ZDRGEV3: ', A, ' Eigenvectors from ', A, 903 $ ' incorrectly normalized.', / ' Bits of error=', 0P, G10.3, 904 $ ',', 3X, 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 905 $ 3( I4, ',' ), I5, ')' ) 906* 907 9997 FORMAT( / 1X, A3, ' -- Complex Generalized eigenvalue problem ', 908 $ 'driver' ) 909* 910 9996 FORMAT( ' Matrix types (see ZDRGEV3 for details): ' ) 911* 912 9995 FORMAT( ' Special Matrices:', 23X, 913 $ '(J''=transposed Jordan block)', 914 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 915 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 916 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 917 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 918 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 919 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 920 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 921 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 922 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 923 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 924 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 925 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 926 $ '23=(small,large) 24=(small,small) 25=(large,large)', 927 $ / ' 26=random O(1) matrices.' ) 928* 929 9993 FORMAT( / ' Tests performed: ', 930 $ / ' 1 = max | ( b A - a B )''*l | / const.,', 931 $ / ' 2 = | |VR(i)| - 1 | / ulp,', 932 $ / ' 3 = max | ( b A - a B )*r | / const.', 933 $ / ' 4 = | |VL(i)| - 1 | / ulp,', 934 $ / ' 5 = 0 if W same no matter if r or l computed,', 935 $ / ' 6 = 0 if l same no matter if l computed,', 936 $ / ' 7 = 0 if r same no matter if r computed,', / 1X ) 937 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 938 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 939 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 940 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 ) 941* 942* End of ZDRGEV3 943* 944 END 945