1*> \brief \b SDRGEV 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 SDRGEV( 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* REAL THRESH 20* .. 21* .. Array Arguments .. 22* LOGICAL DOTYPE( * ) 23* INTEGER ISEED( 4 ), NN( * ) 24* REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ), 25* $ ALPHAR( * ), 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*> SDRGEV checks the nonsymmetric generalized eigenvalue problem driver 38*> routine SGGEV. 39*> 40*> SGGEV 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 SDRGEV 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 SGGEV: 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*> SDRGES 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, SDRGES 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 SDRGES 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 REAL 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 REAL 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 REAL array, 282*> dimension (LDA, max(NN)) 283*> The Schur form matrix computed from A by SGGES. 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 REAL array, 291*> dimension (LDA, max(NN)) 292*> The upper triangular matrix computed from B by SGGES. 293*> \endverbatim 294*> 295*> \param[out] Q 296*> \verbatim 297*> Q is REAL array, 298*> dimension (LDQ, max(NN)) 299*> The (left) eigenvectors matrix computed by SGGEV. 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 REAL array, dimension( LDQ, max(NN) ) 312*> The (right) orthogonal matrix computed by SGGES. 313*> \endverbatim 314*> 315*> \param[out] QE 316*> \verbatim 317*> QE is REAL 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 REAL array, dimension (max(NN)) 330*> \endverbatim 331*> 332*> \param[out] ALPHAI 333*> \verbatim 334*> ALPHAI is REAL array, dimension (max(NN)) 335*> \endverbatim 336*> 337*> \param[out] BETA 338*> \verbatim 339*> BETA is REAL array, dimension (max(NN)) 340*> \verbatim 341*> The generalized eigenvalues of (A,B) computed by SGGEV. 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 REAL array, dimension (max(NN)) 349*> \endverbatim 350*> 351*> \param[out] ALPHI1 352*> \verbatim 353*> ALPHI1 is REAL array, dimension (max(NN)) 354*> \endverbatim 355*> 356*> \param[out] BETA1 357*> \verbatim 358*> BETA1 is REAL array, dimension (max(NN)) 359*> 360*> Like ALPHAR, ALPHAI, BETA, these arrays contain the 361*> eigenvalues of A and B, but those computed when SGGEV 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 REAL 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 REAL 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*> \ingroup single_eig 402* 403* ===================================================================== 404 SUBROUTINE SDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 405 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 406 $ ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1, 407 $ WORK, LWORK, RESULT, INFO ) 408* 409* -- LAPACK test routine -- 410* -- LAPACK is a software package provided by Univ. of Tennessee, -- 411* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 412* 413* .. Scalar Arguments .. 414 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES, 415 $ NTYPES 416 REAL THRESH 417* .. 418* .. Array Arguments .. 419 LOGICAL DOTYPE( * ) 420 INTEGER ISEED( 4 ), NN( * ) 421 REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ), 422 $ ALPHAR( * ), ALPHR1( * ), B( LDA, * ), 423 $ BETA( * ), BETA1( * ), Q( LDQ, * ), 424 $ QE( LDQE, * ), RESULT( * ), S( LDA, * ), 425 $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 426* .. 427* 428* ===================================================================== 429* 430* .. Parameters .. 431 REAL ZERO, ONE 432 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 433 INTEGER MAXTYP 434 PARAMETER ( MAXTYP = 26 ) 435* .. 436* .. Local Scalars .. 437 LOGICAL BADNN 438 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE, 439 $ MAXWRK, MINWRK, MTYPES, N, N1, NERRS, NMATS, 440 $ NMAX, NTESTT 441 REAL SAFMAX, SAFMIN, ULP, ULPINV 442* .. 443* .. Local Arrays .. 444 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ), 445 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 446 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 447 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 448 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 449 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 450 REAL RMAGN( 0: 3 ) 451* .. 452* .. External Functions .. 453 INTEGER ILAENV 454 REAL SLAMCH, SLARND 455 EXTERNAL ILAENV, SLAMCH, SLARND 456* .. 457* .. External Subroutines .. 458 EXTERNAL ALASVM, SGET52, SGGEV, SLABAD, SLACPY, SLARFG, 459 $ SLASET, SLATM4, SORM2R, XERBLA 460* .. 461* .. Intrinsic Functions .. 462 INTRINSIC ABS, MAX, MIN, REAL, SIGN 463* .. 464* .. Data statements .. 465 DATA KCLASS / 15*1, 10*2, 1*3 / 466 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 467 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 468 DATA KADD / 0, 0, 0, 0, 3, 2 / 469 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 470 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 471 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 472 $ 1, 1, -4, 2, -4, 8*8, 0 / 473 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 474 $ 4*5, 4*3, 1 / 475 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 476 $ 4*6, 4*4, 1 / 477 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 478 $ 2, 1 / 479 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 480 $ 2, 1 / 481 DATA KTRIAN / 16*0, 10*1 / 482 DATA IASIGN / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0, 483 $ 5*2, 0 / 484 DATA IBSIGN / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 / 485* .. 486* .. Executable Statements .. 487* 488* Check for errors 489* 490 INFO = 0 491* 492 BADNN = .FALSE. 493 NMAX = 1 494 DO 10 J = 1, NSIZES 495 NMAX = MAX( NMAX, NN( J ) ) 496 IF( NN( J ).LT.0 ) 497 $ BADNN = .TRUE. 498 10 CONTINUE 499* 500 IF( NSIZES.LT.0 ) THEN 501 INFO = -1 502 ELSE IF( BADNN ) THEN 503 INFO = -2 504 ELSE IF( NTYPES.LT.0 ) THEN 505 INFO = -3 506 ELSE IF( THRESH.LT.ZERO ) THEN 507 INFO = -6 508 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 509 INFO = -9 510 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN 511 INFO = -14 512 ELSE IF( LDQE.LE.1 .OR. LDQE.LT.NMAX ) THEN 513 INFO = -17 514 END IF 515* 516* Compute workspace 517* (Note: Comments in the code beginning "Workspace:" describe the 518* minimal amount of workspace needed at that point in the code, 519* as well as the preferred amount for good performance. 520* NB refers to the optimal block size for the immediately 521* following subroutine, as returned by ILAENV. 522* 523 MINWRK = 1 524 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 525 MINWRK = MAX( 1, 8*NMAX, NMAX*( NMAX+1 ) ) 526 MAXWRK = 7*NMAX + NMAX*ILAENV( 1, 'SGEQRF', ' ', NMAX, 1, NMAX, 527 $ 0 ) 528 MAXWRK = MAX( MAXWRK, NMAX*( NMAX+1 ) ) 529 WORK( 1 ) = MAXWRK 530 END IF 531* 532 IF( LWORK.LT.MINWRK ) 533 $ INFO = -25 534* 535 IF( INFO.NE.0 ) THEN 536 CALL XERBLA( 'SDRGEV', -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 SAFMIN = SLAMCH( 'Safe minimum' ) 546 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 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 SLATM4 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* IASIGN: 1 if the diagonal elements of A are to be 601* multiplied by a random magnitude 1 number, =2 if 602* randomly chosen diagonal blocks are to be rotated 603* to form 2x2 blocks. 604* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B. 605* KTRIAN: =0: don't fill in the upper triangle, =1: do. 606* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 607* RMAGN: used to implement KAMAGN and KBMAGN. 608* 609 IF( MTYPES.GT.MAXTYP ) 610 $ GO TO 100 611 IERR = 0 612 IF( KCLASS( JTYPE ).LT.3 ) THEN 613* 614* Generate A (w/o rotation) 615* 616 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 617 IN = 2*( ( N-1 ) / 2 ) + 1 618 IF( IN.NE.N ) 619 $ CALL SLASET( 'Full', N, N, ZERO, ZERO, A, LDA ) 620 ELSE 621 IN = N 622 END IF 623 CALL SLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 624 $ KZ2( KAZERO( JTYPE ) ), IASIGN( JTYPE ), 625 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 626 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 627 $ ISEED, A, LDA ) 628 IADD = KADD( KAZERO( JTYPE ) ) 629 IF( IADD.GT.0 .AND. IADD.LE.N ) 630 $ A( IADD, IADD ) = ONE 631* 632* Generate B (w/o rotation) 633* 634 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 635 IN = 2*( ( N-1 ) / 2 ) + 1 636 IF( IN.NE.N ) 637 $ CALL SLASET( 'Full', N, N, ZERO, ZERO, B, LDA ) 638 ELSE 639 IN = N 640 END IF 641 CALL SLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 642 $ KZ2( KBZERO( JTYPE ) ), IBSIGN( JTYPE ), 643 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 644 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 645 $ ISEED, B, LDA ) 646 IADD = KADD( KBZERO( JTYPE ) ) 647 IF( IADD.NE.0 .AND. IADD.LE.N ) 648 $ B( IADD, IADD ) = ONE 649* 650 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 651* 652* Include rotations 653* 654* Generate Q, Z as Householder transformations times 655* a diagonal matrix. 656* 657 DO 40 JC = 1, N - 1 658 DO 30 JR = JC, N 659 Q( JR, JC ) = SLARND( 3, ISEED ) 660 Z( JR, JC ) = SLARND( 3, ISEED ) 661 30 CONTINUE 662 CALL SLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 663 $ WORK( JC ) ) 664 WORK( 2*N+JC ) = SIGN( ONE, Q( JC, JC ) ) 665 Q( JC, JC ) = ONE 666 CALL SLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 667 $ WORK( N+JC ) ) 668 WORK( 3*N+JC ) = SIGN( ONE, Z( JC, JC ) ) 669 Z( JC, JC ) = ONE 670 40 CONTINUE 671 Q( N, N ) = ONE 672 WORK( N ) = ZERO 673 WORK( 3*N ) = SIGN( ONE, SLARND( 2, ISEED ) ) 674 Z( N, N ) = ONE 675 WORK( 2*N ) = ZERO 676 WORK( 4*N ) = SIGN( ONE, SLARND( 2, ISEED ) ) 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 )*WORK( 3*N+JC )* 683 $ A( JR, JC ) 684 B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )* 685 $ B( JR, JC ) 686 50 CONTINUE 687 60 CONTINUE 688 CALL SORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A, 689 $ LDA, WORK( 2*N+1 ), IERR ) 690 IF( IERR.NE.0 ) 691 $ GO TO 90 692 CALL SORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ), 693 $ A, LDA, WORK( 2*N+1 ), IERR ) 694 IF( IERR.NE.0 ) 695 $ GO TO 90 696 CALL SORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B, 697 $ LDA, WORK( 2*N+1 ), IERR ) 698 IF( IERR.NE.0 ) 699 $ GO TO 90 700 CALL SORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ), 701 $ B, LDA, WORK( 2*N+1 ), IERR ) 702 IF( IERR.NE.0 ) 703 $ GO TO 90 704 END IF 705 ELSE 706* 707* Random matrices 708* 709 DO 80 JC = 1, N 710 DO 70 JR = 1, N 711 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 712 $ SLARND( 2, ISEED ) 713 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 714 $ SLARND( 2, ISEED ) 715 70 CONTINUE 716 80 CONTINUE 717 END IF 718* 719 90 CONTINUE 720* 721 IF( IERR.NE.0 ) THEN 722 WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE, 723 $ IOLDSD 724 INFO = ABS( IERR ) 725 RETURN 726 END IF 727* 728 100 CONTINUE 729* 730 DO 110 I = 1, 7 731 RESULT( I ) = -ONE 732 110 CONTINUE 733* 734* Call SGGEV to compute eigenvalues and eigenvectors. 735* 736 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 737 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 738 CALL SGGEV( 'V', 'V', N, S, LDA, T, LDA, ALPHAR, ALPHAI, 739 $ BETA, Q, LDQ, Z, LDQ, WORK, LWORK, IERR ) 740 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 741 RESULT( 1 ) = ULPINV 742 WRITE( NOUNIT, FMT = 9999 )'SGGEV1', IERR, N, JTYPE, 743 $ IOLDSD 744 INFO = ABS( IERR ) 745 GO TO 190 746 END IF 747* 748* Do the tests (1) and (2) 749* 750 CALL SGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHAR, 751 $ ALPHAI, BETA, WORK, RESULT( 1 ) ) 752 IF( RESULT( 2 ).GT.THRESH ) THEN 753 WRITE( NOUNIT, FMT = 9998 )'Left', 'SGGEV1', 754 $ RESULT( 2 ), N, JTYPE, IOLDSD 755 END IF 756* 757* Do the tests (3) and (4) 758* 759 CALL SGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHAR, 760 $ ALPHAI, BETA, WORK, RESULT( 3 ) ) 761 IF( RESULT( 4 ).GT.THRESH ) THEN 762 WRITE( NOUNIT, FMT = 9998 )'Right', 'SGGEV1', 763 $ RESULT( 4 ), N, JTYPE, IOLDSD 764 END IF 765* 766* Do the test (5) 767* 768 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 769 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 770 CALL SGGEV( 'N', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 771 $ BETA1, Q, LDQ, Z, LDQ, WORK, LWORK, IERR ) 772 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 773 RESULT( 1 ) = ULPINV 774 WRITE( NOUNIT, FMT = 9999 )'SGGEV2', IERR, N, JTYPE, 775 $ IOLDSD 776 INFO = ABS( IERR ) 777 GO TO 190 778 END IF 779* 780 DO 120 J = 1, N 781 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 782 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) ) 783 $ RESULT( 5 ) = ULPINV 784 120 CONTINUE 785* 786* Do the test (6): Compute eigenvalues and left eigenvectors, 787* and test them 788* 789 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 790 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 791 CALL SGGEV( 'V', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 792 $ BETA1, QE, LDQE, Z, LDQ, WORK, LWORK, IERR ) 793 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 794 RESULT( 1 ) = ULPINV 795 WRITE( NOUNIT, FMT = 9999 )'SGGEV3', IERR, N, JTYPE, 796 $ IOLDSD 797 INFO = ABS( IERR ) 798 GO TO 190 799 END IF 800* 801 DO 130 J = 1, N 802 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 803 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) ) 804 $ RESULT( 6 ) = ULPINV 805 130 CONTINUE 806* 807 DO 150 J = 1, N 808 DO 140 JC = 1, N 809 IF( Q( J, JC ).NE.QE( J, JC ) ) 810 $ RESULT( 6 ) = ULPINV 811 140 CONTINUE 812 150 CONTINUE 813* 814* DO the test (7): Compute eigenvalues and right eigenvectors, 815* and test them 816* 817 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 818 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 819 CALL SGGEV( 'N', 'V', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 820 $ BETA1, Q, LDQ, QE, LDQE, WORK, LWORK, IERR ) 821 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 822 RESULT( 1 ) = ULPINV 823 WRITE( NOUNIT, FMT = 9999 )'SGGEV4', IERR, N, JTYPE, 824 $ IOLDSD 825 INFO = ABS( IERR ) 826 GO TO 190 827 END IF 828* 829 DO 160 J = 1, N 830 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 831 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) ) 832 $ RESULT( 7 ) = ULPINV 833 160 CONTINUE 834* 835 DO 180 J = 1, N 836 DO 170 JC = 1, N 837 IF( Z( J, JC ).NE.QE( J, JC ) ) 838 $ RESULT( 7 ) = ULPINV 839 170 CONTINUE 840 180 CONTINUE 841* 842* End of Loop -- Check for RESULT(j) > THRESH 843* 844 190 CONTINUE 845* 846 NTESTT = NTESTT + 7 847* 848* Print out tests which fail. 849* 850 DO 200 JR = 1, 7 851 IF( RESULT( JR ).GE.THRESH ) THEN 852* 853* If this is the first test to fail, 854* print a header to the data file. 855* 856 IF( NERRS.EQ.0 ) THEN 857 WRITE( NOUNIT, FMT = 9997 )'SGV' 858* 859* Matrix types 860* 861 WRITE( NOUNIT, FMT = 9996 ) 862 WRITE( NOUNIT, FMT = 9995 ) 863 WRITE( NOUNIT, FMT = 9994 )'Orthogonal' 864* 865* Tests performed 866* 867 WRITE( NOUNIT, FMT = 9993 ) 868* 869 END IF 870 NERRS = NERRS + 1 871 IF( RESULT( JR ).LT.10000.0 ) THEN 872 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 873 $ RESULT( JR ) 874 ELSE 875 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 876 $ RESULT( JR ) 877 END IF 878 END IF 879 200 CONTINUE 880* 881 210 CONTINUE 882 220 CONTINUE 883* 884* Summary 885* 886 CALL ALASVM( 'SGV', NOUNIT, NERRS, NTESTT, 0 ) 887* 888 WORK( 1 ) = MAXWRK 889* 890 RETURN 891* 892 9999 FORMAT( ' SDRGEV: ', A, ' returned INFO=', I6, '.', / 3X, 'N=', 893 $ I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' ) 894* 895 9998 FORMAT( ' SDRGEV: ', A, ' Eigenvectors from ', A, ' incorrectly ', 896 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 3X, 897 $ 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 4( I4, ',' ), I5, 898 $ ')' ) 899* 900 9997 FORMAT( / 1X, A3, ' -- Real Generalized eigenvalue problem driver' 901 $ ) 902* 903 9996 FORMAT( ' Matrix types (see SDRGEV for details): ' ) 904* 905 9995 FORMAT( ' Special Matrices:', 23X, 906 $ '(J''=transposed Jordan block)', 907 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 908 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 909 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 910 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 911 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 912 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 913 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 914 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 915 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 916 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 917 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 918 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 919 $ '23=(small,large) 24=(small,small) 25=(large,large)', 920 $ / ' 26=random O(1) matrices.' ) 921* 922 9993 FORMAT( / ' Tests performed: ', 923 $ / ' 1 = max | ( b A - a B )''*l | / const.,', 924 $ / ' 2 = | |VR(i)| - 1 | / ulp,', 925 $ / ' 3 = max | ( b A - a B )*r | / const.', 926 $ / ' 4 = | |VL(i)| - 1 | / ulp,', 927 $ / ' 5 = 0 if W same no matter if r or l computed,', 928 $ / ' 6 = 0 if l same no matter if l computed,', 929 $ / ' 7 = 0 if r same no matter if r computed,', / 1X ) 930 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 931 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 932 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 933 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 ) 934* 935* End of SDRGEV 936* 937 END 938