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