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