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