1*> \brief \b DDRGSX 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 DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, 12* BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S, 13* WORK, LWORK, IWORK, LIWORK, BWORK, INFO ) 14* 15* .. Scalar Arguments .. 16* INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN, 17* $ NOUT, NSIZE 18* DOUBLE PRECISION THRESH 19* .. 20* .. Array Arguments .. 21* LOGICAL BWORK( * ) 22* INTEGER IWORK( * ) 23* DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ), 24* $ ALPHAR( * ), B( LDA, * ), BETA( * ), 25* $ BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ), 26* $ WORK( * ), Z( LDA, * ) 27* .. 28* 29* 30*> \par Purpose: 31* ============= 32*> 33*> \verbatim 34*> 35*> DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form) 36*> problem expert driver DGGESX. 37*> 38*> DGGESX 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(1),beta(1)), ..., 44*> (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the 45*> characteristic equation 46*> 47*> det( A - w(j) B ) = 0 48*> 49*> Optionally it also reorders the eigenvalues so that a selected 50*> cluster of eigenvalues appears in the leading diagonal block of the 51*> Schur forms; computes a reciprocal condition number for the average 52*> of the selected eigenvalues; and computes a reciprocal condition 53*> number for the right and left deflating subspaces corresponding to 54*> the selected eigenvalues. 55*> 56*> When DDRGSX is called with NSIZE > 0, five (5) types of built-in 57*> matrix pairs are used to test the routine DGGESX. 58*> 59*> When DDRGSX is called with NSIZE = 0, it reads in test matrix data 60*> to test DGGESX. 61*> 62*> For each matrix pair, the following tests will be performed and 63*> compared with the threshold THRESH except for the tests (7) and (9): 64*> 65*> (1) | A - Q S Z' | / ( |A| n ulp ) 66*> 67*> (2) | B - Q T Z' | / ( |B| n ulp ) 68*> 69*> (3) | I - QQ' | / ( n ulp ) 70*> 71*> (4) | I - ZZ' | / ( n ulp ) 72*> 73*> (5) if A is in Schur form (i.e. quasi-triangular form) 74*> 75*> (6) maximum over j of D(j) where: 76*> 77*> if alpha(j) is real: 78*> |alpha(j) - S(j,j)| |beta(j) - T(j,j)| 79*> D(j) = ------------------------ + ----------------------- 80*> max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|) 81*> 82*> if alpha(j) is complex: 83*> | det( s S - w T ) | 84*> D(j) = --------------------------------------------------- 85*> ulp max( s norm(S), |w| norm(T) )*norm( s S - w T ) 86*> 87*> and S and T are here the 2 x 2 diagonal blocks of S and T 88*> corresponding to the j-th and j+1-th eigenvalues. 89*> 90*> (7) if sorting worked and SDIM is the number of eigenvalues 91*> which were selected. 92*> 93*> (8) the estimated value DIF does not differ from the true values of 94*> Difu and Difl more than a factor 10*THRESH. If the estimate DIF 95*> equals zero the corresponding true values of Difu and Difl 96*> should be less than EPS*norm(A, B). If the true value of Difu 97*> and Difl equal zero, the estimate DIF should be less than 98*> EPS*norm(A, B). 99*> 100*> (9) If INFO = N+3 is returned by DGGESX, the reordering "failed" 101*> and we check that DIF = PL = PR = 0 and that the true value of 102*> Difu and Difl is < EPS*norm(A, B). We count the events when 103*> INFO=N+3. 104*> 105*> For read-in test matrices, the above tests are run except that the 106*> exact value for DIF (and PL) is input data. Additionally, there is 107*> one more test run for read-in test matrices: 108*> 109*> (10) the estimated value PL does not differ from the true value of 110*> PLTRU more than a factor THRESH. If the estimate PL equals 111*> zero the corresponding true value of PLTRU should be less than 112*> EPS*norm(A, B). If the true value of PLTRU equal zero, the 113*> estimate PL should be less than EPS*norm(A, B). 114*> 115*> Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1) 116*> matrix pairs are generated and tested. NSIZE should be kept small. 117*> 118*> SVD (routine DGESVD) is used for computing the true value of DIF_u 119*> and DIF_l when testing the built-in test problems. 120*> 121*> Built-in Test Matrices 122*> ====================== 123*> 124*> All built-in test matrices are the 2 by 2 block of triangular 125*> matrices 126*> 127*> A = [ A11 A12 ] and B = [ B11 B12 ] 128*> [ A22 ] [ B22 ] 129*> 130*> where for different type of A11 and A22 are given as the following. 131*> A12 and B12 are chosen so that the generalized Sylvester equation 132*> 133*> A11*R - L*A22 = -A12 134*> B11*R - L*B22 = -B12 135*> 136*> have prescribed solution R and L. 137*> 138*> Type 1: A11 = J_m(1,-1) and A_22 = J_k(1-a,1). 139*> B11 = I_m, B22 = I_k 140*> where J_k(a,b) is the k-by-k Jordan block with ``a'' on 141*> diagonal and ``b'' on superdiagonal. 142*> 143*> Type 2: A11 = (a_ij) = ( 2(.5-sin(i)) ) and 144*> B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m 145*> A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and 146*> B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k 147*> 148*> Type 3: A11, A22 and B11, B22 are chosen as for Type 2, but each 149*> second diagonal block in A_11 and each third diagonal block 150*> in A_22 are made as 2 by 2 blocks. 151*> 152*> Type 4: A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) ) 153*> for i=1,...,m, j=1,...,m and 154*> A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) ) 155*> for i=m+1,...,k, j=m+1,...,k 156*> 157*> Type 5: (A,B) and have potentially close or common eigenvalues and 158*> very large departure from block diagonality A_11 is chosen 159*> as the m x m leading submatrix of A_1: 160*> | 1 b | 161*> | -b 1 | 162*> | 1+d b | 163*> | -b 1+d | 164*> A_1 = | d 1 | 165*> | -1 d | 166*> | -d 1 | 167*> | -1 -d | 168*> | 1 | 169*> and A_22 is chosen as the k x k leading submatrix of A_2: 170*> | -1 b | 171*> | -b -1 | 172*> | 1-d b | 173*> | -b 1-d | 174*> A_2 = | d 1+b | 175*> | -1-b d | 176*> | -d 1+b | 177*> | -1+b -d | 178*> | 1-d | 179*> and matrix B are chosen as identity matrices (see DLATM5). 180*> 181*> \endverbatim 182* 183* Arguments: 184* ========== 185* 186*> \param[in] NSIZE 187*> \verbatim 188*> NSIZE is INTEGER 189*> The maximum size of the matrices to use. NSIZE >= 0. 190*> If NSIZE = 0, no built-in tests matrices are used, but 191*> read-in test matrices are used to test DGGESX. 192*> \endverbatim 193*> 194*> \param[in] NCMAX 195*> \verbatim 196*> NCMAX is INTEGER 197*> Maximum allowable NMAX for generating Kroneker matrix 198*> in call to DLAKF2 199*> \endverbatim 200*> 201*> \param[in] THRESH 202*> \verbatim 203*> THRESH is DOUBLE PRECISION 204*> A test will count as "failed" if the "error", computed as 205*> described above, exceeds THRESH. Note that the error 206*> is scaled to be O(1), so THRESH should be a reasonably 207*> small multiple of 1, e.g., 10 or 100. In particular, 208*> it should not depend on the precision (single vs. double) 209*> or the size of the matrix. THRESH >= 0. 210*> \endverbatim 211*> 212*> \param[in] NIN 213*> \verbatim 214*> NIN is INTEGER 215*> The FORTRAN unit number for reading in the data file of 216*> problems to solve. 217*> \endverbatim 218*> 219*> \param[in] NOUT 220*> \verbatim 221*> NOUT is INTEGER 222*> The FORTRAN unit number for printing out error messages 223*> (e.g., if a routine returns IINFO not equal to 0.) 224*> \endverbatim 225*> 226*> \param[out] A 227*> \verbatim 228*> A is DOUBLE PRECISION array, dimension (LDA, NSIZE) 229*> Used to store the matrix whose eigenvalues are to be 230*> computed. On exit, A contains the last matrix actually used. 231*> \endverbatim 232*> 233*> \param[in] LDA 234*> \verbatim 235*> LDA is INTEGER 236*> The leading dimension of A, B, AI, BI, Z and Q, 237*> LDA >= max( 1, NSIZE ). For the read-in test, 238*> LDA >= max( 1, N ), N is the size of the test matrices. 239*> \endverbatim 240*> 241*> \param[out] B 242*> \verbatim 243*> B is DOUBLE PRECISION array, dimension (LDA, NSIZE) 244*> Used to store the matrix whose eigenvalues are to be 245*> computed. On exit, B contains the last matrix actually used. 246*> \endverbatim 247*> 248*> \param[out] AI 249*> \verbatim 250*> AI is DOUBLE PRECISION array, dimension (LDA, NSIZE) 251*> Copy of A, modified by DGGESX. 252*> \endverbatim 253*> 254*> \param[out] BI 255*> \verbatim 256*> BI is DOUBLE PRECISION array, dimension (LDA, NSIZE) 257*> Copy of B, modified by DGGESX. 258*> \endverbatim 259*> 260*> \param[out] Z 261*> \verbatim 262*> Z is DOUBLE PRECISION array, dimension (LDA, NSIZE) 263*> Z holds the left Schur vectors computed by DGGESX. 264*> \endverbatim 265*> 266*> \param[out] Q 267*> \verbatim 268*> Q is DOUBLE PRECISION array, dimension (LDA, NSIZE) 269*> Q holds the right Schur vectors computed by DGGESX. 270*> \endverbatim 271*> 272*> \param[out] ALPHAR 273*> \verbatim 274*> ALPHAR is DOUBLE PRECISION array, dimension (NSIZE) 275*> \endverbatim 276*> 277*> \param[out] ALPHAI 278*> \verbatim 279*> ALPHAI is DOUBLE PRECISION array, dimension (NSIZE) 280*> \endverbatim 281*> 282*> \param[out] BETA 283*> \verbatim 284*> BETA is DOUBLE PRECISION array, dimension (NSIZE) 285*> 286*> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues. 287*> \endverbatim 288*> 289*> \param[out] C 290*> \verbatim 291*> C is DOUBLE PRECISION array, dimension (LDC, LDC) 292*> Store the matrix generated by subroutine DLAKF2, this is the 293*> matrix formed by Kronecker products used for estimating 294*> DIF. 295*> \endverbatim 296*> 297*> \param[in] LDC 298*> \verbatim 299*> LDC is INTEGER 300*> The leading dimension of C. LDC >= max(1, LDA*LDA/2 ). 301*> \endverbatim 302*> 303*> \param[out] S 304*> \verbatim 305*> S is DOUBLE PRECISION array, dimension (LDC) 306*> Singular values of C 307*> \endverbatim 308*> 309*> \param[out] WORK 310*> \verbatim 311*> WORK is DOUBLE PRECISION array, dimension (LWORK) 312*> \endverbatim 313*> 314*> \param[in] LWORK 315*> \verbatim 316*> LWORK is INTEGER 317*> The dimension of the array WORK. 318*> LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) ) 319*> \endverbatim 320*> 321*> \param[out] IWORK 322*> \verbatim 323*> IWORK is INTEGER array, dimension (LIWORK) 324*> \endverbatim 325*> 326*> \param[in] LIWORK 327*> \verbatim 328*> LIWORK is INTEGER 329*> The dimension of the array IWORK. LIWORK >= NSIZE + 6. 330*> \endverbatim 331*> 332*> \param[out] BWORK 333*> \verbatim 334*> BWORK is LOGICAL array, dimension (LDA) 335*> \endverbatim 336*> 337*> \param[out] INFO 338*> \verbatim 339*> INFO is INTEGER 340*> = 0: successful exit 341*> < 0: if INFO = -i, the i-th argument had an illegal value. 342*> > 0: A routine returned an error code. 343*> \endverbatim 344* 345* Authors: 346* ======== 347* 348*> \author Univ. of Tennessee 349*> \author Univ. of California Berkeley 350*> \author Univ. of Colorado Denver 351*> \author NAG Ltd. 352* 353*> \ingroup double_eig 354* 355* ===================================================================== 356 SUBROUTINE DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, 357 $ BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S, 358 $ WORK, LWORK, IWORK, LIWORK, BWORK, INFO ) 359* 360* -- LAPACK test routine -- 361* -- LAPACK is a software package provided by Univ. of Tennessee, -- 362* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 363* 364* .. Scalar Arguments .. 365 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN, 366 $ NOUT, NSIZE 367 DOUBLE PRECISION THRESH 368* .. 369* .. Array Arguments .. 370 LOGICAL BWORK( * ) 371 INTEGER IWORK( * ) 372 DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ), 373 $ ALPHAR( * ), B( LDA, * ), BETA( * ), 374 $ BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ), 375 $ WORK( * ), Z( LDA, * ) 376* .. 377* 378* ===================================================================== 379* 380* .. Parameters .. 381 DOUBLE PRECISION ZERO, ONE, TEN 382 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1 ) 383* .. 384* .. Local Scalars .. 385 LOGICAL ILABAD 386 CHARACTER SENSE 387 INTEGER BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK, 388 $ MINWRK, MM, MN2, NERRS, NPTKNT, NTEST, NTESTT, 389 $ PRTYPE, QBA, QBB 390 DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1, 391 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT 392* .. 393* .. Local Arrays .. 394 DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 ) 395* .. 396* .. External Functions .. 397 LOGICAL DLCTSX 398 INTEGER ILAENV 399 DOUBLE PRECISION DLAMCH, DLANGE 400 EXTERNAL DLCTSX, ILAENV, DLAMCH, DLANGE 401* .. 402* .. External Subroutines .. 403 EXTERNAL ALASVM, DGESVD, DGET51, DGET53, DGGESX, DLABAD, 404 $ DLACPY, DLAKF2, DLASET, DLATM5, XERBLA 405* .. 406* .. Intrinsic Functions .. 407 INTRINSIC ABS, MAX, SQRT 408* .. 409* .. Scalars in Common .. 410 LOGICAL FS 411 INTEGER K, M, MPLUSN, N 412* .. 413* .. Common blocks .. 414 COMMON / MN / M, N, MPLUSN, K, FS 415* .. 416* .. Executable Statements .. 417* 418* Check for errors 419* 420 IF( NSIZE.LT.0 ) THEN 421 INFO = -1 422 ELSE IF( THRESH.LT.ZERO ) THEN 423 INFO = -2 424 ELSE IF( NIN.LE.0 ) THEN 425 INFO = -3 426 ELSE IF( NOUT.LE.0 ) THEN 427 INFO = -4 428 ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN 429 INFO = -6 430 ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN 431 INFO = -17 432 ELSE IF( LIWORK.LT.NSIZE+6 ) THEN 433 INFO = -21 434 END IF 435* 436* Compute workspace 437* (Note: Comments in the code beginning "Workspace:" describe the 438* minimal amount of workspace needed at that point in the code, 439* as well as the preferred amount for good performance. 440* NB refers to the optimal block size for the immediately 441* following subroutine, as returned by ILAENV.) 442* 443 MINWRK = 1 444 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 445 MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 ) 446* 447* workspace for sggesx 448* 449 MAXWRK = 9*( NSIZE+1 ) + NSIZE* 450 $ ILAENV( 1, 'DGEQRF', ' ', NSIZE, 1, NSIZE, 0 ) 451 MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE* 452 $ ILAENV( 1, 'DORGQR', ' ', NSIZE, 1, NSIZE, -1 ) ) 453* 454* workspace for dgesvd 455* 456 BDSPAC = 5*NSIZE*NSIZE / 2 457 MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE* 458 $ ILAENV( 1, 'DGEBRD', ' ', NSIZE*NSIZE / 2, 459 $ NSIZE*NSIZE / 2, -1, -1 ) ) 460 MAXWRK = MAX( MAXWRK, BDSPAC ) 461* 462 MAXWRK = MAX( MAXWRK, MINWRK ) 463* 464 WORK( 1 ) = MAXWRK 465 END IF 466* 467 IF( LWORK.LT.MINWRK ) 468 $ INFO = -19 469* 470 IF( INFO.NE.0 ) THEN 471 CALL XERBLA( 'DDRGSX', -INFO ) 472 RETURN 473 END IF 474* 475* Important constants 476* 477 ULP = DLAMCH( 'P' ) 478 ULPINV = ONE / ULP 479 SMLNUM = DLAMCH( 'S' ) / ULP 480 BIGNUM = ONE / SMLNUM 481 CALL DLABAD( SMLNUM, BIGNUM ) 482 THRSH2 = TEN*THRESH 483 NTESTT = 0 484 NERRS = 0 485* 486* Go to the tests for read-in matrix pairs 487* 488 IFUNC = 0 489 IF( NSIZE.EQ.0 ) 490 $ GO TO 70 491* 492* Test the built-in matrix pairs. 493* Loop over different functions (IFUNC) of DGGESX, types (PRTYPE) 494* of test matrices, different size (M+N) 495* 496 PRTYPE = 0 497 QBA = 3 498 QBB = 4 499 WEIGHT = SQRT( ULP ) 500* 501 DO 60 IFUNC = 0, 3 502 DO 50 PRTYPE = 1, 5 503 DO 40 M = 1, NSIZE - 1 504 DO 30 N = 1, NSIZE - M 505* 506 WEIGHT = ONE / WEIGHT 507 MPLUSN = M + N 508* 509* Generate test matrices 510* 511 FS = .TRUE. 512 K = 0 513* 514 CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI, 515 $ LDA ) 516 CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI, 517 $ LDA ) 518* 519 CALL DLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ), 520 $ LDA, AI( 1, M+1 ), LDA, BI, LDA, 521 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA, 522 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB ) 523* 524* Compute the Schur factorization and swapping the 525* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 526* Swapping is accomplished via the function DLCTSX 527* which is supplied below. 528* 529 IF( IFUNC.EQ.0 ) THEN 530 SENSE = 'N' 531 ELSE IF( IFUNC.EQ.1 ) THEN 532 SENSE = 'E' 533 ELSE IF( IFUNC.EQ.2 ) THEN 534 SENSE = 'V' 535 ELSE IF( IFUNC.EQ.3 ) THEN 536 SENSE = 'B' 537 END IF 538* 539 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 540 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 541* 542 CALL DGGESX( 'V', 'V', 'S', DLCTSX, SENSE, MPLUSN, AI, 543 $ LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA, 544 $ Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK, 545 $ IWORK, LIWORK, BWORK, LINFO ) 546* 547 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 548 RESULT( 1 ) = ULPINV 549 WRITE( NOUT, FMT = 9999 )'DGGESX', LINFO, MPLUSN, 550 $ PRTYPE 551 INFO = LINFO 552 GO TO 30 553 END IF 554* 555* Compute the norm(A, B) 556* 557 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, 558 $ MPLUSN ) 559 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 560 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 561 ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, 562 $ WORK ) 563* 564* Do tests (1) to (4) 565* 566 CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, 567 $ LDA, WORK, RESULT( 1 ) ) 568 CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, 569 $ LDA, WORK, RESULT( 2 ) ) 570 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, 571 $ LDA, WORK, RESULT( 3 ) ) 572 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, 573 $ LDA, WORK, RESULT( 4 ) ) 574 NTEST = 4 575* 576* Do tests (5) and (6): check Schur form of A and 577* compare eigenvalues with diagonals. 578* 579 TEMP1 = ZERO 580 RESULT( 5 ) = ZERO 581 RESULT( 6 ) = ZERO 582* 583 DO 10 J = 1, MPLUSN 584 ILABAD = .FALSE. 585 IF( ALPHAI( J ).EQ.ZERO ) THEN 586 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) / 587 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), 588 $ ABS( AI( J, J ) ) )+ 589 $ ABS( BETA( J )-BI( J, J ) ) / 590 $ MAX( SMLNUM, ABS( BETA( J ) ), 591 $ ABS( BI( J, J ) ) ) ) / ULP 592 IF( J.LT.MPLUSN ) THEN 593 IF( AI( J+1, J ).NE.ZERO ) THEN 594 ILABAD = .TRUE. 595 RESULT( 5 ) = ULPINV 596 END IF 597 END IF 598 IF( J.GT.1 ) THEN 599 IF( AI( J, J-1 ).NE.ZERO ) THEN 600 ILABAD = .TRUE. 601 RESULT( 5 ) = ULPINV 602 END IF 603 END IF 604 ELSE 605 IF( ALPHAI( J ).GT.ZERO ) THEN 606 I1 = J 607 ELSE 608 I1 = J - 1 609 END IF 610 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN 611 ILABAD = .TRUE. 612 ELSE IF( I1.LT.MPLUSN-1 ) THEN 613 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN 614 ILABAD = .TRUE. 615 RESULT( 5 ) = ULPINV 616 END IF 617 ELSE IF( I1.GT.1 ) THEN 618 IF( AI( I1, I1-1 ).NE.ZERO ) THEN 619 ILABAD = .TRUE. 620 RESULT( 5 ) = ULPINV 621 END IF 622 END IF 623 IF( .NOT.ILABAD ) THEN 624 CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), 625 $ LDA, BETA( J ), ALPHAR( J ), 626 $ ALPHAI( J ), TEMP2, IINFO ) 627 IF( IINFO.GE.3 ) THEN 628 WRITE( NOUT, FMT = 9997 )IINFO, J, 629 $ MPLUSN, PRTYPE 630 INFO = ABS( IINFO ) 631 END IF 632 ELSE 633 TEMP2 = ULPINV 634 END IF 635 END IF 636 TEMP1 = MAX( TEMP1, TEMP2 ) 637 IF( ILABAD ) THEN 638 WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE 639 END IF 640 10 CONTINUE 641 RESULT( 6 ) = TEMP1 642 NTEST = NTEST + 2 643* 644* Test (7) (if sorting worked) 645* 646 RESULT( 7 ) = ZERO 647 IF( LINFO.EQ.MPLUSN+3 ) THEN 648 RESULT( 7 ) = ULPINV 649 ELSE IF( MM.NE.N ) THEN 650 RESULT( 7 ) = ULPINV 651 END IF 652 NTEST = NTEST + 1 653* 654* Test (8): compare the estimated value DIF and its 655* value. first, compute the exact DIF. 656* 657 RESULT( 8 ) = ZERO 658 MN2 = MM*( MPLUSN-MM )*2 659 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN 660* 661* Note: for either following two causes, there are 662* almost same number of test cases fail the test. 663* 664 CALL DLAKF2( MM, MPLUSN-MM, AI, LDA, 665 $ AI( MM+1, MM+1 ), BI, 666 $ BI( MM+1, MM+1 ), C, LDC ) 667* 668 CALL DGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK, 669 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2, 670 $ INFO ) 671 DIFTRU = S( MN2 ) 672* 673 IF( DIFEST( 2 ).EQ.ZERO ) THEN 674 IF( DIFTRU.GT.ABNRM*ULP ) 675 $ RESULT( 8 ) = ULPINV 676 ELSE IF( DIFTRU.EQ.ZERO ) THEN 677 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 678 $ RESULT( 8 ) = ULPINV 679 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 680 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 681 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), 682 $ DIFEST( 2 ) / DIFTRU ) 683 END IF 684 NTEST = NTEST + 1 685 END IF 686* 687* Test (9) 688* 689 RESULT( 9 ) = ZERO 690 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 691 IF( DIFTRU.GT.ABNRM*ULP ) 692 $ RESULT( 9 ) = ULPINV 693 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 694 $ RESULT( 9 ) = ULPINV 695 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 696 $ RESULT( 9 ) = ULPINV 697 NTEST = NTEST + 1 698 END IF 699* 700 NTESTT = NTESTT + NTEST 701* 702* Print out tests which fail. 703* 704 DO 20 J = 1, 9 705 IF( RESULT( J ).GE.THRESH ) THEN 706* 707* If this is the first test to fail, 708* print a header to the data file. 709* 710 IF( NERRS.EQ.0 ) THEN 711 WRITE( NOUT, FMT = 9995 )'DGX' 712* 713* Matrix types 714* 715 WRITE( NOUT, FMT = 9993 ) 716* 717* Tests performed 718* 719 WRITE( NOUT, FMT = 9992 )'orthogonal', '''', 720 $ 'transpose', ( '''', I = 1, 4 ) 721* 722 END IF 723 NERRS = NERRS + 1 724 IF( RESULT( J ).LT.10000.0D0 ) THEN 725 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE, 726 $ WEIGHT, M, J, RESULT( J ) 727 ELSE 728 WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE, 729 $ WEIGHT, M, J, RESULT( J ) 730 END IF 731 END IF 732 20 CONTINUE 733* 734 30 CONTINUE 735 40 CONTINUE 736 50 CONTINUE 737 60 CONTINUE 738* 739 GO TO 150 740* 741 70 CONTINUE 742* 743* Read in data from file to check accuracy of condition estimation 744* Read input data until N=0 745* 746 NPTKNT = 0 747* 748 80 CONTINUE 749 READ( NIN, FMT = *, END = 140 )MPLUSN 750 IF( MPLUSN.EQ.0 ) 751 $ GO TO 140 752 READ( NIN, FMT = *, END = 140 )N 753 DO 90 I = 1, MPLUSN 754 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN ) 755 90 CONTINUE 756 DO 100 I = 1, MPLUSN 757 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN ) 758 100 CONTINUE 759 READ( NIN, FMT = * )PLTRU, DIFTRU 760* 761 NPTKNT = NPTKNT + 1 762 FS = .TRUE. 763 K = 0 764 M = MPLUSN - N 765* 766 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 767 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 768* 769* Compute the Schur factorization while swapping the 770* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 771* 772 CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA, 773 $ MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST, 774 $ WORK, LWORK, IWORK, LIWORK, BWORK, LINFO ) 775* 776 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 777 RESULT( 1 ) = ULPINV 778 WRITE( NOUT, FMT = 9998 )'DGGESX', LINFO, MPLUSN, NPTKNT 779 GO TO 130 780 END IF 781* 782* Compute the norm(A, B) 783* (should this be norm of (A,B) or (AI,BI)?) 784* 785 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN ) 786 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 787 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 788 ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK ) 789* 790* Do tests (1) to (4) 791* 792 CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK, 793 $ RESULT( 1 ) ) 794 CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK, 795 $ RESULT( 2 ) ) 796 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK, 797 $ RESULT( 3 ) ) 798 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK, 799 $ RESULT( 4 ) ) 800* 801* Do tests (5) and (6): check Schur form of A and compare 802* eigenvalues with diagonals. 803* 804 NTEST = 6 805 TEMP1 = ZERO 806 RESULT( 5 ) = ZERO 807 RESULT( 6 ) = ZERO 808* 809 DO 110 J = 1, MPLUSN 810 ILABAD = .FALSE. 811 IF( ALPHAI( J ).EQ.ZERO ) THEN 812 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) / 813 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J, 814 $ J ) ) )+ABS( BETA( J )-BI( J, J ) ) / 815 $ MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) ) 816 $ / ULP 817 IF( J.LT.MPLUSN ) THEN 818 IF( AI( J+1, J ).NE.ZERO ) THEN 819 ILABAD = .TRUE. 820 RESULT( 5 ) = ULPINV 821 END IF 822 END IF 823 IF( J.GT.1 ) THEN 824 IF( AI( J, J-1 ).NE.ZERO ) THEN 825 ILABAD = .TRUE. 826 RESULT( 5 ) = ULPINV 827 END IF 828 END IF 829 ELSE 830 IF( ALPHAI( J ).GT.ZERO ) THEN 831 I1 = J 832 ELSE 833 I1 = J - 1 834 END IF 835 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN 836 ILABAD = .TRUE. 837 ELSE IF( I1.LT.MPLUSN-1 ) THEN 838 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN 839 ILABAD = .TRUE. 840 RESULT( 5 ) = ULPINV 841 END IF 842 ELSE IF( I1.GT.1 ) THEN 843 IF( AI( I1, I1-1 ).NE.ZERO ) THEN 844 ILABAD = .TRUE. 845 RESULT( 5 ) = ULPINV 846 END IF 847 END IF 848 IF( .NOT.ILABAD ) THEN 849 CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA, 850 $ BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2, 851 $ IINFO ) 852 IF( IINFO.GE.3 ) THEN 853 WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT 854 INFO = ABS( IINFO ) 855 END IF 856 ELSE 857 TEMP2 = ULPINV 858 END IF 859 END IF 860 TEMP1 = MAX( TEMP1, TEMP2 ) 861 IF( ILABAD ) THEN 862 WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT 863 END IF 864 110 CONTINUE 865 RESULT( 6 ) = TEMP1 866* 867* Test (7) (if sorting worked) <--------- need to be checked. 868* 869 NTEST = 7 870 RESULT( 7 ) = ZERO 871 IF( LINFO.EQ.MPLUSN+3 ) 872 $ RESULT( 7 ) = ULPINV 873* 874* Test (8): compare the estimated value of DIF and its true value. 875* 876 NTEST = 8 877 RESULT( 8 ) = ZERO 878 IF( DIFEST( 2 ).EQ.ZERO ) THEN 879 IF( DIFTRU.GT.ABNRM*ULP ) 880 $ RESULT( 8 ) = ULPINV 881 ELSE IF( DIFTRU.EQ.ZERO ) THEN 882 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 883 $ RESULT( 8 ) = ULPINV 884 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 885 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 886 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU ) 887 END IF 888* 889* Test (9) 890* 891 NTEST = 9 892 RESULT( 9 ) = ZERO 893 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 894 IF( DIFTRU.GT.ABNRM*ULP ) 895 $ RESULT( 9 ) = ULPINV 896 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 897 $ RESULT( 9 ) = ULPINV 898 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 899 $ RESULT( 9 ) = ULPINV 900 END IF 901* 902* Test (10): compare the estimated value of PL and it true value. 903* 904 NTEST = 10 905 RESULT( 10 ) = ZERO 906 IF( PL( 1 ).EQ.ZERO ) THEN 907 IF( PLTRU.GT.ABNRM*ULP ) 908 $ RESULT( 10 ) = ULPINV 909 ELSE IF( PLTRU.EQ.ZERO ) THEN 910 IF( PL( 1 ).GT.ABNRM*ULP ) 911 $ RESULT( 10 ) = ULPINV 912 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR. 913 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN 914 RESULT( 10 ) = ULPINV 915 END IF 916* 917 NTESTT = NTESTT + NTEST 918* 919* Print out tests which fail. 920* 921 DO 120 J = 1, NTEST 922 IF( RESULT( J ).GE.THRESH ) THEN 923* 924* If this is the first test to fail, 925* print a header to the data file. 926* 927 IF( NERRS.EQ.0 ) THEN 928 WRITE( NOUT, FMT = 9995 )'DGX' 929* 930* Matrix types 931* 932 WRITE( NOUT, FMT = 9994 ) 933* 934* Tests performed 935* 936 WRITE( NOUT, FMT = 9992 )'orthogonal', '''', 937 $ 'transpose', ( '''', I = 1, 4 ) 938* 939 END IF 940 NERRS = NERRS + 1 941 IF( RESULT( J ).LT.10000.0D0 ) THEN 942 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J ) 943 ELSE 944 WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J ) 945 END IF 946 END IF 947* 948 120 CONTINUE 949* 950 130 CONTINUE 951 GO TO 80 952 140 CONTINUE 953* 954 150 CONTINUE 955* 956* Summary 957* 958 CALL ALASVM( 'DGX', NOUT, NERRS, NTESTT, 0 ) 959* 960 WORK( 1 ) = MAXWRK 961* 962 RETURN 963* 964 9999 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 965 $ I6, ', JTYPE=', I6, ')' ) 966* 967 9998 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 968 $ I6, ', Input Example #', I2, ')' ) 969* 970 9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', I1, ' for eigenvalue ', 971 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' ) 972* 973 9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', I6, '.', 974 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' ) 975* 976 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form', 977 $ ' problem driver' ) 978* 979 9994 FORMAT( 'Input Example' ) 980* 981 9993 FORMAT( ' Matrix types: ', / 982 $ ' 1: A is a block diagonal matrix of Jordan blocks ', 983 $ 'and B is the identity ', / ' matrix, ', 984 $ / ' 2: A and B are upper triangular matrices, ', 985 $ / ' 3: A and B are as type 2, but each second diagonal ', 986 $ 'block in A_11 and ', / 987 $ ' each third diaongal block in A_22 are 2x2 blocks,', 988 $ / ' 4: A and B are block diagonal matrices, ', 989 $ / ' 5: (A,B) has potentially close or common ', 990 $ 'eigenvalues.', / ) 991* 992 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ', 993 $ 'Q and Z are ', A, ',', / 19X, 994 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)', 995 $ / ' 1 = | A - Q S Z', A, 996 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A, 997 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A, 998 $ ' | / ( n ulp ) 4 = | I - ZZ', A, 999 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ', 1000 $ 'Schur form S', / ' 6 = difference between (alpha,beta)', 1001 $ ' and diagonals of (S,T)', / 1002 $ ' 7 = 1/ULP if SDIM is not the correct number of ', 1003 $ 'selected eigenvalues', / 1004 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ', 1005 $ 'DIFTRU/DIFEST > 10*THRESH', 1006 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ', 1007 $ 'when reordering fails', / 1008 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ', 1009 $ 'PLTRU/PLEST > THRESH', / 1010 $ ' ( Test 10 is only for input examples )', / ) 1011 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3, 1012 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 ) 1013 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3, 1014 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.3 ) 1015 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 1016 $ ' result ', I2, ' is', 0P, F8.2 ) 1017 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 1018 $ ' result ', I2, ' is', 1P, D10.3 ) 1019* 1020* End of DDRGSX 1021* 1022 END 1023