 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, 00002 $ BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK, 00003 $ RWORK, IWORK, LIWORK, BWORK, INFO ) 00004 * 00005 * -- LAPACK test routine (version 3.3.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN, 00011 $ NOUT, NSIZE 00012 DOUBLE PRECISION THRESH 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL BWORK( * ) 00016 INTEGER IWORK( * ) 00017 DOUBLE PRECISION RWORK( * ), S( * ) 00018 COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ), 00019 $ B( LDA, * ), BETA( * ), BI( LDA, * ), 00020 $ C( LDC, * ), Q( LDA, * ), WORK( * ), 00021 $ Z( LDA, * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * ZDRGSX checks the nonsymmetric generalized eigenvalue (Schur form) 00028 * problem expert driver ZGGESX. 00029 * 00030 * ZGGES factors A and B as Q*S*Z' and Q*T*Z' , where ' means conjugate 00031 * transpose, S and T are upper triangular (i.e., in generalized Schur 00032 * form), and Q and Z are unitary. It also computes the generalized 00033 * eigenvalues (alpha(j),beta(j)), j=1,...,n. Thus, 00034 * w(j) = alpha(j)/beta(j) is a root of the characteristic equation 00035 * 00036 * det( A - w(j) B ) = 0 00037 * 00038 * Optionally it also reorders the eigenvalues so that a selected 00039 * cluster of eigenvalues appears in the leading diagonal block of the 00040 * Schur forms; computes a reciprocal condition number for the average 00041 * of the selected eigenvalues; and computes a reciprocal condition 00042 * number for the right and left deflating subspaces corresponding to 00043 * the selected eigenvalues. 00044 * 00045 * When ZDRGSX is called with NSIZE > 0, five (5) types of built-in 00046 * matrix pairs are used to test the routine ZGGESX. 00047 * 00048 * When ZDRGSX is called with NSIZE = 0, it reads in test matrix data 00049 * to test ZGGESX. 00050 * (need more details on what kind of read-in data are needed). 00051 * 00052 * For each matrix pair, the following tests will be performed and 00053 * compared with the threshhold THRESH except for the tests (7) and (9): 00054 * 00055 * (1) | A - Q S Z' | / ( |A| n ulp ) 00056 * 00057 * (2) | B - Q T Z' | / ( |B| n ulp ) 00058 * 00059 * (3) | I - QQ' | / ( n ulp ) 00060 * 00061 * (4) | I - ZZ' | / ( n ulp ) 00062 * 00063 * (5) if A is in Schur form (i.e. triangular form) 00064 * 00065 * (6) maximum over j of D(j) where: 00066 * 00067 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)| 00068 * D(j) = ------------------------ + ----------------------- 00069 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|) 00070 * 00071 * (7) if sorting worked and SDIM is the number of eigenvalues 00072 * which were selected. 00073 * 00074 * (8) the estimated value DIF does not differ from the true values of 00075 * Difu and Difl more than a factor 10*THRESH. If the estimate DIF 00076 * equals zero the corresponding true values of Difu and Difl 00077 * should be less than EPS*norm(A, B). If the true value of Difu 00078 * and Difl equal zero, the estimate DIF should be less than 00079 * EPS*norm(A, B). 00080 * 00081 * (9) If INFO = N+3 is returned by ZGGESX, the reordering "failed" 00082 * and we check that DIF = PL = PR = 0 and that the true value of 00083 * Difu and Difl is < EPS*norm(A, B). We count the events when 00084 * INFO=N+3. 00085 * 00086 * For read-in test matrices, the same tests are run except that the 00087 * exact value for DIF (and PL) is input data. Additionally, there is 00088 * one more test run for read-in test matrices: 00089 * 00090 * (10) the estimated value PL does not differ from the true value of 00091 * PLTRU more than a factor THRESH. If the estimate PL equals 00092 * zero the corresponding true value of PLTRU should be less than 00093 * EPS*norm(A, B). If the true value of PLTRU equal zero, the 00094 * estimate PL should be less than EPS*norm(A, B). 00095 * 00096 * Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1) 00097 * matrix pairs are generated and tested. NSIZE should be kept small. 00098 * 00099 * SVD (routine ZGESVD) is used for computing the true value of DIF_u 00100 * and DIF_l when testing the built-in test problems. 00101 * 00102 * Built-in Test Matrices 00103 * ====================== 00104 * 00105 * All built-in test matrices are the 2 by 2 block of triangular 00106 * matrices 00107 * 00108 * A = [ A11 A12 ] and B = [ B11 B12 ] 00109 * [ A22 ] [ B22 ] 00110 * 00111 * where for different type of A11 and A22 are given as the following. 00112 * A12 and B12 are chosen so that the generalized Sylvester equation 00113 * 00114 * A11*R - L*A22 = -A12 00115 * B11*R - L*B22 = -B12 00116 * 00117 * have prescribed solution R and L. 00118 * 00119 * Type 1: A11 = J_m(1,-1) and A_22 = J_k(1-a,1). 00120 * B11 = I_m, B22 = I_k 00121 * where J_k(a,b) is the k-by-k Jordan block with ``a'' on 00122 * diagonal and ``b'' on superdiagonal. 00123 * 00124 * Type 2: A11 = (a_ij) = ( 2(.5-sin(i)) ) and 00125 * B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m 00126 * A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and 00127 * B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k 00128 * 00129 * Type 3: A11, A22 and B11, B22 are chosen as for Type 2, but each 00130 * second diagonal block in A_11 and each third diagonal block 00131 * in A_22 are made as 2 by 2 blocks. 00132 * 00133 * Type 4: A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) ) 00134 * for i=1,...,m, j=1,...,m and 00135 * A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) ) 00136 * for i=m+1,...,k, j=m+1,...,k 00137 * 00138 * Type 5: (A,B) and have potentially close or common eigenvalues and 00139 * very large departure from block diagonality A_11 is chosen 00140 * as the m x m leading submatrix of A_1: 00141 * | 1 b | 00142 * | -b 1 | 00143 * | 1+d b | 00144 * | -b 1+d | 00145 * A_1 = | d 1 | 00146 * | -1 d | 00147 * | -d 1 | 00148 * | -1 -d | 00149 * | 1 | 00150 * and A_22 is chosen as the k x k leading submatrix of A_2: 00151 * | -1 b | 00152 * | -b -1 | 00153 * | 1-d b | 00154 * | -b 1-d | 00155 * A_2 = | d 1+b | 00156 * | -1-b d | 00157 * | -d 1+b | 00158 * | -1+b -d | 00159 * | 1-d | 00160 * and matrix B are chosen as identity matrices (see DLATM5). 00161 * 00162 * 00163 * Arguments 00164 * ========= 00165 * 00166 * NSIZE (input) INTEGER 00167 * The maximum size of the matrices to use. NSIZE >= 0. 00168 * If NSIZE = 0, no built-in tests matrices are used, but 00169 * read-in test matrices are used to test DGGESX. 00170 * 00171 * NCMAX (input) INTEGER 00172 * Maximum allowable NMAX for generating Kroneker matrix 00173 * in call to ZLAKF2 00174 * 00175 * THRESH (input) DOUBLE PRECISION 00176 * A test will count as "failed" if the "error", computed as 00177 * described above, exceeds THRESH. Note that the error 00178 * is scaled to be O(1), so THRESH should be a reasonably 00179 * small multiple of 1, e.g., 10 or 100. In particular, 00180 * it should not depend on the precision (single vs. double) 00181 * or the size of the matrix. THRESH >= 0. 00182 * 00183 * NIN (input) INTEGER 00184 * The FORTRAN unit number for reading in the data file of 00185 * problems to solve. 00186 * 00187 * NOUT (input) INTEGER 00188 * The FORTRAN unit number for printing out error messages 00189 * (e.g., if a routine returns INFO not equal to 0.) 00190 * 00191 * A (workspace) COMPLEX*16 array, dimension (LDA, NSIZE) 00192 * Used to store the matrix whose eigenvalues are to be 00193 * computed. On exit, A contains the last matrix actually used. 00194 * 00195 * LDA (input) INTEGER 00196 * The leading dimension of A, B, AI, BI, Z and Q, 00197 * LDA >= max( 1, NSIZE ). For the read-in test, 00198 * LDA >= max( 1, N ), N is the size of the test matrices. 00199 * 00200 * B (workspace) COMPLEX*16 array, dimension (LDA, NSIZE) 00201 * Used to store the matrix whose eigenvalues are to be 00202 * computed. On exit, B contains the last matrix actually used. 00203 * 00204 * AI (workspace) COMPLEX*16 array, dimension (LDA, NSIZE) 00205 * Copy of A, modified by ZGGESX. 00206 * 00207 * BI (workspace) COMPLEX*16 array, dimension (LDA, NSIZE) 00208 * Copy of B, modified by ZGGESX. 00209 * 00210 * Z (workspace) COMPLEX*16 array, dimension (LDA, NSIZE) 00211 * Z holds the left Schur vectors computed by ZGGESX. 00212 * 00213 * Q (workspace) COMPLEX*16 array, dimension (LDA, NSIZE) 00214 * Q holds the right Schur vectors computed by ZGGESX. 00215 * 00216 * ALPHA (workspace) COMPLEX*16 array, dimension (NSIZE) 00217 * BETA (workspace) COMPLEX*16 array, dimension (NSIZE) 00218 * On exit, ALPHA/BETA are the eigenvalues. 00219 * 00220 * C (workspace) COMPLEX*16 array, dimension (LDC, LDC) 00221 * Store the matrix generated by subroutine ZLAKF2, this is the 00222 * matrix formed by Kronecker products used for estimating 00223 * DIF. 00224 * 00225 * LDC (input) INTEGER 00226 * The leading dimension of C. LDC >= max(1, LDA*LDA/2 ). 00227 * 00228 * S (workspace) DOUBLE PRECISION array, dimension (LDC) 00229 * Singular values of C 00230 * 00231 * WORK (workspace) COMPLEX*16 array, dimension (LWORK) 00232 * 00233 * LWORK (input) INTEGER 00234 * The dimension of the array WORK. LWORK >= 3*NSIZE*NSIZE/2 00235 * 00236 * RWORK (workspace) DOUBLE PRECISION array, 00237 * dimension (5*NSIZE*NSIZE/2 - 4) 00238 * 00239 * IWORK (workspace) INTEGER array, dimension (LIWORK) 00240 * 00241 * LIWORK (input) INTEGER 00242 * The dimension of the array IWORK. LIWORK >= NSIZE + 2. 00243 * 00244 * BWORK (workspace) LOGICAL array, dimension (NSIZE) 00245 * 00246 * INFO (output) INTEGER 00247 * = 0: successful exit 00248 * < 0: if INFO = -i, the i-th argument had an illegal value. 00249 * > 0: A routine returned an error code. 00250 * 00251 * ===================================================================== 00252 * 00253 * .. Parameters .. 00254 DOUBLE PRECISION ZERO, ONE, TEN 00255 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1 ) 00256 COMPLEX*16 CZERO 00257 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 00258 * .. 00259 * .. Local Scalars .. 00260 LOGICAL ILABAD 00261 CHARACTER SENSE 00262 INTEGER BDSPAC, I, IFUNC, J, LINFO, MAXWRK, MINWRK, MM, 00263 $ MN2, NERRS, NPTKNT, NTEST, NTESTT, PRTYPE, QBA, 00264 $ QBB 00265 DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1, 00266 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT 00267 COMPLEX*16 X 00268 * .. 00269 * .. Local Arrays .. 00270 DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 ) 00271 * .. 00272 * .. External Functions .. 00273 LOGICAL ZLCTSX 00274 INTEGER ILAENV 00275 DOUBLE PRECISION DLAMCH, ZLANGE 00276 EXTERNAL ZLCTSX, ILAENV, DLAMCH, ZLANGE 00277 * .. 00278 * .. External Subroutines .. 00279 EXTERNAL ALASVM, DLABAD, XERBLA, ZGESVD, ZGET51, ZGGESX, 00280 $ ZLACPY, ZLAKF2, ZLASET, ZLATM5 00281 * .. 00282 * .. Scalars in Common .. 00283 LOGICAL FS 00284 INTEGER K, M, MPLUSN, N 00285 * .. 00286 * .. Common blocks .. 00287 COMMON / MN / M, N, MPLUSN, K, FS 00288 * .. 00289 * .. Intrinsic Functions .. 00290 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT 00291 * .. 00292 * .. Statement Functions .. 00293 DOUBLE PRECISION ABS1 00294 * .. 00295 * .. Statement Function definitions .. 00296 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) 00297 * .. 00298 * .. Executable Statements .. 00299 * 00300 * Check for errors 00301 * 00302 INFO = 0 00303 IF( NSIZE.LT.0 ) THEN 00304 INFO = -1 00305 ELSE IF( THRESH.LT.ZERO ) THEN 00306 INFO = -2 00307 ELSE IF( NIN.LE.0 ) THEN 00308 INFO = -3 00309 ELSE IF( NOUT.LE.0 ) THEN 00310 INFO = -4 00311 ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN 00312 INFO = -6 00313 ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN 00314 INFO = -15 00315 ELSE IF( LIWORK.LT.NSIZE+2 ) THEN 00316 INFO = -21 00317 END IF 00318 * 00319 * Compute workspace 00320 * (Note: Comments in the code beginning "Workspace:" describe the 00321 * minimal amount of workspace needed at that point in the code, 00322 * as well as the preferred amount for good performance. 00323 * NB refers to the optimal block size for the immediately 00324 * following subroutine, as returned by ILAENV.) 00325 * 00326 MINWRK = 1 00327 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 00328 MINWRK = 3*NSIZE*NSIZE / 2 00329 * 00330 * workspace for cggesx 00331 * 00332 MAXWRK = NSIZE*( 1+ILAENV( 1, 'ZGEQRF', ' ', NSIZE, 1, NSIZE, 00333 $ 0 ) ) 00334 MAXWRK = MAX( MAXWRK, NSIZE*( 1+ILAENV( 1, 'ZUNGQR', ' ', 00335 $ NSIZE, 1, NSIZE, -1 ) ) ) 00336 * 00337 * workspace for zgesvd 00338 * 00339 BDSPAC = 3*NSIZE*NSIZE / 2 00340 MAXWRK = MAX( MAXWRK, NSIZE*NSIZE* 00341 $ ( 1+ILAENV( 1, 'ZGEBRD', ' ', NSIZE*NSIZE / 2, 00342 $ NSIZE*NSIZE / 2, -1, -1 ) ) ) 00343 MAXWRK = MAX( MAXWRK, BDSPAC ) 00344 * 00345 MAXWRK = MAX( MAXWRK, MINWRK ) 00346 * 00347 WORK( 1 ) = MAXWRK 00348 END IF 00349 * 00350 IF( LWORK.LT.MINWRK ) 00351 $ INFO = -18 00352 * 00353 IF( INFO.NE.0 ) THEN 00354 CALL XERBLA( 'ZDRGSX', -INFO ) 00355 RETURN 00356 END IF 00357 * 00358 * Important constants 00359 * 00360 ULP = DLAMCH( 'P' ) 00361 ULPINV = ONE / ULP 00362 SMLNUM = DLAMCH( 'S' ) / ULP 00363 BIGNUM = ONE / SMLNUM 00364 CALL DLABAD( SMLNUM, BIGNUM ) 00365 THRSH2 = TEN*THRESH 00366 NTESTT = 0 00367 NERRS = 0 00368 * 00369 * Go to the tests for read-in matrix pairs 00370 * 00371 IFUNC = 0 00372 IF( NSIZE.EQ.0 ) 00373 $ GO TO 70 00374 * 00375 * Test the built-in matrix pairs. 00376 * Loop over different functions (IFUNC) of ZGGESX, types (PRTYPE) 00377 * of test matrices, different size (M+N) 00378 * 00379 PRTYPE = 0 00380 QBA = 3 00381 QBB = 4 00382 WEIGHT = SQRT( ULP ) 00383 * 00384 DO 60 IFUNC = 0, 3 00385 DO 50 PRTYPE = 1, 5 00386 DO 40 M = 1, NSIZE - 1 00387 DO 30 N = 1, NSIZE - M 00388 * 00389 WEIGHT = ONE / WEIGHT 00390 MPLUSN = M + N 00391 * 00392 * Generate test matrices 00393 * 00394 FS = .TRUE. 00395 K = 0 00396 * 00397 CALL ZLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, AI, 00398 $ LDA ) 00399 CALL ZLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, BI, 00400 $ LDA ) 00401 * 00402 CALL ZLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ), 00403 $ LDA, AI( 1, M+1 ), LDA, BI, LDA, 00404 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA, 00405 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB ) 00406 * 00407 * Compute the Schur factorization and swapping the 00408 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 00409 * Swapping is accomplished via the function ZLCTSX 00410 * which is supplied below. 00411 * 00412 IF( IFUNC.EQ.0 ) THEN 00413 SENSE = 'N' 00414 ELSE IF( IFUNC.EQ.1 ) THEN 00415 SENSE = 'E' 00416 ELSE IF( IFUNC.EQ.2 ) THEN 00417 SENSE = 'V' 00418 ELSE IF( IFUNC.EQ.3 ) THEN 00419 SENSE = 'B' 00420 END IF 00421 * 00422 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00423 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00424 * 00425 CALL ZGGESX( 'V', 'V', 'S', ZLCTSX, SENSE, MPLUSN, AI, 00426 $ LDA, BI, LDA, MM, ALPHA, BETA, Q, LDA, Z, 00427 $ LDA, PL, DIFEST, WORK, LWORK, RWORK, 00428 $ IWORK, LIWORK, BWORK, LINFO ) 00429 * 00430 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 00431 RESULT( 1 ) = ULPINV 00432 WRITE( NOUT, FMT = 9999 )'ZGGESX', LINFO, MPLUSN, 00433 $ PRTYPE 00434 INFO = LINFO 00435 GO TO 30 00436 END IF 00437 * 00438 * Compute the norm(A, B) 00439 * 00440 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, 00441 $ MPLUSN ) 00442 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00443 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00444 ABNRM = ZLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, 00445 $ RWORK ) 00446 * 00447 * Do tests (1) to (4) 00448 * 00449 RESULT( 2 ) = ZERO 00450 CALL ZGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, 00451 $ LDA, WORK, RWORK, RESULT( 1 ) ) 00452 CALL ZGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, 00453 $ LDA, WORK, RWORK, RESULT( 2 ) ) 00454 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, 00455 $ LDA, WORK, RWORK, RESULT( 3 ) ) 00456 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, 00457 $ LDA, WORK, RWORK, RESULT( 4 ) ) 00458 NTEST = 4 00459 * 00460 * Do tests (5) and (6): check Schur form of A and 00461 * compare eigenvalues with diagonals. 00462 * 00463 TEMP1 = ZERO 00464 RESULT( 5 ) = ZERO 00465 RESULT( 6 ) = ZERO 00466 * 00467 DO 10 J = 1, MPLUSN 00468 ILABAD = .FALSE. 00469 TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) / 00470 $ MAX( SMLNUM, ABS1( ALPHA( J ) ), 00471 $ ABS1( AI( J, J ) ) )+ 00472 $ ABS1( BETA( J )-BI( J, J ) ) / 00473 $ MAX( SMLNUM, ABS1( BETA( J ) ), 00474 $ ABS1( BI( J, J ) ) ) ) / ULP 00475 IF( J.LT.MPLUSN ) THEN 00476 IF( AI( J+1, J ).NE.ZERO ) THEN 00477 ILABAD = .TRUE. 00478 RESULT( 5 ) = ULPINV 00479 END IF 00480 END IF 00481 IF( J.GT.1 ) THEN 00482 IF( AI( J, J-1 ).NE.ZERO ) THEN 00483 ILABAD = .TRUE. 00484 RESULT( 5 ) = ULPINV 00485 END IF 00486 END IF 00487 TEMP1 = MAX( TEMP1, TEMP2 ) 00488 IF( ILABAD ) THEN 00489 WRITE( NOUT, FMT = 9997 )J, MPLUSN, PRTYPE 00490 END IF 00491 10 CONTINUE 00492 RESULT( 6 ) = TEMP1 00493 NTEST = NTEST + 2 00494 * 00495 * Test (7) (if sorting worked) 00496 * 00497 RESULT( 7 ) = ZERO 00498 IF( LINFO.EQ.MPLUSN+3 ) THEN 00499 RESULT( 7 ) = ULPINV 00500 ELSE IF( MM.NE.N ) THEN 00501 RESULT( 7 ) = ULPINV 00502 END IF 00503 NTEST = NTEST + 1 00504 * 00505 * Test (8): compare the estimated value DIF and its 00506 * value. first, compute the exact DIF. 00507 * 00508 RESULT( 8 ) = ZERO 00509 MN2 = MM*( MPLUSN-MM )*2 00510 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN 00511 * 00512 * Note: for either following two cases, there are 00513 * almost same number of test cases fail the test. 00514 * 00515 CALL ZLAKF2( MM, MPLUSN-MM, AI, LDA, 00516 $ AI( MM+1, MM+1 ), BI, 00517 $ BI( MM+1, MM+1 ), C, LDC ) 00518 * 00519 CALL ZGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK, 00520 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2, 00521 $ RWORK, INFO ) 00522 DIFTRU = S( MN2 ) 00523 * 00524 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00525 IF( DIFTRU.GT.ABNRM*ULP ) 00526 $ RESULT( 8 ) = ULPINV 00527 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00528 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00529 $ RESULT( 8 ) = ULPINV 00530 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00531 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00532 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), 00533 $ DIFEST( 2 ) / DIFTRU ) 00534 END IF 00535 NTEST = NTEST + 1 00536 END IF 00537 * 00538 * Test (9) 00539 * 00540 RESULT( 9 ) = ZERO 00541 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00542 IF( DIFTRU.GT.ABNRM*ULP ) 00543 $ RESULT( 9 ) = ULPINV 00544 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00545 $ RESULT( 9 ) = ULPINV 00546 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00547 $ RESULT( 9 ) = ULPINV 00548 NTEST = NTEST + 1 00549 END IF 00550 * 00551 NTESTT = NTESTT + NTEST 00552 * 00553 * Print out tests which fail. 00554 * 00555 DO 20 J = 1, 9 00556 IF( RESULT( J ).GE.THRESH ) THEN 00557 * 00558 * If this is the first test to fail, 00559 * print a header to the data file. 00560 * 00561 IF( NERRS.EQ.0 ) THEN 00562 WRITE( NOUT, FMT = 9996 )'ZGX' 00563 * 00564 * Matrix types 00565 * 00566 WRITE( NOUT, FMT = 9994 ) 00567 * 00568 * Tests performed 00569 * 00570 WRITE( NOUT, FMT = 9993 )'unitary', '''', 00571 $ 'transpose', ( '''', I = 1, 4 ) 00572 * 00573 END IF 00574 NERRS = NERRS + 1 00575 IF( RESULT( J ).LT.10000.0D0 ) THEN 00576 WRITE( NOUT, FMT = 9992 )MPLUSN, PRTYPE, 00577 $ WEIGHT, M, J, RESULT( J ) 00578 ELSE 00579 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE, 00580 $ WEIGHT, M, J, RESULT( J ) 00581 END IF 00582 END IF 00583 20 CONTINUE 00584 * 00585 30 CONTINUE 00586 40 CONTINUE 00587 50 CONTINUE 00588 60 CONTINUE 00589 * 00590 GO TO 150 00591 * 00592 70 CONTINUE 00593 * 00594 * Read in data from file to check accuracy of condition estimation 00595 * Read input data until N=0 00596 * 00597 NPTKNT = 0 00598 * 00599 80 CONTINUE 00600 READ( NIN, FMT = *, END = 140 )MPLUSN 00601 IF( MPLUSN.EQ.0 ) 00602 $ GO TO 140 00603 READ( NIN, FMT = *, END = 140 )N 00604 DO 90 I = 1, MPLUSN 00605 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN ) 00606 90 CONTINUE 00607 DO 100 I = 1, MPLUSN 00608 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN ) 00609 100 CONTINUE 00610 READ( NIN, FMT = * )PLTRU, DIFTRU 00611 * 00612 NPTKNT = NPTKNT + 1 00613 FS = .TRUE. 00614 K = 0 00615 M = MPLUSN - N 00616 * 00617 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00618 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00619 * 00620 * Compute the Schur factorization while swaping the 00621 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 00622 * 00623 CALL ZGGESX( 'V', 'V', 'S', ZLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA, 00624 $ MM, ALPHA, BETA, Q, LDA, Z, LDA, PL, DIFEST, WORK, 00625 $ LWORK, RWORK, IWORK, LIWORK, BWORK, LINFO ) 00626 * 00627 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 00628 RESULT( 1 ) = ULPINV 00629 WRITE( NOUT, FMT = 9998 )'ZGGESX', LINFO, MPLUSN, NPTKNT 00630 GO TO 130 00631 END IF 00632 * 00633 * Compute the norm(A, B) 00634 * (should this be norm of (A,B) or (AI,BI)?) 00635 * 00636 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN ) 00637 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00638 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00639 ABNRM = ZLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, RWORK ) 00640 * 00641 * Do tests (1) to (4) 00642 * 00643 CALL ZGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK, 00644 $ RWORK, RESULT( 1 ) ) 00645 CALL ZGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK, 00646 $ RWORK, RESULT( 2 ) ) 00647 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK, 00648 $ RWORK, RESULT( 3 ) ) 00649 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK, 00650 $ RWORK, RESULT( 4 ) ) 00651 * 00652 * Do tests (5) and (6): check Schur form of A and compare 00653 * eigenvalues with diagonals. 00654 * 00655 NTEST = 6 00656 TEMP1 = ZERO 00657 RESULT( 5 ) = ZERO 00658 RESULT( 6 ) = ZERO 00659 * 00660 DO 110 J = 1, MPLUSN 00661 ILABAD = .FALSE. 00662 TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) / 00663 $ MAX( SMLNUM, ABS1( ALPHA( J ) ), ABS1( AI( J, J ) ) )+ 00664 $ ABS1( BETA( J )-BI( J, J ) ) / 00665 $ MAX( SMLNUM, ABS1( BETA( J ) ), ABS1( BI( J, J ) ) ) ) 00666 $ / ULP 00667 IF( J.LT.MPLUSN ) THEN 00668 IF( AI( J+1, J ).NE.ZERO ) THEN 00669 ILABAD = .TRUE. 00670 RESULT( 5 ) = ULPINV 00671 END IF 00672 END IF 00673 IF( J.GT.1 ) THEN 00674 IF( AI( J, J-1 ).NE.ZERO ) THEN 00675 ILABAD = .TRUE. 00676 RESULT( 5 ) = ULPINV 00677 END IF 00678 END IF 00679 TEMP1 = MAX( TEMP1, TEMP2 ) 00680 IF( ILABAD ) THEN 00681 WRITE( NOUT, FMT = 9997 )J, MPLUSN, NPTKNT 00682 END IF 00683 110 CONTINUE 00684 RESULT( 6 ) = TEMP1 00685 * 00686 * Test (7) (if sorting worked) <--------- need to be checked. 00687 * 00688 NTEST = 7 00689 RESULT( 7 ) = ZERO 00690 IF( LINFO.EQ.MPLUSN+3 ) 00691 $ RESULT( 7 ) = ULPINV 00692 * 00693 * Test (8): compare the estimated value of DIF and its true value. 00694 * 00695 NTEST = 8 00696 RESULT( 8 ) = ZERO 00697 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00698 IF( DIFTRU.GT.ABNRM*ULP ) 00699 $ RESULT( 8 ) = ULPINV 00700 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00701 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00702 $ RESULT( 8 ) = ULPINV 00703 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00704 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00705 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU ) 00706 END IF 00707 * 00708 * Test (9) 00709 * 00710 NTEST = 9 00711 RESULT( 9 ) = ZERO 00712 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00713 IF( DIFTRU.GT.ABNRM*ULP ) 00714 $ RESULT( 9 ) = ULPINV 00715 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00716 $ RESULT( 9 ) = ULPINV 00717 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00718 $ RESULT( 9 ) = ULPINV 00719 END IF 00720 * 00721 * Test (10): compare the estimated value of PL and it true value. 00722 * 00723 NTEST = 10 00724 RESULT( 10 ) = ZERO 00725 IF( PL( 1 ).EQ.ZERO ) THEN 00726 IF( PLTRU.GT.ABNRM*ULP ) 00727 $ RESULT( 10 ) = ULPINV 00728 ELSE IF( PLTRU.EQ.ZERO ) THEN 00729 IF( PL( 1 ).GT.ABNRM*ULP ) 00730 $ RESULT( 10 ) = ULPINV 00731 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR. 00732 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN 00733 RESULT( 10 ) = ULPINV 00734 END IF 00735 * 00736 NTESTT = NTESTT + NTEST 00737 * 00738 * Print out tests which fail. 00739 * 00740 DO 120 J = 1, NTEST 00741 IF( RESULT( J ).GE.THRESH ) THEN 00742 * 00743 * If this is the first test to fail, 00744 * print a header to the data file. 00745 * 00746 IF( NERRS.EQ.0 ) THEN 00747 WRITE( NOUT, FMT = 9996 )'ZGX' 00748 * 00749 * Matrix types 00750 * 00751 WRITE( NOUT, FMT = 9995 ) 00752 * 00753 * Tests performed 00754 * 00755 WRITE( NOUT, FMT = 9993 )'unitary', '''', 'transpose', 00756 $ ( '''', I = 1, 4 ) 00757 * 00758 END IF 00759 NERRS = NERRS + 1 00760 IF( RESULT( J ).LT.10000.0D0 ) THEN 00761 WRITE( NOUT, FMT = 9990 )NPTKNT, MPLUSN, J, RESULT( J ) 00762 ELSE 00763 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J ) 00764 END IF 00765 END IF 00766 * 00767 120 CONTINUE 00768 * 00769 130 CONTINUE 00770 GO TO 80 00771 140 CONTINUE 00772 * 00773 150 CONTINUE 00774 * 00775 * Summary 00776 * 00777 CALL ALASVM( 'ZGX', NOUT, NERRS, NTESTT, 0 ) 00778 * 00779 WORK( 1 ) = MAXWRK 00780 * 00781 RETURN 00782 * 00783 9999 FORMAT( ' ZDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00784 $ I6, ', JTYPE=', I6, ')' ) 00785 * 00786 9998 FORMAT( ' ZDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00787 $ I6, ', Input Example #', I2, ')' ) 00788 * 00789 9997 FORMAT( ' ZDRGSX: S not in Schur form at eigenvalue ', I6, '.', 00790 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' ) 00791 * 00792 9996 FORMAT( / 1X, A3, ' -- Complex Expert Generalized Schur form', 00793 $ ' problem driver' ) 00794 * 00795 9995 FORMAT( 'Input Example' ) 00796 * 00797 9994 FORMAT( ' Matrix types: ', / 00798 $ ' 1: A is a block diagonal matrix of Jordan blocks ', 00799 $ 'and B is the identity ', / ' matrix, ', 00800 $ / ' 2: A and B are upper triangular matrices, ', 00801 $ / ' 3: A and B are as type 2, but each second diagonal ', 00802 $ 'block in A_11 and ', / 00803 $ ' each third diaongal block in A_22 are 2x2 blocks,', 00804 $ / ' 4: A and B are block diagonal matrices, ', 00805 $ / ' 5: (A,B) has potentially close or common ', 00806 $ 'eigenvalues.', / ) 00807 * 00808 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ', 00809 $ 'Q and Z are ', A, ',', / 19X, 00810 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)', 00811 $ / ' 1 = | A - Q S Z', A, 00812 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A, 00813 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A, 00814 $ ' | / ( n ulp ) 4 = | I - ZZ', A, 00815 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ', 00816 $ 'Schur form S', / ' 6 = difference between (alpha,beta)', 00817 $ ' and diagonals of (S,T)', / 00818 $ ' 7 = 1/ULP if SDIM is not the correct number of ', 00819 $ 'selected eigenvalues', / 00820 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ', 00821 $ 'DIFTRU/DIFEST > 10*THRESH', 00822 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ', 00823 $ 'when reordering fails', / 00824 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ', 00825 $ 'PLTRU/PLEST > THRESH', / 00826 $ ' ( Test 10 is only for input examples )', / ) 00827 9992 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3, 00828 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 ) 00829 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3, 00830 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.3 ) 00831 9990 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00832 $ ' result ', I2, ' is', 0P, F8.2 ) 00833 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00834 $ ' result ', I2, ' is', 1P, D10.3 ) 00835 * 00836 * End of ZDRGSX 00837 * 00838 END
1.7.3 
