1*> \brief \b DDRGEV3 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 DDRGEV3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 12* NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 13* ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1, 14* WORK, LWORK, 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 A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 25* $ ALPHI1( * ), ALPHR1( * ), B( LDA, * ), 26* $ BETA( * ), BETA1( * ), Q( LDQ, * ), 27* $ QE( LDQE, * ), RESULT( * ), S( LDA, * ), 28* $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> DDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver 38*> routine DGGEV3. 39*> 40*> DGGEV3 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 DDRGEV3 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 DGGEV3: 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*> DDRGEV3 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, DDRGEV3 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 DDRGEV3 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 DOUBLE PRECISION array, 257*> dimension(LDA, max(NN)) 258*> Used to hold the original A matrix. Used as input only 259*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 260*> DOTYPE(MAXTYP+1)=.TRUE. 261*> \endverbatim 262*> 263*> \param[in] LDA 264*> \verbatim 265*> LDA is INTEGER 266*> The leading dimension of A, B, S, and T. 267*> It must be at least 1 and at least max( NN ). 268*> \endverbatim 269*> 270*> \param[in,out] B 271*> \verbatim 272*> B is DOUBLE PRECISION array, 273*> dimension(LDA, max(NN)) 274*> Used to hold the original B matrix. Used as input only 275*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 276*> DOTYPE(MAXTYP+1)=.TRUE. 277*> \endverbatim 278*> 279*> \param[out] S 280*> \verbatim 281*> S is DOUBLE PRECISION array, 282*> dimension (LDA, max(NN)) 283*> The Schur form matrix computed from A by DGGEV3. On exit, S 284*> contains the Schur form matrix corresponding to the matrix 285*> in A. 286*> \endverbatim 287*> 288*> \param[out] T 289*> \verbatim 290*> T is DOUBLE PRECISION array, 291*> dimension (LDA, max(NN)) 292*> The upper triangular matrix computed from B by DGGEV3. 293*> \endverbatim 294*> 295*> \param[out] Q 296*> \verbatim 297*> Q is DOUBLE PRECISION array, 298*> dimension (LDQ, max(NN)) 299*> The (left) eigenvectors matrix computed by DGGEV3. 300*> \endverbatim 301*> 302*> \param[in] LDQ 303*> \verbatim 304*> LDQ is INTEGER 305*> The leading dimension of Q and Z. It must 306*> be at least 1 and at least max( NN ). 307*> \endverbatim 308*> 309*> \param[out] Z 310*> \verbatim 311*> Z is DOUBLE PRECISION array, dimension( LDQ, max(NN) ) 312*> The (right) orthogonal matrix computed by DGGEV3. 313*> \endverbatim 314*> 315*> \param[out] QE 316*> \verbatim 317*> QE is DOUBLE PRECISION array, dimension( LDQ, max(NN) ) 318*> QE holds the computed right or left eigenvectors. 319*> \endverbatim 320*> 321*> \param[in] LDQE 322*> \verbatim 323*> LDQE is INTEGER 324*> The leading dimension of QE. LDQE >= max(1,max(NN)). 325*> \endverbatim 326*> 327*> \param[out] ALPHAR 328*> \verbatim 329*> ALPHAR is DOUBLE PRECISION array, dimension (max(NN)) 330*> \endverbatim 331*> 332*> \param[out] ALPHAI 333*> \verbatim 334*> ALPHAI is DOUBLE PRECISION array, dimension (max(NN)) 335*> \endverbatim 336*> 337*> \param[out] BETA 338*> \verbatim 339*> BETA is DOUBLE PRECISION array, dimension (max(NN)) 340*> 341*> The generalized eigenvalues of (A,B) computed by DGGEV3. 342*> ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th 343*> generalized eigenvalue of A and B. 344*> \endverbatim 345*> 346*> \param[out] ALPHR1 347*> \verbatim 348*> ALPHR1 is DOUBLE PRECISION array, dimension (max(NN)) 349*> \endverbatim 350*> 351*> \param[out] ALPHI1 352*> \verbatim 353*> ALPHI1 is DOUBLE PRECISION array, dimension (max(NN)) 354*> \endverbatim 355*> 356*> \param[out] BETA1 357*> \verbatim 358*> BETA1 is DOUBLE PRECISION array, dimension (max(NN)) 359*> 360*> Like ALPHAR, ALPHAI, BETA, these arrays contain the 361*> eigenvalues of A and B, but those computed when DGGEV3 only 362*> computes a partial eigendecomposition, i.e. not the 363*> eigenvalues and left and right eigenvectors. 364*> \endverbatim 365*> 366*> \param[out] WORK 367*> \verbatim 368*> WORK is DOUBLE PRECISION array, dimension (LWORK) 369*> \endverbatim 370*> 371*> \param[in] LWORK 372*> \verbatim 373*> LWORK is INTEGER 374*> The number of entries in WORK. LWORK >= MAX( 8*N, N*(N+1) ). 375*> \endverbatim 376*> 377*> \param[out] RESULT 378*> \verbatim 379*> RESULT is DOUBLE PRECISION array, dimension (2) 380*> The values computed by the tests described above. 381*> The values are currently limited to 1/ulp, to avoid overflow. 382*> \endverbatim 383*> 384*> \param[out] INFO 385*> \verbatim 386*> INFO is INTEGER 387*> = 0: successful exit 388*> < 0: if INFO = -i, the i-th argument had an illegal value. 389*> > 0: A routine returned an error code. INFO is the 390*> absolute value of the INFO value returned. 391*> \endverbatim 392* 393* Authors: 394* ======== 395* 396*> \author Univ. of Tennessee 397*> \author Univ. of California Berkeley 398*> \author Univ. of Colorado Denver 399*> \author NAG Ltd. 400* 401*> \date February 2015 402* 403*> \ingroup double_eig 404* 405* ===================================================================== 406 SUBROUTINE DDRGEV3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 407 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 408 $ ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1, 409 $ WORK, LWORK, RESULT, INFO ) 410* 411* -- LAPACK test routine (version 3.6.1) -- 412* -- LAPACK is a software package provided by Univ. of Tennessee, -- 413* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 414* February 2015 415* 416* .. Scalar Arguments .. 417 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES, 418 $ NTYPES 419 DOUBLE PRECISION THRESH 420* .. 421* .. Array Arguments .. 422 LOGICAL DOTYPE( * ) 423 INTEGER ISEED( 4 ), NN( * ) 424 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 425 $ ALPHI1( * ), ALPHR1( * ), B( LDA, * ), 426 $ BETA( * ), BETA1( * ), Q( LDQ, * ), 427 $ QE( LDQE, * ), RESULT( * ), S( LDA, * ), 428 $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 429* .. 430* 431* ===================================================================== 432* 433* .. Parameters .. 434 DOUBLE PRECISION ZERO, ONE 435 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 436 INTEGER MAXTYP 437 PARAMETER ( MAXTYP = 26 ) 438* .. 439* .. Local Scalars .. 440 LOGICAL BADNN 441 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE, 442 $ MAXWRK, MINWRK, MTYPES, N, N1, NERRS, NMATS, 443 $ NMAX, NTESTT 444 DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV 445* .. 446* .. Local Arrays .. 447 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ), 448 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 449 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 450 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 451 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 452 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 453 DOUBLE PRECISION RMAGN( 0: 3 ) 454* .. 455* .. External Functions .. 456 INTEGER ILAENV 457 DOUBLE PRECISION DLAMCH, DLARND 458 EXTERNAL ILAENV, DLAMCH, DLARND 459* .. 460* .. External Subroutines .. 461 EXTERNAL ALASVM, DGET52, DGGEV3, DLABAD, DLACPY, DLARFG, 462 $ DLASET, DLATM4, DORM2R, XERBLA 463* .. 464* .. Intrinsic Functions .. 465 INTRINSIC ABS, DBLE, MAX, MIN, SIGN 466* .. 467* .. Data statements .. 468 DATA KCLASS / 15*1, 10*2, 1*3 / 469 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 470 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 471 DATA KADD / 0, 0, 0, 0, 3, 2 / 472 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 473 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 474 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 475 $ 1, 1, -4, 2, -4, 8*8, 0 / 476 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 477 $ 4*5, 4*3, 1 / 478 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 479 $ 4*6, 4*4, 1 / 480 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 481 $ 2, 1 / 482 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 483 $ 2, 1 / 484 DATA KTRIAN / 16*0, 10*1 / 485 DATA IASIGN / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0, 486 $ 5*2, 0 / 487 DATA IBSIGN / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 / 488* .. 489* .. Executable Statements .. 490* 491* Check for errors 492* 493 INFO = 0 494* 495 BADNN = .FALSE. 496 NMAX = 1 497 DO 10 J = 1, NSIZES 498 NMAX = MAX( NMAX, NN( J ) ) 499 IF( NN( J ).LT.0 ) 500 $ BADNN = .TRUE. 501 10 CONTINUE 502* 503 IF( NSIZES.LT.0 ) THEN 504 INFO = -1 505 ELSE IF( BADNN ) THEN 506 INFO = -2 507 ELSE IF( NTYPES.LT.0 ) THEN 508 INFO = -3 509 ELSE IF( THRESH.LT.ZERO ) THEN 510 INFO = -6 511 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 512 INFO = -9 513 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN 514 INFO = -14 515 ELSE IF( LDQE.LE.1 .OR. LDQE.LT.NMAX ) THEN 516 INFO = -17 517 END IF 518* 519* Compute workspace 520* (Note: Comments in the code beginning "Workspace:" describe the 521* minimal amount of workspace needed at that point in the code, 522* as well as the preferred amount for good performance. 523* NB refers to the optimal block size for the immediately 524* following subroutine, as returned by ILAENV. 525* 526 MINWRK = 1 527 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 528 MINWRK = MAX( 1, 8*NMAX, NMAX*( NMAX+1 ) ) 529 MAXWRK = 7*NMAX + NMAX*ILAENV( 1, 'DGEQRF', ' ', NMAX, 1, NMAX, 530 $ 0 ) 531 MAXWRK = MAX( MAXWRK, NMAX*( NMAX+1 ) ) 532 WORK( 1 ) = MAXWRK 533 END IF 534* 535 IF( LWORK.LT.MINWRK ) 536 $ INFO = -25 537* 538 IF( INFO.NE.0 ) THEN 539 CALL XERBLA( 'DDRGEV3', -INFO ) 540 RETURN 541 END IF 542* 543* Quick return if possible 544* 545 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 546 $ RETURN 547* 548 SAFMIN = DLAMCH( 'Safe minimum' ) 549 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 550 SAFMIN = SAFMIN / ULP 551 SAFMAX = ONE / SAFMIN 552 CALL DLABAD( SAFMIN, SAFMAX ) 553 ULPINV = ONE / ULP 554* 555* The values RMAGN(2:3) depend on N, see below. 556* 557 RMAGN( 0 ) = ZERO 558 RMAGN( 1 ) = ONE 559* 560* Loop over sizes, types 561* 562 NTESTT = 0 563 NERRS = 0 564 NMATS = 0 565* 566 DO 220 JSIZE = 1, NSIZES 567 N = NN( JSIZE ) 568 N1 = MAX( 1, N ) 569 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 ) 570 RMAGN( 3 ) = SAFMIN*ULPINV*N1 571* 572 IF( NSIZES.NE.1 ) THEN 573 MTYPES = MIN( MAXTYP, NTYPES ) 574 ELSE 575 MTYPES = MIN( MAXTYP+1, NTYPES ) 576 END IF 577* 578 DO 210 JTYPE = 1, MTYPES 579 IF( .NOT.DOTYPE( JTYPE ) ) 580 $ GO TO 210 581 NMATS = NMATS + 1 582* 583* Save ISEED in case of an error. 584* 585 DO 20 J = 1, 4 586 IOLDSD( J ) = ISEED( J ) 587 20 CONTINUE 588* 589* Generate test matrices A and B 590* 591* Description of control parameters: 592* 593* KZLASS: =1 means w/o rotation, =2 means w/ rotation, 594* =3 means random. 595* KATYPE: the "type" to be passed to DLATM4 for computing A. 596* KAZERO: the pattern of zeros on the diagonal for A: 597* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 598* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 599* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 600* non-zero entries.) 601* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 602* =2: large, =3: small. 603* IASIGN: 1 if the diagonal elements of A are to be 604* multiplied by a random magnitude 1 number, =2 if 605* randomly chosen diagonal blocks are to be rotated 606* to form 2x2 blocks. 607* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B. 608* KTRIAN: =0: don't fill in the upper triangle, =1: do. 609* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 610* RMAGN: used to implement KAMAGN and KBMAGN. 611* 612 IF( MTYPES.GT.MAXTYP ) 613 $ GO TO 100 614 IERR = 0 615 IF( KCLASS( JTYPE ).LT.3 ) THEN 616* 617* Generate A (w/o rotation) 618* 619 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 620 IN = 2*( ( N-1 ) / 2 ) + 1 621 IF( IN.NE.N ) 622 $ CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA ) 623 ELSE 624 IN = N 625 END IF 626 CALL DLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 627 $ KZ2( KAZERO( JTYPE ) ), IASIGN( JTYPE ), 628 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 629 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 630 $ ISEED, A, LDA ) 631 IADD = KADD( KAZERO( JTYPE ) ) 632 IF( IADD.GT.0 .AND. IADD.LE.N ) 633 $ A( IADD, IADD ) = ONE 634* 635* Generate B (w/o rotation) 636* 637 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 638 IN = 2*( ( N-1 ) / 2 ) + 1 639 IF( IN.NE.N ) 640 $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDA ) 641 ELSE 642 IN = N 643 END IF 644 CALL DLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 645 $ KZ2( KBZERO( JTYPE ) ), IBSIGN( JTYPE ), 646 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 647 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 648 $ ISEED, B, LDA ) 649 IADD = KADD( KBZERO( JTYPE ) ) 650 IF( IADD.NE.0 .AND. IADD.LE.N ) 651 $ B( IADD, IADD ) = ONE 652* 653 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 654* 655* Include rotations 656* 657* Generate Q, Z as Householder transformations times 658* a diagonal matrix. 659* 660 DO 40 JC = 1, N - 1 661 DO 30 JR = JC, N 662 Q( JR, JC ) = DLARND( 3, ISEED ) 663 Z( JR, JC ) = DLARND( 3, ISEED ) 664 30 CONTINUE 665 CALL DLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 666 $ WORK( JC ) ) 667 WORK( 2*N+JC ) = SIGN( ONE, Q( JC, JC ) ) 668 Q( JC, JC ) = ONE 669 CALL DLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 670 $ WORK( N+JC ) ) 671 WORK( 3*N+JC ) = SIGN( ONE, Z( JC, JC ) ) 672 Z( JC, JC ) = ONE 673 40 CONTINUE 674 Q( N, N ) = ONE 675 WORK( N ) = ZERO 676 WORK( 3*N ) = SIGN( ONE, DLARND( 2, ISEED ) ) 677 Z( N, N ) = ONE 678 WORK( 2*N ) = ZERO 679 WORK( 4*N ) = SIGN( ONE, DLARND( 2, ISEED ) ) 680* 681* Apply the diagonal matrices 682* 683 DO 60 JC = 1, N 684 DO 50 JR = 1, N 685 A( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )* 686 $ A( JR, JC ) 687 B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )* 688 $ B( JR, JC ) 689 50 CONTINUE 690 60 CONTINUE 691 CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A, 692 $ LDA, WORK( 2*N+1 ), IERR ) 693 IF( IERR.NE.0 ) 694 $ GO TO 90 695 CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ), 696 $ A, LDA, WORK( 2*N+1 ), IERR ) 697 IF( IERR.NE.0 ) 698 $ GO TO 90 699 CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B, 700 $ LDA, WORK( 2*N+1 ), IERR ) 701 IF( IERR.NE.0 ) 702 $ GO TO 90 703 CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ), 704 $ B, LDA, WORK( 2*N+1 ), IERR ) 705 IF( IERR.NE.0 ) 706 $ GO TO 90 707 END IF 708 ELSE 709* 710* Random matrices 711* 712 DO 80 JC = 1, N 713 DO 70 JR = 1, N 714 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 715 $ DLARND( 2, ISEED ) 716 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 717 $ DLARND( 2, ISEED ) 718 70 CONTINUE 719 80 CONTINUE 720 END IF 721* 722 90 CONTINUE 723* 724 IF( IERR.NE.0 ) THEN 725 WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE, 726 $ IOLDSD 727 INFO = ABS( IERR ) 728 RETURN 729 END IF 730* 731 100 CONTINUE 732* 733 DO 110 I = 1, 7 734 RESULT( I ) = -ONE 735 110 CONTINUE 736* 737* Call DGGEV3 to compute eigenvalues and eigenvectors. 738* 739 CALL DLACPY( ' ', N, N, A, LDA, S, LDA ) 740 CALL DLACPY( ' ', N, N, B, LDA, T, LDA ) 741 CALL DGGEV3( 'V', 'V', N, S, LDA, T, LDA, ALPHAR, ALPHAI, 742 $ BETA, Q, LDQ, Z, LDQ, WORK, LWORK, IERR ) 743 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 744 RESULT( 1 ) = ULPINV 745 WRITE( NOUNIT, FMT = 9999 )'DGGEV31', IERR, N, JTYPE, 746 $ IOLDSD 747 INFO = ABS( IERR ) 748 GO TO 190 749 END IF 750* 751* Do the tests (1) and (2) 752* 753 CALL DGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHAR, 754 $ ALPHAI, BETA, WORK, RESULT( 1 ) ) 755 IF( RESULT( 2 ).GT.THRESH ) THEN 756 WRITE( NOUNIT, FMT = 9998 )'Left', 'DGGEV31', 757 $ RESULT( 2 ), N, JTYPE, IOLDSD 758 END IF 759* 760* Do the tests (3) and (4) 761* 762 CALL DGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHAR, 763 $ ALPHAI, BETA, WORK, RESULT( 3 ) ) 764 IF( RESULT( 4 ).GT.THRESH ) THEN 765 WRITE( NOUNIT, FMT = 9998 )'Right', 'DGGEV31', 766 $ RESULT( 4 ), N, JTYPE, IOLDSD 767 END IF 768* 769* Do the test (5) 770* 771 CALL DLACPY( ' ', N, N, A, LDA, S, LDA ) 772 CALL DLACPY( ' ', N, N, B, LDA, T, LDA ) 773 CALL DGGEV3( 'N', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 774 $ BETA1, Q, LDQ, Z, LDQ, WORK, LWORK, IERR ) 775 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 776 RESULT( 1 ) = ULPINV 777 WRITE( NOUNIT, FMT = 9999 )'DGGEV32', IERR, N, JTYPE, 778 $ IOLDSD 779 INFO = ABS( IERR ) 780 GO TO 190 781 END IF 782* 783 DO 120 J = 1, N 784 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 785 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) )RESULT( 5 ) 786 $ = ULPINV 787 120 CONTINUE 788* 789* Do the test (6): Compute eigenvalues and left eigenvectors, 790* and test them 791* 792 CALL DLACPY( ' ', N, N, A, LDA, S, LDA ) 793 CALL DLACPY( ' ', N, N, B, LDA, T, LDA ) 794 CALL DGGEV3( 'V', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 795 $ BETA1, QE, LDQE, Z, LDQ, WORK, LWORK, IERR ) 796 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 797 RESULT( 1 ) = ULPINV 798 WRITE( NOUNIT, FMT = 9999 )'DGGEV33', IERR, N, JTYPE, 799 $ IOLDSD 800 INFO = ABS( IERR ) 801 GO TO 190 802 END IF 803* 804 DO 130 J = 1, N 805 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 806 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) )RESULT( 6 ) 807 $ = ULPINV 808 130 CONTINUE 809* 810 DO 150 J = 1, N 811 DO 140 JC = 1, N 812 IF( Q( J, JC ).NE.QE( J, JC ) ) 813 $ RESULT( 6 ) = ULPINV 814 140 CONTINUE 815 150 CONTINUE 816* 817* DO the test (7): Compute eigenvalues and right eigenvectors, 818* and test them 819* 820 CALL DLACPY( ' ', N, N, A, LDA, S, LDA ) 821 CALL DLACPY( ' ', N, N, B, LDA, T, LDA ) 822 CALL DGGEV3( 'N', 'V', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 823 $ BETA1, Q, LDQ, QE, LDQE, WORK, LWORK, IERR ) 824 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 825 RESULT( 1 ) = ULPINV 826 WRITE( NOUNIT, FMT = 9999 )'DGGEV34', IERR, N, JTYPE, 827 $ IOLDSD 828 INFO = ABS( IERR ) 829 GO TO 190 830 END IF 831* 832 DO 160 J = 1, N 833 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 834 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) )RESULT( 7 ) 835 $ = ULPINV 836 160 CONTINUE 837* 838 DO 180 J = 1, N 839 DO 170 JC = 1, N 840 IF( Z( J, JC ).NE.QE( J, JC ) ) 841 $ RESULT( 7 ) = ULPINV 842 170 CONTINUE 843 180 CONTINUE 844* 845* End of Loop -- Check for RESULT(j) > THRESH 846* 847 190 CONTINUE 848* 849 NTESTT = NTESTT + 7 850* 851* Print out tests which fail. 852* 853 DO 200 JR = 1, 7 854 IF( RESULT( JR ).GE.THRESH ) THEN 855* 856* If this is the first test to fail, 857* print a header to the data file. 858* 859 IF( NERRS.EQ.0 ) THEN 860 WRITE( NOUNIT, FMT = 9997 )'DGV' 861* 862* Matrix types 863* 864 WRITE( NOUNIT, FMT = 9996 ) 865 WRITE( NOUNIT, FMT = 9995 ) 866 WRITE( NOUNIT, FMT = 9994 )'Orthogonal' 867* 868* Tests performed 869* 870 WRITE( NOUNIT, FMT = 9993 ) 871* 872 END IF 873 NERRS = NERRS + 1 874 IF( RESULT( JR ).LT.10000.0D0 ) THEN 875 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 876 $ RESULT( JR ) 877 ELSE 878 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 879 $ RESULT( JR ) 880 END IF 881 END IF 882 200 CONTINUE 883* 884 210 CONTINUE 885 220 CONTINUE 886* 887* Summary 888* 889 CALL ALASVM( 'DGV', NOUNIT, NERRS, NTESTT, 0 ) 890* 891 WORK( 1 ) = MAXWRK 892* 893 RETURN 894* 895 9999 FORMAT( ' DDRGEV3: ', A, ' returned INFO=', I6, '.', / 3X, 'N=', 896 $ I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' ) 897* 898 9998 FORMAT( ' DDRGEV3: ', A, ' Eigenvectors from ', A, 899 $ ' incorrectly normalized.', / ' Bits of error=', 0P, G10.3, 900 $ ',', 3X, 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 901 $ 4( I4, ',' ), I5, ')' ) 902* 903 9997 FORMAT( / 1X, A3, ' -- Real Generalized eigenvalue problem driver' 904 $ ) 905* 906 9996 FORMAT( ' Matrix types (see DDRGEV3 for details): ' ) 907* 908 9995 FORMAT( ' Special Matrices:', 23X, 909 $ '(J''=transposed Jordan block)', 910 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 911 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 912 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 913 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 914 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 915 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 916 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 917 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 918 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 919 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 920 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 921 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 922 $ '23=(small,large) 24=(small,small) 25=(large,large)', 923 $ / ' 26=random O(1) matrices.' ) 924* 925 9993 FORMAT( / ' Tests performed: ', 926 $ / ' 1 = max | ( b A - a B )''*l | / const.,', 927 $ / ' 2 = | |VR(i)| - 1 | / ulp,', 928 $ / ' 3 = max | ( b A - a B )*r | / const.', 929 $ / ' 4 = | |VL(i)| - 1 | / ulp,', 930 $ / ' 5 = 0 if W same no matter if r or l computed,', 931 $ / ' 6 = 0 if l same no matter if l computed,', 932 $ / ' 7 = 0 if r same no matter if r computed,', / 1X ) 933 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 934 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 935 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 936 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 ) 937* 938* End of DDRGEV3 939* 940 END 941