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