|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, 00002 $ H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN, 00003 $ RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK, 00004 $ LWORK, RWORK, BWORK, INFO ) 00005 * 00006 * -- LAPACK test routine (version 3.1) -- 00007 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 LOGICAL COMP 00012 INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, 00013 $ NSLCT 00014 REAL RCDEIN, RCDVIN, THRESH 00015 * .. 00016 * .. Array Arguments .. 00017 LOGICAL BWORK( * ) 00018 INTEGER ISEED( 4 ), ISLCT( * ) 00019 REAL RESULT( 17 ), RWORK( * ) 00020 COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ), 00021 $ VS( LDVS, * ), VS1( LDVS, * ), W( * ), 00022 $ WORK( * ), WT( * ), WTMP( * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * CGET24 checks the nonsymmetric eigenvalue (Schur form) problem 00029 * expert driver CGEESX. 00030 * 00031 * If COMP = .FALSE., the first 13 of the following tests will be 00032 * be performed on the input matrix A, and also tests 14 and 15 00033 * if LWORK is sufficiently large. 00034 * If COMP = .TRUE., all 17 test will be performed. 00035 * 00036 * (1) 0 if T is in Schur form, 1/ulp otherwise 00037 * (no sorting of eigenvalues) 00038 * 00039 * (2) | A - VS T VS' | / ( n |A| ulp ) 00040 * 00041 * Here VS is the matrix of Schur eigenvectors, and T is in Schur 00042 * form (no sorting of eigenvalues). 00043 * 00044 * (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues). 00045 * 00046 * (4) 0 if W are eigenvalues of T 00047 * 1/ulp otherwise 00048 * (no sorting of eigenvalues) 00049 * 00050 * (5) 0 if T(with VS) = T(without VS), 00051 * 1/ulp otherwise 00052 * (no sorting of eigenvalues) 00053 * 00054 * (6) 0 if eigenvalues(with VS) = eigenvalues(without VS), 00055 * 1/ulp otherwise 00056 * (no sorting of eigenvalues) 00057 * 00058 * (7) 0 if T is in Schur form, 1/ulp otherwise 00059 * (with sorting of eigenvalues) 00060 * 00061 * (8) | A - VS T VS' | / ( n |A| ulp ) 00062 * 00063 * Here VS is the matrix of Schur eigenvectors, and T is in Schur 00064 * form (with sorting of eigenvalues). 00065 * 00066 * (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues). 00067 * 00068 * (10) 0 if W are eigenvalues of T 00069 * 1/ulp otherwise 00070 * If workspace sufficient, also compare W with and 00071 * without reciprocal condition numbers 00072 * (with sorting of eigenvalues) 00073 * 00074 * (11) 0 if T(with VS) = T(without VS), 00075 * 1/ulp otherwise 00076 * If workspace sufficient, also compare T with and without 00077 * reciprocal condition numbers 00078 * (with sorting of eigenvalues) 00079 * 00080 * (12) 0 if eigenvalues(with VS) = eigenvalues(without VS), 00081 * 1/ulp otherwise 00082 * If workspace sufficient, also compare VS with and without 00083 * reciprocal condition numbers 00084 * (with sorting of eigenvalues) 00085 * 00086 * (13) if sorting worked and SDIM is the number of 00087 * eigenvalues which were SELECTed 00088 * If workspace sufficient, also compare SDIM with and 00089 * without reciprocal condition numbers 00090 * 00091 * (14) if RCONDE the same no matter if VS and/or RCONDV computed 00092 * 00093 * (15) if RCONDV the same no matter if VS and/or RCONDE computed 00094 * 00095 * (16) |RCONDE - RCDEIN| / cond(RCONDE) 00096 * 00097 * RCONDE is the reciprocal average eigenvalue condition number 00098 * computed by CGEESX and RCDEIN (the precomputed true value) 00099 * is supplied as input. cond(RCONDE) is the condition number 00100 * of RCONDE, and takes errors in computing RCONDE into account, 00101 * so that the resulting quantity should be O(ULP). cond(RCONDE) 00102 * is essentially given by norm(A)/RCONDV. 00103 * 00104 * (17) |RCONDV - RCDVIN| / cond(RCONDV) 00105 * 00106 * RCONDV is the reciprocal right invariant subspace condition 00107 * number computed by CGEESX and RCDVIN (the precomputed true 00108 * value) is supplied as input. cond(RCONDV) is the condition 00109 * number of RCONDV, and takes errors in computing RCONDV into 00110 * account, so that the resulting quantity should be O(ULP). 00111 * cond(RCONDV) is essentially given by norm(A)/RCONDE. 00112 * 00113 * Arguments 00114 * ========= 00115 * 00116 * COMP (input) LOGICAL 00117 * COMP describes which input tests to perform: 00118 * = .FALSE. if the computed condition numbers are not to 00119 * be tested against RCDVIN and RCDEIN 00120 * = .TRUE. if they are to be compared 00121 * 00122 * JTYPE (input) INTEGER 00123 * Type of input matrix. Used to label output if error occurs. 00124 * 00125 * ISEED (input) INTEGER array, dimension (4) 00126 * If COMP = .FALSE., the random number generator seed 00127 * used to produce matrix. 00128 * If COMP = .TRUE., ISEED(1) = the number of the example. 00129 * Used to label output if error occurs. 00130 * 00131 * THRESH (input) REAL 00132 * A test will count as "failed" if the "error", computed as 00133 * described above, exceeds THRESH. Note that the error 00134 * is scaled to be O(1), so THRESH should be a reasonably 00135 * small multiple of 1, e.g., 10 or 100. In particular, 00136 * it should not depend on the precision (single vs. double) 00137 * or the size of the matrix. It must be at least zero. 00138 * 00139 * NOUNIT (input) INTEGER 00140 * The FORTRAN unit number for printing out error messages 00141 * (e.g., if a routine returns INFO not equal to 0.) 00142 * 00143 * N (input) INTEGER 00144 * The dimension of A. N must be at least 0. 00145 * 00146 * A (input/output) COMPLEX array, dimension (LDA, N) 00147 * Used to hold the matrix whose eigenvalues are to be 00148 * computed. 00149 * 00150 * LDA (input) INTEGER 00151 * The leading dimension of A, and H. LDA must be at 00152 * least 1 and at least N. 00153 * 00154 * H (workspace) COMPLEX array, dimension (LDA, N) 00155 * Another copy of the test matrix A, modified by CGEESX. 00156 * 00157 * HT (workspace) COMPLEX array, dimension (LDA, N) 00158 * Yet another copy of the test matrix A, modified by CGEESX. 00159 * 00160 * W (workspace) COMPLEX array, dimension (N) 00161 * The computed eigenvalues of A. 00162 * 00163 * WT (workspace) COMPLEX array, dimension (N) 00164 * Like W, this array contains the eigenvalues of A, 00165 * but those computed when CGEESX only computes a partial 00166 * eigendecomposition, i.e. not Schur vectors 00167 * 00168 * WTMP (workspace) COMPLEX array, dimension (N) 00169 * Like W, this array contains the eigenvalues of A, 00170 * but sorted by increasing real or imaginary part. 00171 * 00172 * VS (workspace) COMPLEX array, dimension (LDVS, N) 00173 * VS holds the computed Schur vectors. 00174 * 00175 * LDVS (input) INTEGER 00176 * Leading dimension of VS. Must be at least max(1, N). 00177 * 00178 * VS1 (workspace) COMPLEX array, dimension (LDVS, N) 00179 * VS1 holds another copy of the computed Schur vectors. 00180 * 00181 * RCDEIN (input) REAL 00182 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal 00183 * condition number for the average of selected eigenvalues. 00184 * 00185 * RCDVIN (input) REAL 00186 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal 00187 * condition number for the selected right invariant subspace. 00188 * 00189 * NSLCT (input) INTEGER 00190 * When COMP = .TRUE. the number of selected eigenvalues 00191 * corresponding to the precomputed values RCDEIN and RCDVIN. 00192 * 00193 * ISLCT (input) INTEGER array, dimension (NSLCT) 00194 * When COMP = .TRUE. ISLCT selects the eigenvalues of the 00195 * input matrix corresponding to the precomputed values RCDEIN 00196 * and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the 00197 * eigenvalue with the J-th largest real or imaginary part is 00198 * selected. The real part is used if ISRT = 0, and the 00199 * imaginary part if ISRT = 1. 00200 * Not referenced if COMP = .FALSE. 00201 * 00202 * ISRT (input) INTEGER 00203 * When COMP = .TRUE., ISRT describes how ISLCT is used to 00204 * choose a subset of the spectrum. 00205 * Not referenced if COMP = .FALSE. 00206 * 00207 * RESULT (output) REAL array, dimension (17) 00208 * The values computed by the 17 tests described above. 00209 * The values are currently limited to 1/ulp, to avoid 00210 * overflow. 00211 * 00212 * WORK (workspace) COMPLEX array, dimension (2*N*N) 00213 * 00214 * LWORK (input) INTEGER 00215 * The number of entries in WORK to be passed to CGEESX. This 00216 * must be at least 2*N, and N*(N+1)/2 if tests 14--16 are to 00217 * be performed. 00218 * 00219 * RWORK (workspace) REAL array, dimension (N) 00220 * 00221 * BWORK (workspace) LOGICAL array, dimension (N) 00222 * 00223 * INFO (output) INTEGER 00224 * If 0, successful exit. 00225 * If <0, input parameter -INFO had an incorrect value. 00226 * If >0, CGEESX returned an error code, the absolute 00227 * value of which is returned. 00228 * 00229 * ===================================================================== 00230 * 00231 * .. Parameters .. 00232 COMPLEX CZERO, CONE 00233 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00234 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00235 REAL ZERO, ONE 00236 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00237 REAL EPSIN 00238 PARAMETER ( EPSIN = 5.9605E-8 ) 00239 * .. 00240 * .. Local Scalars .. 00241 CHARACTER SORT 00242 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB, 00243 $ SDIM, SDIM1 00244 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV, 00245 $ SMLNUM, TOL, TOLIN, ULP, ULPINV, V, VRICMP, 00246 $ VRIMIN, WNORM 00247 COMPLEX CTMP 00248 * .. 00249 * .. Local Arrays .. 00250 INTEGER IPNT( 20 ) 00251 * .. 00252 * .. External Functions .. 00253 LOGICAL CSLECT 00254 REAL CLANGE, SLAMCH 00255 EXTERNAL CSLECT, CLANGE, SLAMCH 00256 * .. 00257 * .. External Subroutines .. 00258 EXTERNAL CCOPY, CGEESX, CGEMM, CLACPY, CUNT01, XERBLA 00259 * .. 00260 * .. Intrinsic Functions .. 00261 INTRINSIC ABS, AIMAG, MAX, MIN, REAL 00262 * .. 00263 * .. Arrays in Common .. 00264 LOGICAL SELVAL( 20 ) 00265 REAL SELWI( 20 ), SELWR( 20 ) 00266 * .. 00267 * .. Scalars in Common .. 00268 INTEGER SELDIM, SELOPT 00269 * .. 00270 * .. Common blocks .. 00271 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI 00272 * .. 00273 * .. Executable Statements .. 00274 * 00275 * Check for errors 00276 * 00277 INFO = 0 00278 IF( THRESH.LT.ZERO ) THEN 00279 INFO = -3 00280 ELSE IF( NOUNIT.LE.0 ) THEN 00281 INFO = -5 00282 ELSE IF( N.LT.0 ) THEN 00283 INFO = -6 00284 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN 00285 INFO = -8 00286 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN 00287 INFO = -15 00288 ELSE IF( LWORK.LT.2*N ) THEN 00289 INFO = -24 00290 END IF 00291 * 00292 IF( INFO.NE.0 ) THEN 00293 CALL XERBLA( 'CGET24', -INFO ) 00294 RETURN 00295 END IF 00296 * 00297 * Quick return if nothing to do 00298 * 00299 DO 10 I = 1, 17 00300 RESULT( I ) = -ONE 00301 10 CONTINUE 00302 * 00303 IF( N.EQ.0 ) 00304 $ RETURN 00305 * 00306 * Important constants 00307 * 00308 SMLNUM = SLAMCH( 'Safe minimum' ) 00309 ULP = SLAMCH( 'Precision' ) 00310 ULPINV = ONE / ULP 00311 * 00312 * Perform tests (1)-(13) 00313 * 00314 SELOPT = 0 00315 DO 90 ISORT = 0, 1 00316 IF( ISORT.EQ.0 ) THEN 00317 SORT = 'N' 00318 RSUB = 0 00319 ELSE 00320 SORT = 'S' 00321 RSUB = 6 00322 END IF 00323 * 00324 * Compute Schur form and Schur vectors, and test them 00325 * 00326 CALL CLACPY( 'F', N, N, A, LDA, H, LDA ) 00327 CALL CGEESX( 'V', SORT, CSLECT, 'N', N, H, LDA, SDIM, W, VS, 00328 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK, 00329 $ IINFO ) 00330 IF( IINFO.NE.0 ) THEN 00331 RESULT( 1+RSUB ) = ULPINV 00332 IF( JTYPE.NE.22 ) THEN 00333 WRITE( NOUNIT, FMT = 9998 )'CGEESX1', IINFO, N, JTYPE, 00334 $ ISEED 00335 ELSE 00336 WRITE( NOUNIT, FMT = 9999 )'CGEESX1', IINFO, N, 00337 $ ISEED( 1 ) 00338 END IF 00339 INFO = ABS( IINFO ) 00340 RETURN 00341 END IF 00342 IF( ISORT.EQ.0 ) THEN 00343 CALL CCOPY( N, W, 1, WTMP, 1 ) 00344 END IF 00345 * 00346 * Do Test (1) or Test (7) 00347 * 00348 RESULT( 1+RSUB ) = ZERO 00349 DO 30 J = 1, N - 1 00350 DO 20 I = J + 1, N 00351 IF( H( I, J ).NE.CZERO ) 00352 $ RESULT( 1+RSUB ) = ULPINV 00353 20 CONTINUE 00354 30 CONTINUE 00355 * 00356 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP) 00357 * 00358 * Copy A to VS1, used as workspace 00359 * 00360 CALL CLACPY( ' ', N, N, A, LDA, VS1, LDVS ) 00361 * 00362 * Compute Q*H and store in HT. 00363 * 00364 CALL CGEMM( 'No transpose', 'No transpose', N, N, N, CONE, VS, 00365 $ LDVS, H, LDA, CZERO, HT, LDA ) 00366 * 00367 * Compute A - Q*H*Q' 00368 * 00369 CALL CGEMM( 'No transpose', 'Conjugate transpose', N, N, N, 00370 $ -CONE, HT, LDA, VS, LDVS, CONE, VS1, LDVS ) 00371 * 00372 ANORM = MAX( CLANGE( '1', N, N, A, LDA, RWORK ), SMLNUM ) 00373 WNORM = CLANGE( '1', N, N, VS1, LDVS, RWORK ) 00374 * 00375 IF( ANORM.GT.WNORM ) THEN 00376 RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP ) 00377 ELSE 00378 IF( ANORM.LT.ONE ) THEN 00379 RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / 00380 $ ( N*ULP ) 00381 ELSE 00382 RESULT( 2+RSUB ) = MIN( WNORM / ANORM, REAL( N ) ) / 00383 $ ( N*ULP ) 00384 END IF 00385 END IF 00386 * 00387 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP ) 00388 * 00389 CALL CUNT01( 'Columns', N, N, VS, LDVS, WORK, LWORK, RWORK, 00390 $ RESULT( 3+RSUB ) ) 00391 * 00392 * Do Test (4) or Test (10) 00393 * 00394 RESULT( 4+RSUB ) = ZERO 00395 DO 40 I = 1, N 00396 IF( H( I, I ).NE.W( I ) ) 00397 $ RESULT( 4+RSUB ) = ULPINV 00398 40 CONTINUE 00399 * 00400 * Do Test (5) or Test (11) 00401 * 00402 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00403 CALL CGEESX( 'N', SORT, CSLECT, 'N', N, HT, LDA, SDIM, WT, VS, 00404 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK, 00405 $ IINFO ) 00406 IF( IINFO.NE.0 ) THEN 00407 RESULT( 5+RSUB ) = ULPINV 00408 IF( JTYPE.NE.22 ) THEN 00409 WRITE( NOUNIT, FMT = 9998 )'CGEESX2', IINFO, N, JTYPE, 00410 $ ISEED 00411 ELSE 00412 WRITE( NOUNIT, FMT = 9999 )'CGEESX2', IINFO, N, 00413 $ ISEED( 1 ) 00414 END IF 00415 INFO = ABS( IINFO ) 00416 GO TO 220 00417 END IF 00418 * 00419 RESULT( 5+RSUB ) = ZERO 00420 DO 60 J = 1, N 00421 DO 50 I = 1, N 00422 IF( H( I, J ).NE.HT( I, J ) ) 00423 $ RESULT( 5+RSUB ) = ULPINV 00424 50 CONTINUE 00425 60 CONTINUE 00426 * 00427 * Do Test (6) or Test (12) 00428 * 00429 RESULT( 6+RSUB ) = ZERO 00430 DO 70 I = 1, N 00431 IF( W( I ).NE.WT( I ) ) 00432 $ RESULT( 6+RSUB ) = ULPINV 00433 70 CONTINUE 00434 * 00435 * Do Test (13) 00436 * 00437 IF( ISORT.EQ.1 ) THEN 00438 RESULT( 13 ) = ZERO 00439 KNTEIG = 0 00440 DO 80 I = 1, N 00441 IF( CSLECT( W( I ) ) ) 00442 $ KNTEIG = KNTEIG + 1 00443 IF( I.LT.N ) THEN 00444 IF( CSLECT( W( I+1 ) ) .AND. 00445 $ ( .NOT.CSLECT( W( I ) ) ) )RESULT( 13 ) = ULPINV 00446 END IF 00447 80 CONTINUE 00448 IF( SDIM.NE.KNTEIG ) 00449 $ RESULT( 13 ) = ULPINV 00450 END IF 00451 * 00452 90 CONTINUE 00453 * 00454 * If there is enough workspace, perform tests (14) and (15) 00455 * as well as (10) through (13) 00456 * 00457 IF( LWORK.GE.( N*( N+1 ) ) / 2 ) THEN 00458 * 00459 * Compute both RCONDE and RCONDV with VS 00460 * 00461 SORT = 'S' 00462 RESULT( 14 ) = ZERO 00463 RESULT( 15 ) = ZERO 00464 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00465 CALL CGEESX( 'V', SORT, CSLECT, 'B', N, HT, LDA, SDIM1, WT, 00466 $ VS1, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, 00467 $ BWORK, IINFO ) 00468 IF( IINFO.NE.0 ) THEN 00469 RESULT( 14 ) = ULPINV 00470 RESULT( 15 ) = ULPINV 00471 IF( JTYPE.NE.22 ) THEN 00472 WRITE( NOUNIT, FMT = 9998 )'CGEESX3', IINFO, N, JTYPE, 00473 $ ISEED 00474 ELSE 00475 WRITE( NOUNIT, FMT = 9999 )'CGEESX3', IINFO, N, 00476 $ ISEED( 1 ) 00477 END IF 00478 INFO = ABS( IINFO ) 00479 GO TO 220 00480 END IF 00481 * 00482 * Perform tests (10), (11), (12), and (13) 00483 * 00484 DO 110 I = 1, N 00485 IF( W( I ).NE.WT( I ) ) 00486 $ RESULT( 10 ) = ULPINV 00487 DO 100 J = 1, N 00488 IF( H( I, J ).NE.HT( I, J ) ) 00489 $ RESULT( 11 ) = ULPINV 00490 IF( VS( I, J ).NE.VS1( I, J ) ) 00491 $ RESULT( 12 ) = ULPINV 00492 100 CONTINUE 00493 110 CONTINUE 00494 IF( SDIM.NE.SDIM1 ) 00495 $ RESULT( 13 ) = ULPINV 00496 * 00497 * Compute both RCONDE and RCONDV without VS, and compare 00498 * 00499 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00500 CALL CGEESX( 'N', SORT, CSLECT, 'B', N, HT, LDA, SDIM1, WT, 00501 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK, 00502 $ BWORK, IINFO ) 00503 IF( IINFO.NE.0 ) THEN 00504 RESULT( 14 ) = ULPINV 00505 RESULT( 15 ) = ULPINV 00506 IF( JTYPE.NE.22 ) THEN 00507 WRITE( NOUNIT, FMT = 9998 )'CGEESX4', IINFO, N, JTYPE, 00508 $ ISEED 00509 ELSE 00510 WRITE( NOUNIT, FMT = 9999 )'CGEESX4', IINFO, N, 00511 $ ISEED( 1 ) 00512 END IF 00513 INFO = ABS( IINFO ) 00514 GO TO 220 00515 END IF 00516 * 00517 * Perform tests (14) and (15) 00518 * 00519 IF( RCNDE1.NE.RCONDE ) 00520 $ RESULT( 14 ) = ULPINV 00521 IF( RCNDV1.NE.RCONDV ) 00522 $ RESULT( 15 ) = ULPINV 00523 * 00524 * Perform tests (10), (11), (12), and (13) 00525 * 00526 DO 130 I = 1, N 00527 IF( W( I ).NE.WT( I ) ) 00528 $ RESULT( 10 ) = ULPINV 00529 DO 120 J = 1, N 00530 IF( H( I, J ).NE.HT( I, J ) ) 00531 $ RESULT( 11 ) = ULPINV 00532 IF( VS( I, J ).NE.VS1( I, J ) ) 00533 $ RESULT( 12 ) = ULPINV 00534 120 CONTINUE 00535 130 CONTINUE 00536 IF( SDIM.NE.SDIM1 ) 00537 $ RESULT( 13 ) = ULPINV 00538 * 00539 * Compute RCONDE with VS, and compare 00540 * 00541 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00542 CALL CGEESX( 'V', SORT, CSLECT, 'E', N, HT, LDA, SDIM1, WT, 00543 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK, 00544 $ BWORK, IINFO ) 00545 IF( IINFO.NE.0 ) THEN 00546 RESULT( 14 ) = ULPINV 00547 IF( JTYPE.NE.22 ) THEN 00548 WRITE( NOUNIT, FMT = 9998 )'CGEESX5', IINFO, N, JTYPE, 00549 $ ISEED 00550 ELSE 00551 WRITE( NOUNIT, FMT = 9999 )'CGEESX5', IINFO, N, 00552 $ ISEED( 1 ) 00553 END IF 00554 INFO = ABS( IINFO ) 00555 GO TO 220 00556 END IF 00557 * 00558 * Perform test (14) 00559 * 00560 IF( RCNDE1.NE.RCONDE ) 00561 $ RESULT( 14 ) = ULPINV 00562 * 00563 * Perform tests (10), (11), (12), and (13) 00564 * 00565 DO 150 I = 1, N 00566 IF( W( I ).NE.WT( I ) ) 00567 $ RESULT( 10 ) = ULPINV 00568 DO 140 J = 1, N 00569 IF( H( I, J ).NE.HT( I, J ) ) 00570 $ RESULT( 11 ) = ULPINV 00571 IF( VS( I, J ).NE.VS1( I, J ) ) 00572 $ RESULT( 12 ) = ULPINV 00573 140 CONTINUE 00574 150 CONTINUE 00575 IF( SDIM.NE.SDIM1 ) 00576 $ RESULT( 13 ) = ULPINV 00577 * 00578 * Compute RCONDE without VS, and compare 00579 * 00580 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00581 CALL CGEESX( 'N', SORT, CSLECT, 'E', N, HT, LDA, SDIM1, WT, 00582 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK, 00583 $ BWORK, IINFO ) 00584 IF( IINFO.NE.0 ) THEN 00585 RESULT( 14 ) = ULPINV 00586 IF( JTYPE.NE.22 ) THEN 00587 WRITE( NOUNIT, FMT = 9998 )'CGEESX6', IINFO, N, JTYPE, 00588 $ ISEED 00589 ELSE 00590 WRITE( NOUNIT, FMT = 9999 )'CGEESX6', IINFO, N, 00591 $ ISEED( 1 ) 00592 END IF 00593 INFO = ABS( IINFO ) 00594 GO TO 220 00595 END IF 00596 * 00597 * Perform test (14) 00598 * 00599 IF( RCNDE1.NE.RCONDE ) 00600 $ RESULT( 14 ) = ULPINV 00601 * 00602 * Perform tests (10), (11), (12), and (13) 00603 * 00604 DO 170 I = 1, N 00605 IF( W( I ).NE.WT( I ) ) 00606 $ RESULT( 10 ) = ULPINV 00607 DO 160 J = 1, N 00608 IF( H( I, J ).NE.HT( I, J ) ) 00609 $ RESULT( 11 ) = ULPINV 00610 IF( VS( I, J ).NE.VS1( I, J ) ) 00611 $ RESULT( 12 ) = ULPINV 00612 160 CONTINUE 00613 170 CONTINUE 00614 IF( SDIM.NE.SDIM1 ) 00615 $ RESULT( 13 ) = ULPINV 00616 * 00617 * Compute RCONDV with VS, and compare 00618 * 00619 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00620 CALL CGEESX( 'V', SORT, CSLECT, 'V', N, HT, LDA, SDIM1, WT, 00621 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK, 00622 $ BWORK, IINFO ) 00623 IF( IINFO.NE.0 ) THEN 00624 RESULT( 15 ) = ULPINV 00625 IF( JTYPE.NE.22 ) THEN 00626 WRITE( NOUNIT, FMT = 9998 )'CGEESX7', IINFO, N, JTYPE, 00627 $ ISEED 00628 ELSE 00629 WRITE( NOUNIT, FMT = 9999 )'CGEESX7', IINFO, N, 00630 $ ISEED( 1 ) 00631 END IF 00632 INFO = ABS( IINFO ) 00633 GO TO 220 00634 END IF 00635 * 00636 * Perform test (15) 00637 * 00638 IF( RCNDV1.NE.RCONDV ) 00639 $ RESULT( 15 ) = ULPINV 00640 * 00641 * Perform tests (10), (11), (12), and (13) 00642 * 00643 DO 190 I = 1, N 00644 IF( W( I ).NE.WT( I ) ) 00645 $ RESULT( 10 ) = ULPINV 00646 DO 180 J = 1, N 00647 IF( H( I, J ).NE.HT( I, J ) ) 00648 $ RESULT( 11 ) = ULPINV 00649 IF( VS( I, J ).NE.VS1( I, J ) ) 00650 $ RESULT( 12 ) = ULPINV 00651 180 CONTINUE 00652 190 CONTINUE 00653 IF( SDIM.NE.SDIM1 ) 00654 $ RESULT( 13 ) = ULPINV 00655 * 00656 * Compute RCONDV without VS, and compare 00657 * 00658 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00659 CALL CGEESX( 'N', SORT, CSLECT, 'V', N, HT, LDA, SDIM1, WT, 00660 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK, 00661 $ BWORK, IINFO ) 00662 IF( IINFO.NE.0 ) THEN 00663 RESULT( 15 ) = ULPINV 00664 IF( JTYPE.NE.22 ) THEN 00665 WRITE( NOUNIT, FMT = 9998 )'CGEESX8', IINFO, N, JTYPE, 00666 $ ISEED 00667 ELSE 00668 WRITE( NOUNIT, FMT = 9999 )'CGEESX8', IINFO, N, 00669 $ ISEED( 1 ) 00670 END IF 00671 INFO = ABS( IINFO ) 00672 GO TO 220 00673 END IF 00674 * 00675 * Perform test (15) 00676 * 00677 IF( RCNDV1.NE.RCONDV ) 00678 $ RESULT( 15 ) = ULPINV 00679 * 00680 * Perform tests (10), (11), (12), and (13) 00681 * 00682 DO 210 I = 1, N 00683 IF( W( I ).NE.WT( I ) ) 00684 $ RESULT( 10 ) = ULPINV 00685 DO 200 J = 1, N 00686 IF( H( I, J ).NE.HT( I, J ) ) 00687 $ RESULT( 11 ) = ULPINV 00688 IF( VS( I, J ).NE.VS1( I, J ) ) 00689 $ RESULT( 12 ) = ULPINV 00690 200 CONTINUE 00691 210 CONTINUE 00692 IF( SDIM.NE.SDIM1 ) 00693 $ RESULT( 13 ) = ULPINV 00694 * 00695 END IF 00696 * 00697 220 CONTINUE 00698 * 00699 * If there are precomputed reciprocal condition numbers, compare 00700 * computed values with them. 00701 * 00702 IF( COMP ) THEN 00703 * 00704 * First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that 00705 * the logical function CSLECT selects the eigenvalues specified 00706 * by NSLCT, ISLCT and ISRT. 00707 * 00708 SELDIM = N 00709 SELOPT = 1 00710 EPS = MAX( ULP, EPSIN ) 00711 DO 230 I = 1, N 00712 IPNT( I ) = I 00713 SELVAL( I ) = .FALSE. 00714 SELWR( I ) = REAL( WTMP( I ) ) 00715 SELWI( I ) = AIMAG( WTMP( I ) ) 00716 230 CONTINUE 00717 DO 250 I = 1, N - 1 00718 KMIN = I 00719 IF( ISRT.EQ.0 ) THEN 00720 VRIMIN = REAL( WTMP( I ) ) 00721 ELSE 00722 VRIMIN = AIMAG( WTMP( I ) ) 00723 END IF 00724 DO 240 J = I + 1, N 00725 IF( ISRT.EQ.0 ) THEN 00726 VRICMP = REAL( WTMP( J ) ) 00727 ELSE 00728 VRICMP = AIMAG( WTMP( J ) ) 00729 END IF 00730 IF( VRICMP.LT.VRIMIN ) THEN 00731 KMIN = J 00732 VRIMIN = VRICMP 00733 END IF 00734 240 CONTINUE 00735 CTMP = WTMP( KMIN ) 00736 WTMP( KMIN ) = WTMP( I ) 00737 WTMP( I ) = CTMP 00738 ITMP = IPNT( I ) 00739 IPNT( I ) = IPNT( KMIN ) 00740 IPNT( KMIN ) = ITMP 00741 250 CONTINUE 00742 DO 260 I = 1, NSLCT 00743 SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE. 00744 260 CONTINUE 00745 * 00746 * Compute condition numbers 00747 * 00748 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00749 CALL CGEESX( 'N', 'S', CSLECT, 'B', N, HT, LDA, SDIM1, WT, VS1, 00750 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK, 00751 $ IINFO ) 00752 IF( IINFO.NE.0 ) THEN 00753 RESULT( 16 ) = ULPINV 00754 RESULT( 17 ) = ULPINV 00755 WRITE( NOUNIT, FMT = 9999 )'CGEESX9', IINFO, N, ISEED( 1 ) 00756 INFO = ABS( IINFO ) 00757 GO TO 270 00758 END IF 00759 * 00760 * Compare condition number for average of selected eigenvalues 00761 * taking its condition number into account 00762 * 00763 ANORM = CLANGE( '1', N, N, A, LDA, RWORK ) 00764 V = MAX( REAL( N )*EPS*ANORM, SMLNUM ) 00765 IF( ANORM.EQ.ZERO ) 00766 $ V = ONE 00767 IF( V.GT.RCONDV ) THEN 00768 TOL = ONE 00769 ELSE 00770 TOL = V / RCONDV 00771 END IF 00772 IF( V.GT.RCDVIN ) THEN 00773 TOLIN = ONE 00774 ELSE 00775 TOLIN = V / RCDVIN 00776 END IF 00777 TOL = MAX( TOL, SMLNUM / EPS ) 00778 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00779 IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN 00780 RESULT( 16 ) = ULPINV 00781 ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN 00782 RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL ) 00783 ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN 00784 RESULT( 16 ) = ULPINV 00785 ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN 00786 RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN ) 00787 ELSE 00788 RESULT( 16 ) = ONE 00789 END IF 00790 * 00791 * Compare condition numbers for right invariant subspace 00792 * taking its condition number into account 00793 * 00794 IF( V.GT.RCONDV*RCONDE ) THEN 00795 TOL = RCONDV 00796 ELSE 00797 TOL = V / RCONDE 00798 END IF 00799 IF( V.GT.RCDVIN*RCDEIN ) THEN 00800 TOLIN = RCDVIN 00801 ELSE 00802 TOLIN = V / RCDEIN 00803 END IF 00804 TOL = MAX( TOL, SMLNUM / EPS ) 00805 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00806 IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN 00807 RESULT( 17 ) = ULPINV 00808 ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN 00809 RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL ) 00810 ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN 00811 RESULT( 17 ) = ULPINV 00812 ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN 00813 RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN ) 00814 ELSE 00815 RESULT( 17 ) = ONE 00816 END IF 00817 * 00818 270 CONTINUE 00819 * 00820 END IF 00821 * 00822 9999 FORMAT( ' CGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00823 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 ) 00824 9998 FORMAT( ' CGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00825 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00826 * 00827 RETURN 00828 * 00829 * End of CGET24 00830 * 00831 END
1.7.3 
