![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZDRGVX 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 ZDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI, 00012 * ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE, 00013 * S, DTRU, DIF, DIFTRU, WORK, LWORK, RWORK, 00014 * IWORK, LIWORK, RESULT, BWORK, INFO ) 00015 * 00016 * .. Scalar Arguments .. 00017 * INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT, 00018 * $ NSIZE 00019 * DOUBLE PRECISION THRESH 00020 * .. 00021 * .. Array Arguments .. 00022 * LOGICAL BWORK( * ) 00023 * INTEGER IWORK( * ) 00024 * DOUBLE PRECISION DIF( * ), DIFTRU( * ), DTRU( * ), LSCALE( * ), 00025 * $ RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * ) 00026 * COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ), 00027 * $ B( LDA, * ), BETA( * ), BI( LDA, * ), 00028 * $ VL( LDA, * ), VR( LDA, * ), WORK( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> ZDRGVX checks the nonsymmetric generalized eigenvalue problem 00038 *> expert driver ZGGEVX. 00039 *> 00040 *> ZGGEVX computes the generalized eigenvalues, (optionally) the left 00041 *> and/or right eigenvectors, (optionally) computes a balancing 00042 *> transformation to improve the conditioning, and (optionally) 00043 *> reciprocal condition numbers for the eigenvalues and eigenvectors. 00044 *> 00045 *> When ZDRGVX is called with NSIZE > 0, two types of test matrix pairs 00046 *> are generated by the subroutine DLATM6 and test the driver ZGGEVX. 00047 *> The test matrices have the known exact condition numbers for 00048 *> eigenvalues. For the condition numbers of the eigenvectors 00049 *> corresponding the first and last eigenvalues are also know 00050 *> ``exactly'' (see ZLATM6). 00051 *> For each matrix pair, the following tests will be performed and 00052 *> compared with the threshhold THRESH. 00053 *> 00054 *> (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of 00055 *> 00056 *> | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) ) 00057 *> 00058 *> where l**H is the conjugate tranpose of l. 00059 *> 00060 *> (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of 00061 *> 00062 *> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) ) 00063 *> 00064 *> (3) The condition number S(i) of eigenvalues computed by ZGGEVX 00065 *> differs less than a factor THRESH from the exact S(i) (see 00066 *> ZLATM6). 00067 *> 00068 *> (4) DIF(i) computed by ZTGSNA differs less than a factor 10*THRESH 00069 *> from the exact value (for the 1st and 5th vectors only). 00070 *> 00071 *> Test Matrices 00072 *> ============= 00073 *> 00074 *> Two kinds of test matrix pairs 00075 *> (A, B) = inverse(YH) * (Da, Db) * inverse(X) 00076 *> are used in the tests: 00077 *> 00078 *> 1: Da = 1+a 0 0 0 0 Db = 1 0 0 0 0 00079 *> 0 2+a 0 0 0 0 1 0 0 0 00080 *> 0 0 3+a 0 0 0 0 1 0 0 00081 *> 0 0 0 4+a 0 0 0 0 1 0 00082 *> 0 0 0 0 5+a , 0 0 0 0 1 , and 00083 *> 00084 *> 2: Da = 1 -1 0 0 0 Db = 1 0 0 0 0 00085 *> 1 1 0 0 0 0 1 0 0 0 00086 *> 0 0 1 0 0 0 0 1 0 0 00087 *> 0 0 0 1+a 1+b 0 0 0 1 0 00088 *> 0 0 0 -1-b 1+a , 0 0 0 0 1 . 00089 *> 00090 *> In both cases the same inverse(YH) and inverse(X) are used to compute 00091 *> (A, B), giving the exact eigenvectors to (A,B) as (YH, X): 00092 *> 00093 *> YH: = 1 0 -y y -y X = 1 0 -x -x x 00094 *> 0 1 -y y -y 0 1 x -x -x 00095 *> 0 0 1 0 0 0 0 1 0 0 00096 *> 0 0 0 1 0 0 0 0 1 0 00097 *> 0 0 0 0 1, 0 0 0 0 1 , where 00098 *> 00099 *> a, b, x and y will have all values independently of each other from 00100 *> { sqrt(sqrt(ULP)), 0.1, 1, 10, 1/sqrt(sqrt(ULP)) }. 00101 *> \endverbatim 00102 * 00103 * Arguments: 00104 * ========== 00105 * 00106 *> \param[in] NSIZE 00107 *> \verbatim 00108 *> NSIZE is INTEGER 00109 *> The number of sizes of matrices to use. NSIZE must be at 00110 *> least zero. If it is zero, no randomly generated matrices 00111 *> are tested, but any test matrices read from NIN will be 00112 *> tested. If it is not zero, then N = 5. 00113 *> \endverbatim 00114 *> 00115 *> \param[in] THRESH 00116 *> \verbatim 00117 *> THRESH is DOUBLE PRECISION 00118 *> A test will count as "failed" if the "error", computed as 00119 *> described above, exceeds THRESH. Note that the error 00120 *> is scaled to be O(1), so THRESH should be a reasonably 00121 *> small multiple of 1, e.g., 10 or 100. In particular, 00122 *> it should not depend on the precision (single vs. double) 00123 *> or the size of the matrix. It must be at least zero. 00124 *> \endverbatim 00125 *> 00126 *> \param[in] NIN 00127 *> \verbatim 00128 *> NIN is INTEGER 00129 *> The FORTRAN unit number for reading in the data file of 00130 *> problems to solve. 00131 *> \endverbatim 00132 *> 00133 *> \param[in] NOUT 00134 *> \verbatim 00135 *> NOUT is INTEGER 00136 *> The FORTRAN unit number for printing out error messages 00137 *> (e.g., if a routine returns IINFO not equal to 0.) 00138 *> \endverbatim 00139 *> 00140 *> \param[out] A 00141 *> \verbatim 00142 *> A is COMPLEX*16 array, dimension (LDA, NSIZE) 00143 *> Used to hold the matrix whose eigenvalues are to be 00144 *> computed. On exit, A contains the last matrix actually used. 00145 *> \endverbatim 00146 *> 00147 *> \param[in] LDA 00148 *> \verbatim 00149 *> LDA is INTEGER 00150 *> The leading dimension of A, B, AI, BI, Ao, and Bo. 00151 *> It must be at least 1 and at least NSIZE. 00152 *> \endverbatim 00153 *> 00154 *> \param[out] B 00155 *> \verbatim 00156 *> B is COMPLEX*16 array, dimension (LDA, NSIZE) 00157 *> Used to hold the matrix whose eigenvalues are to be 00158 *> computed. On exit, B contains the last matrix actually used. 00159 *> \endverbatim 00160 *> 00161 *> \param[out] AI 00162 *> \verbatim 00163 *> AI is COMPLEX*16 array, dimension (LDA, NSIZE) 00164 *> Copy of A, modified by ZGGEVX. 00165 *> \endverbatim 00166 *> 00167 *> \param[out] BI 00168 *> \verbatim 00169 *> BI is COMPLEX*16 array, dimension (LDA, NSIZE) 00170 *> Copy of B, modified by ZGGEVX. 00171 *> \endverbatim 00172 *> 00173 *> \param[out] ALPHA 00174 *> \verbatim 00175 *> ALPHA is COMPLEX*16 array, dimension (NSIZE) 00176 *> \endverbatim 00177 *> 00178 *> \param[out] BETA 00179 *> \verbatim 00180 *> BETA is COMPLEX*16 array, dimension (NSIZE) 00181 *> 00182 *> On exit, ALPHA/BETA are the eigenvalues. 00183 *> \endverbatim 00184 *> 00185 *> \param[out] VL 00186 *> \verbatim 00187 *> VL is COMPLEX*16 array, dimension (LDA, NSIZE) 00188 *> VL holds the left eigenvectors computed by ZGGEVX. 00189 *> \endverbatim 00190 *> 00191 *> \param[out] VR 00192 *> \verbatim 00193 *> VR is COMPLEX*16 array, dimension (LDA, NSIZE) 00194 *> VR holds the right eigenvectors computed by ZGGEVX. 00195 *> \endverbatim 00196 *> 00197 *> \param[out] ILO 00198 *> \verbatim 00199 *> ILO is INTEGER 00200 *> \endverbatim 00201 *> 00202 *> \param[out] IHI 00203 *> \verbatim 00204 *> IHI is INTEGER 00205 *> \endverbatim 00206 *> 00207 *> \param[out] LSCALE 00208 *> \verbatim 00209 *> LSCALE is DOUBLE PRECISION array, dimension (N) 00210 *> \endverbatim 00211 *> 00212 *> \param[out] RSCALE 00213 *> \verbatim 00214 *> RSCALE is DOUBLE PRECISION array, dimension (N) 00215 *> \endverbatim 00216 *> 00217 *> \param[out] S 00218 *> \verbatim 00219 *> S is DOUBLE PRECISION array, dimension (N) 00220 *> \endverbatim 00221 *> 00222 *> \param[out] DTRU 00223 *> \verbatim 00224 *> DTRU is DOUBLE PRECISION array, dimension (N) 00225 *> \endverbatim 00226 *> 00227 *> \param[out] DIF 00228 *> \verbatim 00229 *> DIF is DOUBLE PRECISION array, dimension (N) 00230 *> \endverbatim 00231 *> 00232 *> \param[out] DIFTRU 00233 *> \verbatim 00234 *> DIFTRU is DOUBLE PRECISION array, dimension (N) 00235 *> \endverbatim 00236 *> 00237 *> \param[out] WORK 00238 *> \verbatim 00239 *> WORK is COMPLEX*16 array, dimension (LWORK) 00240 *> \endverbatim 00241 *> 00242 *> \param[in] LWORK 00243 *> \verbatim 00244 *> LWORK is INTEGER 00245 *> Leading dimension of WORK. LWORK >= 2*N*N + 2*N 00246 *> \endverbatim 00247 *> 00248 *> \param[out] RWORK 00249 *> \verbatim 00250 *> RWORK is DOUBLE PRECISION array, dimension (6*N) 00251 *> \endverbatim 00252 *> 00253 *> \param[out] IWORK 00254 *> \verbatim 00255 *> IWORK is INTEGER array, dimension (LIWORK) 00256 *> \endverbatim 00257 *> 00258 *> \param[in] LIWORK 00259 *> \verbatim 00260 *> LIWORK is INTEGER 00261 *> Leading dimension of IWORK. LIWORK >= N+2. 00262 *> \endverbatim 00263 *> 00264 *> \param[out] RESULT 00265 *> \verbatim 00266 *> RESULT is DOUBLE PRECISION array, dimension (4) 00267 *> \endverbatim 00268 *> 00269 *> \param[out] BWORK 00270 *> \verbatim 00271 *> BWORK is LOGICAL array, dimension (N) 00272 *> \endverbatim 00273 *> 00274 *> \param[out] INFO 00275 *> \verbatim 00276 *> INFO is INTEGER 00277 *> = 0: successful exit 00278 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00279 *> > 0: A routine returned an error code. 00280 *> \endverbatim 00281 * 00282 * Authors: 00283 * ======== 00284 * 00285 *> \author Univ. of Tennessee 00286 *> \author Univ. of California Berkeley 00287 *> \author Univ. of Colorado Denver 00288 *> \author NAG Ltd. 00289 * 00290 *> \date November 2011 00291 * 00292 *> \ingroup complex16_eig 00293 * 00294 * ===================================================================== 00295 SUBROUTINE ZDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI, 00296 $ ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE, 00297 $ S, DTRU, DIF, DIFTRU, WORK, LWORK, RWORK, 00298 $ IWORK, LIWORK, RESULT, BWORK, INFO ) 00299 * 00300 * -- LAPACK test routine (version 3.4.0) -- 00301 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00302 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00303 * November 2011 00304 * 00305 * .. Scalar Arguments .. 00306 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT, 00307 $ NSIZE 00308 DOUBLE PRECISION THRESH 00309 * .. 00310 * .. Array Arguments .. 00311 LOGICAL BWORK( * ) 00312 INTEGER IWORK( * ) 00313 DOUBLE PRECISION DIF( * ), DIFTRU( * ), DTRU( * ), LSCALE( * ), 00314 $ RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * ) 00315 COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ), 00316 $ B( LDA, * ), BETA( * ), BI( LDA, * ), 00317 $ VL( LDA, * ), VR( LDA, * ), WORK( * ) 00318 * .. 00319 * 00320 * ===================================================================== 00321 * 00322 * .. Parameters .. 00323 DOUBLE PRECISION ZERO, ONE, TEN, TNTH, HALF 00324 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1, 00325 $ TNTH = 1.0D-1, HALF = 0.5D+0 ) 00326 * .. 00327 * .. Local Scalars .. 00328 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO, 00329 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT 00330 DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2, 00331 $ ULP, ULPINV 00332 * .. 00333 * .. Local Arrays .. 00334 COMPLEX*16 WEIGHT( 5 ) 00335 * .. 00336 * .. External Functions .. 00337 INTEGER ILAENV 00338 DOUBLE PRECISION DLAMCH, ZLANGE 00339 EXTERNAL ILAENV, DLAMCH, ZLANGE 00340 * .. 00341 * .. External Subroutines .. 00342 EXTERNAL ALASVM, XERBLA, ZGET52, ZGGEVX, ZLACPY, ZLATM6 00343 * .. 00344 * .. Intrinsic Functions .. 00345 INTRINSIC ABS, DCMPLX, MAX, SQRT 00346 * .. 00347 * .. Executable Statements .. 00348 * 00349 * Check for errors 00350 * 00351 INFO = 0 00352 * 00353 NMAX = 5 00354 * 00355 IF( NSIZE.LT.0 ) THEN 00356 INFO = -1 00357 ELSE IF( THRESH.LT.ZERO ) THEN 00358 INFO = -2 00359 ELSE IF( NIN.LE.0 ) THEN 00360 INFO = -3 00361 ELSE IF( NOUT.LE.0 ) THEN 00362 INFO = -4 00363 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN 00364 INFO = -6 00365 ELSE IF( LIWORK.LT.NMAX+2 ) THEN 00366 INFO = -26 00367 END IF 00368 * 00369 * Compute workspace 00370 * (Note: Comments in the code beginning "Workspace:" describe the 00371 * minimal amount of workspace needed at that point in the code, 00372 * as well as the preferred amount for good performance. 00373 * NB refers to the optimal block size for the immediately 00374 * following subroutine, as returned by ILAENV.) 00375 * 00376 MINWRK = 1 00377 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 00378 MINWRK = 2*NMAX*( NMAX+1 ) 00379 MAXWRK = NMAX*( 1+ILAENV( 1, 'ZGEQRF', ' ', NMAX, 1, NMAX, 00380 $ 0 ) ) 00381 MAXWRK = MAX( MAXWRK, 2*NMAX*( NMAX+1 ) ) 00382 WORK( 1 ) = MAXWRK 00383 END IF 00384 * 00385 IF( LWORK.LT.MINWRK ) 00386 $ INFO = -23 00387 * 00388 IF( INFO.NE.0 ) THEN 00389 CALL XERBLA( 'ZDRGVX', -INFO ) 00390 RETURN 00391 END IF 00392 * 00393 N = 5 00394 ULP = DLAMCH( 'P' ) 00395 ULPINV = ONE / ULP 00396 THRSH2 = TEN*THRESH 00397 NERRS = 0 00398 NPTKNT = 0 00399 NTESTT = 0 00400 * 00401 IF( NSIZE.EQ.0 ) 00402 $ GO TO 90 00403 * 00404 * Parameters used for generating test matrices. 00405 * 00406 WEIGHT( 1 ) = DCMPLX( TNTH, ZERO ) 00407 WEIGHT( 2 ) = DCMPLX( HALF, ZERO ) 00408 WEIGHT( 3 ) = ONE 00409 WEIGHT( 4 ) = ONE / WEIGHT( 2 ) 00410 WEIGHT( 5 ) = ONE / WEIGHT( 1 ) 00411 * 00412 DO 80 IPTYPE = 1, 2 00413 DO 70 IWA = 1, 5 00414 DO 60 IWB = 1, 5 00415 DO 50 IWX = 1, 5 00416 DO 40 IWY = 1, 5 00417 * 00418 * generated a pair of test matrix 00419 * 00420 CALL ZLATM6( IPTYPE, 5, A, LDA, B, VR, LDA, VL, 00421 $ LDA, WEIGHT( IWA ), WEIGHT( IWB ), 00422 $ WEIGHT( IWX ), WEIGHT( IWY ), DTRU, 00423 $ DIFTRU ) 00424 * 00425 * Compute eigenvalues/eigenvectors of (A, B). 00426 * Compute eigenvalue/eigenvector condition numbers 00427 * using computed eigenvectors. 00428 * 00429 CALL ZLACPY( 'F', N, N, A, LDA, AI, LDA ) 00430 CALL ZLACPY( 'F', N, N, B, LDA, BI, LDA ) 00431 * 00432 CALL ZGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI, 00433 $ LDA, ALPHA, BETA, VL, LDA, VR, LDA, 00434 $ ILO, IHI, LSCALE, RSCALE, ANORM, 00435 $ BNORM, S, DIF, WORK, LWORK, RWORK, 00436 $ IWORK, BWORK, LINFO ) 00437 IF( LINFO.NE.0 ) THEN 00438 WRITE( NOUT, FMT = 9999 )'ZGGEVX', LINFO, N, 00439 $ IPTYPE, IWA, IWB, IWX, IWY 00440 GO TO 30 00441 END IF 00442 * 00443 * Compute the norm(A, B) 00444 * 00445 CALL ZLACPY( 'Full', N, N, AI, LDA, WORK, N ) 00446 CALL ZLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ), 00447 $ N ) 00448 ABNORM = ZLANGE( 'Fro', N, 2*N, WORK, N, RWORK ) 00449 * 00450 * Tests (1) and (2) 00451 * 00452 RESULT( 1 ) = ZERO 00453 CALL ZGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, 00454 $ ALPHA, BETA, WORK, RWORK, 00455 $ RESULT( 1 ) ) 00456 IF( RESULT( 2 ).GT.THRESH ) THEN 00457 WRITE( NOUT, FMT = 9998 )'Left', 'ZGGEVX', 00458 $ RESULT( 2 ), N, IPTYPE, IWA, IWB, IWX, IWY 00459 END IF 00460 * 00461 RESULT( 2 ) = ZERO 00462 CALL ZGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, 00463 $ ALPHA, BETA, WORK, RWORK, 00464 $ RESULT( 2 ) ) 00465 IF( RESULT( 3 ).GT.THRESH ) THEN 00466 WRITE( NOUT, FMT = 9998 )'Right', 'ZGGEVX', 00467 $ RESULT( 3 ), N, IPTYPE, IWA, IWB, IWX, IWY 00468 END IF 00469 * 00470 * Test (3) 00471 * 00472 RESULT( 3 ) = ZERO 00473 DO 10 I = 1, N 00474 IF( S( I ).EQ.ZERO ) THEN 00475 IF( DTRU( I ).GT.ABNORM*ULP ) 00476 $ RESULT( 3 ) = ULPINV 00477 ELSE IF( DTRU( I ).EQ.ZERO ) THEN 00478 IF( S( I ).GT.ABNORM*ULP ) 00479 $ RESULT( 3 ) = ULPINV 00480 ELSE 00481 RWORK( I ) = MAX( ABS( DTRU( I ) / S( I ) ), 00482 $ ABS( S( I ) / DTRU( I ) ) ) 00483 RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) ) 00484 END IF 00485 10 CONTINUE 00486 * 00487 * Test (4) 00488 * 00489 RESULT( 4 ) = ZERO 00490 IF( DIF( 1 ).EQ.ZERO ) THEN 00491 IF( DIFTRU( 1 ).GT.ABNORM*ULP ) 00492 $ RESULT( 4 ) = ULPINV 00493 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN 00494 IF( DIF( 1 ).GT.ABNORM*ULP ) 00495 $ RESULT( 4 ) = ULPINV 00496 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN 00497 IF( DIFTRU( 5 ).GT.ABNORM*ULP ) 00498 $ RESULT( 4 ) = ULPINV 00499 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN 00500 IF( DIF( 5 ).GT.ABNORM*ULP ) 00501 $ RESULT( 4 ) = ULPINV 00502 ELSE 00503 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ), 00504 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) ) 00505 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ), 00506 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) ) 00507 RESULT( 4 ) = MAX( RATIO1, RATIO2 ) 00508 END IF 00509 * 00510 NTESTT = NTESTT + 4 00511 * 00512 * Print out tests which fail. 00513 * 00514 DO 20 J = 1, 4 00515 IF( ( RESULT( J ).GE.THRSH2 .AND. J.GE.4 ) .OR. 00516 $ ( RESULT( J ).GE.THRESH .AND. J.LE.3 ) ) 00517 $ THEN 00518 * 00519 * If this is the first test to fail, 00520 * print a header to the data file. 00521 * 00522 IF( NERRS.EQ.0 ) THEN 00523 WRITE( NOUT, FMT = 9997 )'ZXV' 00524 * 00525 * Print out messages for built-in examples 00526 * 00527 * Matrix types 00528 * 00529 WRITE( NOUT, FMT = 9995 ) 00530 WRITE( NOUT, FMT = 9994 ) 00531 WRITE( NOUT, FMT = 9993 ) 00532 * 00533 * Tests performed 00534 * 00535 WRITE( NOUT, FMT = 9992 )'''', 00536 $ 'transpose', '''' 00537 * 00538 END IF 00539 NERRS = NERRS + 1 00540 IF( RESULT( J ).LT.10000.0D0 ) THEN 00541 WRITE( NOUT, FMT = 9991 )IPTYPE, IWA, 00542 $ IWB, IWX, IWY, J, RESULT( J ) 00543 ELSE 00544 WRITE( NOUT, FMT = 9990 )IPTYPE, IWA, 00545 $ IWB, IWX, IWY, J, RESULT( J ) 00546 END IF 00547 END IF 00548 20 CONTINUE 00549 * 00550 30 CONTINUE 00551 * 00552 40 CONTINUE 00553 50 CONTINUE 00554 60 CONTINUE 00555 70 CONTINUE 00556 80 CONTINUE 00557 * 00558 GO TO 150 00559 * 00560 90 CONTINUE 00561 * 00562 * Read in data from file to check accuracy of condition estimation 00563 * Read input data until N=0 00564 * 00565 READ( NIN, FMT = *, END = 150 )N 00566 IF( N.EQ.0 ) 00567 $ GO TO 150 00568 DO 100 I = 1, N 00569 READ( NIN, FMT = * )( A( I, J ), J = 1, N ) 00570 100 CONTINUE 00571 DO 110 I = 1, N 00572 READ( NIN, FMT = * )( B( I, J ), J = 1, N ) 00573 110 CONTINUE 00574 READ( NIN, FMT = * )( DTRU( I ), I = 1, N ) 00575 READ( NIN, FMT = * )( DIFTRU( I ), I = 1, N ) 00576 * 00577 NPTKNT = NPTKNT + 1 00578 * 00579 * Compute eigenvalues/eigenvectors of (A, B). 00580 * Compute eigenvalue/eigenvector condition numbers 00581 * using computed eigenvectors. 00582 * 00583 CALL ZLACPY( 'F', N, N, A, LDA, AI, LDA ) 00584 CALL ZLACPY( 'F', N, N, B, LDA, BI, LDA ) 00585 * 00586 CALL ZGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI, LDA, ALPHA, BETA, 00587 $ VL, LDA, VR, LDA, ILO, IHI, LSCALE, RSCALE, ANORM, 00588 $ BNORM, S, DIF, WORK, LWORK, RWORK, IWORK, BWORK, 00589 $ LINFO ) 00590 * 00591 IF( LINFO.NE.0 ) THEN 00592 WRITE( NOUT, FMT = 9987 )'ZGGEVX', LINFO, N, NPTKNT 00593 GO TO 140 00594 END IF 00595 * 00596 * Compute the norm(A, B) 00597 * 00598 CALL ZLACPY( 'Full', N, N, AI, LDA, WORK, N ) 00599 CALL ZLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ), N ) 00600 ABNORM = ZLANGE( 'Fro', N, 2*N, WORK, N, RWORK ) 00601 * 00602 * Tests (1) and (2) 00603 * 00604 RESULT( 1 ) = ZERO 00605 CALL ZGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, ALPHA, BETA, 00606 $ WORK, RWORK, RESULT( 1 ) ) 00607 IF( RESULT( 2 ).GT.THRESH ) THEN 00608 WRITE( NOUT, FMT = 9986 )'Left', 'ZGGEVX', RESULT( 2 ), N, 00609 $ NPTKNT 00610 END IF 00611 * 00612 RESULT( 2 ) = ZERO 00613 CALL ZGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, ALPHA, BETA, 00614 $ WORK, RWORK, RESULT( 2 ) ) 00615 IF( RESULT( 3 ).GT.THRESH ) THEN 00616 WRITE( NOUT, FMT = 9986 )'Right', 'ZGGEVX', RESULT( 3 ), N, 00617 $ NPTKNT 00618 END IF 00619 * 00620 * Test (3) 00621 * 00622 RESULT( 3 ) = ZERO 00623 DO 120 I = 1, N 00624 IF( S( I ).EQ.ZERO ) THEN 00625 IF( DTRU( I ).GT.ABNORM*ULP ) 00626 $ RESULT( 3 ) = ULPINV 00627 ELSE IF( DTRU( I ).EQ.ZERO ) THEN 00628 IF( S( I ).GT.ABNORM*ULP ) 00629 $ RESULT( 3 ) = ULPINV 00630 ELSE 00631 RWORK( I ) = MAX( ABS( DTRU( I ) / S( I ) ), 00632 $ ABS( S( I ) / DTRU( I ) ) ) 00633 RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) ) 00634 END IF 00635 120 CONTINUE 00636 * 00637 * Test (4) 00638 * 00639 RESULT( 4 ) = ZERO 00640 IF( DIF( 1 ).EQ.ZERO ) THEN 00641 IF( DIFTRU( 1 ).GT.ABNORM*ULP ) 00642 $ RESULT( 4 ) = ULPINV 00643 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN 00644 IF( DIF( 1 ).GT.ABNORM*ULP ) 00645 $ RESULT( 4 ) = ULPINV 00646 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN 00647 IF( DIFTRU( 5 ).GT.ABNORM*ULP ) 00648 $ RESULT( 4 ) = ULPINV 00649 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN 00650 IF( DIF( 5 ).GT.ABNORM*ULP ) 00651 $ RESULT( 4 ) = ULPINV 00652 ELSE 00653 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ), 00654 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) ) 00655 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ), 00656 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) ) 00657 RESULT( 4 ) = MAX( RATIO1, RATIO2 ) 00658 END IF 00659 * 00660 NTESTT = NTESTT + 4 00661 * 00662 * Print out tests which fail. 00663 * 00664 DO 130 J = 1, 4 00665 IF( RESULT( J ).GE.THRSH2 ) THEN 00666 * 00667 * If this is the first test to fail, 00668 * print a header to the data file. 00669 * 00670 IF( NERRS.EQ.0 ) THEN 00671 WRITE( NOUT, FMT = 9997 )'ZXV' 00672 * 00673 * Print out messages for built-in examples 00674 * 00675 * Matrix types 00676 * 00677 WRITE( NOUT, FMT = 9996 ) 00678 * 00679 * Tests performed 00680 * 00681 WRITE( NOUT, FMT = 9992 )'''', 'transpose', '''' 00682 * 00683 END IF 00684 NERRS = NERRS + 1 00685 IF( RESULT( J ).LT.10000.0D0 ) THEN 00686 WRITE( NOUT, FMT = 9989 )NPTKNT, N, J, RESULT( J ) 00687 ELSE 00688 WRITE( NOUT, FMT = 9988 )NPTKNT, N, J, RESULT( J ) 00689 END IF 00690 END IF 00691 130 CONTINUE 00692 * 00693 140 CONTINUE 00694 * 00695 GO TO 90 00696 150 CONTINUE 00697 * 00698 * Summary 00699 * 00700 CALL ALASVM( 'ZXV', NOUT, NERRS, NTESTT, 0 ) 00701 * 00702 WORK( 1 ) = MAXWRK 00703 * 00704 RETURN 00705 * 00706 9999 FORMAT( ' ZDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00707 $ I6, ', JTYPE=', I6, ')' ) 00708 * 00709 9998 FORMAT( ' ZDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ', 00710 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 00711 $ 'N=', I6, ', JTYPE=', I6, ', IWA=', I5, ', IWB=', I5, 00712 $ ', IWX=', I5, ', IWY=', I5 ) 00713 * 00714 9997 FORMAT( / 1X, A3, ' -- Complex Expert Eigenvalue/vector', 00715 $ ' problem driver' ) 00716 * 00717 9996 FORMAT( 'Input Example' ) 00718 * 00719 9995 FORMAT( ' Matrix types: ', / ) 00720 * 00721 9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ', 00722 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ', 00723 $ / ' YH and X are left and right eigenvectors. ', / ) 00724 * 00725 9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ', 00726 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ', 00727 $ / ' YH and X are left and right eigenvectors. ', / ) 00728 * 00729 9992 FORMAT( / ' Tests performed: ', / 4X, 00730 $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4X, 00731 $ ' r is a right eigenvector and ', A, ' means ', A, '.', 00732 $ / ' 1 = max | ( b A - a B )', A, ' l | / const.', 00733 $ / ' 2 = max | ( b A - a B ) r | / const.', 00734 $ / ' 3 = max ( Sest/Stru, Stru/Sest ) ', 00735 $ ' over all eigenvalues', / 00736 $ ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ', 00737 $ ' over the 1st and 5th eigenvectors', / ) 00738 * 00739 9991 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=', 00740 $ I2, ', IWY=', I2, ', result ', I2, ' is', 0P, F8.2 ) 00741 * 00742 9990 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=', 00743 $ I2, ', IWY=', I2, ', result ', I2, ' is', 1P, D10.3 ) 00744 * 00745 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00746 $ ' result ', I2, ' is', 0P, F8.2 ) 00747 * 00748 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00749 $ ' result ', I2, ' is', 1P, D10.3 ) 00750 * 00751 9987 FORMAT( ' ZDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00752 $ I6, ', Input example #', I2, ')' ) 00753 * 00754 9986 FORMAT( ' ZDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ', 00755 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 00756 $ 'N=', I6, ', Input Example #', I2, ')' ) 00757 * 00758 * End of ZDRGVX 00759 * 00760 END