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