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