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