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