 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGET23( COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N, 00002 $ A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR, 00003 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN, 00004 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT, 00005 $ WORK, LWORK, IWORK, INFO ) 00006 * 00007 * -- LAPACK test routine (version 3.1) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * November 2006 00010 * 00011 * .. Scalar Arguments .. 00012 LOGICAL COMP 00013 CHARACTER BALANC 00014 INTEGER INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N, 00015 $ NOUNIT 00016 DOUBLE PRECISION THRESH 00017 * .. 00018 * .. Array Arguments .. 00019 INTEGER ISEED( 4 ), IWORK( * ) 00020 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ), 00021 $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ), 00022 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ), 00023 $ RESULT( 11 ), SCALE( * ), SCALE1( * ), 00024 $ VL( LDVL, * ), VR( LDVR, * ), WI( * ), 00025 $ WI1( * ), WORK( * ), WR( * ), WR1( * ) 00026 * .. 00027 * 00028 * Purpose 00029 * ======= 00030 * 00031 * DGET23 checks the nonsymmetric eigenvalue problem driver SGEEVX. 00032 * If COMP = .FALSE., the first 8 of the following tests will be 00033 * performed on the input matrix A, and also test 9 if LWORK is 00034 * sufficiently large. 00035 * if COMP is .TRUE. all 11 tests will be performed. 00036 * 00037 * (1) | A * VR - VR * W | / ( n |A| ulp ) 00038 * 00039 * Here VR is the matrix of unit right eigenvectors. 00040 * W is a block diagonal matrix, with a 1x1 block for each 00041 * real eigenvalue and a 2x2 block for each complex conjugate 00042 * pair. If eigenvalues j and j+1 are a complex conjugate pair, 00043 * so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 00044 * 2 x 2 block corresponding to the pair will be: 00045 * 00046 * ( wr wi ) 00047 * ( -wi wr ) 00048 * 00049 * Such a block multiplying an n x 2 matrix ( ur ui ) on the 00050 * right will be the same as multiplying ur + i*ui by wr + i*wi. 00051 * 00052 * (2) | A**H * VL - VL * W**H | / ( n |A| ulp ) 00053 * 00054 * Here VL is the matrix of unit left eigenvectors, A**H is the 00055 * conjugate transpose of A, and W is as above. 00056 * 00057 * (3) | |VR(i)| - 1 | / ulp and largest component real 00058 * 00059 * VR(i) denotes the i-th column of VR. 00060 * 00061 * (4) | |VL(i)| - 1 | / ulp and largest component real 00062 * 00063 * VL(i) denotes the i-th column of VL. 00064 * 00065 * (5) 0 if W(full) = W(partial), 1/ulp otherwise 00066 * 00067 * W(full) denotes the eigenvalues computed when VR, VL, RCONDV 00068 * and RCONDE are also computed, and W(partial) denotes the 00069 * eigenvalues computed when only some of VR, VL, RCONDV, and 00070 * RCONDE are computed. 00071 * 00072 * (6) 0 if VR(full) = VR(partial), 1/ulp otherwise 00073 * 00074 * VR(full) denotes the right eigenvectors computed when VL, RCONDV 00075 * and RCONDE are computed, and VR(partial) denotes the result 00076 * when only some of VL and RCONDV are computed. 00077 * 00078 * (7) 0 if VL(full) = VL(partial), 1/ulp otherwise 00079 * 00080 * VL(full) denotes the left eigenvectors computed when VR, RCONDV 00081 * and RCONDE are computed, and VL(partial) denotes the result 00082 * when only some of VR and RCONDV are computed. 00083 * 00084 * (8) 0 if SCALE, ILO, IHI, ABNRM (full) = 00085 * SCALE, ILO, IHI, ABNRM (partial) 00086 * 1/ulp otherwise 00087 * 00088 * SCALE, ILO, IHI and ABNRM describe how the matrix is balanced. 00089 * (full) is when VR, VL, RCONDE and RCONDV are also computed, and 00090 * (partial) is when some are not computed. 00091 * 00092 * (9) 0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise 00093 * 00094 * RCONDV(full) denotes the reciprocal condition numbers of the 00095 * right eigenvectors computed when VR, VL and RCONDE are also 00096 * computed. RCONDV(partial) denotes the reciprocal condition 00097 * numbers when only some of VR, VL and RCONDE are computed. 00098 * 00099 * (10) |RCONDV - RCDVIN| / cond(RCONDV) 00100 * 00101 * RCONDV is the reciprocal right eigenvector condition number 00102 * computed by DGEEVX and RCDVIN (the precomputed true value) 00103 * is supplied as input. cond(RCONDV) is the condition number of 00104 * RCONDV, and takes errors in computing RCONDV into account, so 00105 * that the resulting quantity should be O(ULP). cond(RCONDV) is 00106 * essentially given by norm(A)/RCONDE. 00107 * 00108 * (11) |RCONDE - RCDEIN| / cond(RCONDE) 00109 * 00110 * RCONDE is the reciprocal eigenvalue condition number 00111 * computed by DGEEVX and RCDEIN (the precomputed true value) 00112 * is supplied as input. cond(RCONDE) is the condition number 00113 * of RCONDE, and takes errors in computing RCONDE into account, 00114 * so that the resulting quantity should be O(ULP). cond(RCONDE) 00115 * is essentially given by norm(A)/RCONDV. 00116 * 00117 * Arguments 00118 * ========= 00119 * 00120 * COMP (input) LOGICAL 00121 * COMP describes which input tests to perform: 00122 * = .FALSE. if the computed condition numbers are not to 00123 * be tested against RCDVIN and RCDEIN 00124 * = .TRUE. if they are to be compared 00125 * 00126 * BALANC (input) CHARACTER 00127 * Describes the balancing option to be tested. 00128 * = 'N' for no permuting or diagonal scaling 00129 * = 'P' for permuting but no diagonal scaling 00130 * = 'S' for no permuting but diagonal scaling 00131 * = 'B' for permuting and diagonal scaling 00132 * 00133 * JTYPE (input) INTEGER 00134 * Type of input matrix. Used to label output if error occurs. 00135 * 00136 * THRESH (input) DOUBLE PRECISION 00137 * A test will count as "failed" if the "error", computed as 00138 * described above, exceeds THRESH. Note that the error 00139 * is scaled to be O(1), so THRESH should be a reasonably 00140 * small multiple of 1, e.g., 10 or 100. In particular, 00141 * it should not depend on the precision (single vs. double) 00142 * or the size of the matrix. It must be at least zero. 00143 * 00144 * ISEED (input) INTEGER array, dimension (4) 00145 * If COMP = .FALSE., the random number generator seed 00146 * used to produce matrix. 00147 * If COMP = .TRUE., ISEED(1) = the number of the example. 00148 * Used to label output if error occurs. 00149 * 00150 * NOUNIT (input) INTEGER 00151 * The FORTRAN unit number for printing out error messages 00152 * (e.g., if a routine returns INFO not equal to 0.) 00153 * 00154 * N (input) INTEGER 00155 * The dimension of A. N must be at least 0. 00156 * 00157 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00158 * Used to hold the matrix whose eigenvalues are to be 00159 * computed. 00160 * 00161 * LDA (input) INTEGER 00162 * The leading dimension of A, and H. LDA must be at 00163 * least 1 and at least N. 00164 * 00165 * H (workspace) DOUBLE PRECISION array, dimension (LDA,N) 00166 * Another copy of the test matrix A, modified by DGEEVX. 00167 * 00168 * WR (workspace) DOUBLE PRECISION array, dimension (N) 00169 * WI (workspace) DOUBLE PRECISION array, dimension (N) 00170 * The real and imaginary parts of the eigenvalues of A. 00171 * On exit, WR + WI*i are the eigenvalues of the matrix in A. 00172 * 00173 * WR1 (workspace) DOUBLE PRECISION array, dimension (N) 00174 * WI1 (workspace) DOUBLE PRECISION array, dimension (N) 00175 * Like WR, WI, these arrays contain the eigenvalues of A, 00176 * but those computed when DGEEVX only computes a partial 00177 * eigendecomposition, i.e. not the eigenvalues and left 00178 * and right eigenvectors. 00179 * 00180 * VL (workspace) DOUBLE PRECISION array, dimension (LDVL,N) 00181 * VL holds the computed left eigenvectors. 00182 * 00183 * LDVL (input) INTEGER 00184 * Leading dimension of VL. Must be at least max(1,N). 00185 * 00186 * VR (workspace) DOUBLE PRECISION array, dimension (LDVR,N) 00187 * VR holds the computed right eigenvectors. 00188 * 00189 * LDVR (input) INTEGER 00190 * Leading dimension of VR. Must be at least max(1,N). 00191 * 00192 * LRE (workspace) DOUBLE PRECISION array, dimension (LDLRE,N) 00193 * LRE holds the computed right or left eigenvectors. 00194 * 00195 * LDLRE (input) INTEGER 00196 * Leading dimension of LRE. Must be at least max(1,N). 00197 * 00198 * RCONDV (workspace) DOUBLE PRECISION array, dimension (N) 00199 * RCONDV holds the computed reciprocal condition numbers 00200 * for eigenvectors. 00201 * 00202 * RCNDV1 (workspace) DOUBLE PRECISION array, dimension (N) 00203 * RCNDV1 holds more computed reciprocal condition numbers 00204 * for eigenvectors. 00205 * 00206 * RCDVIN (input) DOUBLE PRECISION array, dimension (N) 00207 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal 00208 * condition numbers for eigenvectors to be compared with 00209 * RCONDV. 00210 * 00211 * RCONDE (workspace) DOUBLE PRECISION array, dimension (N) 00212 * RCONDE holds the computed reciprocal condition numbers 00213 * for eigenvalues. 00214 * 00215 * RCNDE1 (workspace) DOUBLE PRECISION array, dimension (N) 00216 * RCNDE1 holds more computed reciprocal condition numbers 00217 * for eigenvalues. 00218 * 00219 * RCDEIN (input) DOUBLE PRECISION array, dimension (N) 00220 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal 00221 * condition numbers for eigenvalues to be compared with 00222 * RCONDE. 00223 * 00224 * SCALE (workspace) DOUBLE PRECISION array, dimension (N) 00225 * Holds information describing balancing of matrix. 00226 * 00227 * SCALE1 (workspace) DOUBLE PRECISION array, dimension (N) 00228 * Holds information describing balancing of matrix. 00229 * 00230 * RESULT (output) DOUBLE PRECISION array, dimension (11) 00231 * The values computed by the 11 tests described above. 00232 * The values are currently limited to 1/ulp, to avoid 00233 * overflow. 00234 * 00235 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00236 * 00237 * LWORK (input) INTEGER 00238 * The number of entries in WORK. This must be at least 00239 * 3*N, and 6*N+N**2 if tests 9, 10 or 11 are to be performed. 00240 * 00241 * IWORK (workspace) INTEGER array, dimension (2*N) 00242 * 00243 * INFO (output) INTEGER 00244 * If 0, successful exit. 00245 * If <0, input parameter -INFO had an incorrect value. 00246 * If >0, DGEEVX returned an error code, the absolute 00247 * value of which is returned. 00248 * 00249 * ===================================================================== 00250 * 00251 * 00252 * .. Parameters .. 00253 DOUBLE PRECISION ZERO, ONE, TWO 00254 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00255 DOUBLE PRECISION EPSIN 00256 PARAMETER ( EPSIN = 5.9605D-8 ) 00257 * .. 00258 * .. Local Scalars .. 00259 LOGICAL BALOK, NOBAL 00260 CHARACTER SENSE 00261 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM, 00262 $ J, JJ, KMIN 00263 DOUBLE PRECISION ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN, 00264 $ ULP, ULPINV, V, VIMIN, VMAX, VMX, VRMIN, VRMX, 00265 $ VTST 00266 * .. 00267 * .. Local Arrays .. 00268 CHARACTER SENS( 2 ) 00269 DOUBLE PRECISION DUM( 1 ), RES( 2 ) 00270 * .. 00271 * .. External Functions .. 00272 LOGICAL LSAME 00273 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 00274 EXTERNAL LSAME, DLAMCH, DLAPY2, DNRM2 00275 * .. 00276 * .. External Subroutines .. 00277 EXTERNAL DGEEVX, DGET22, DLACPY, XERBLA 00278 * .. 00279 * .. Intrinsic Functions .. 00280 INTRINSIC ABS, DBLE, MAX, MIN 00281 * .. 00282 * .. Data statements .. 00283 DATA SENS / 'N', 'V' / 00284 * .. 00285 * .. Executable Statements .. 00286 * 00287 * Check for errors 00288 * 00289 NOBAL = LSAME( BALANC, 'N' ) 00290 BALOK = NOBAL .OR. LSAME( BALANC, 'P' ) .OR. 00291 $ LSAME( BALANC, 'S' ) .OR. LSAME( BALANC, 'B' ) 00292 INFO = 0 00293 IF( .NOT.BALOK ) THEN 00294 INFO = -2 00295 ELSE IF( THRESH.LT.ZERO ) THEN 00296 INFO = -4 00297 ELSE IF( NOUNIT.LE.0 ) THEN 00298 INFO = -6 00299 ELSE IF( N.LT.0 ) THEN 00300 INFO = -7 00301 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN 00302 INFO = -9 00303 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.N ) THEN 00304 INFO = -16 00305 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.N ) THEN 00306 INFO = -18 00307 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.N ) THEN 00308 INFO = -20 00309 ELSE IF( LWORK.LT.3*N .OR. ( COMP .AND. LWORK.LT.6*N+N*N ) ) THEN 00310 INFO = -31 00311 END IF 00312 * 00313 IF( INFO.NE.0 ) THEN 00314 CALL XERBLA( 'DGET23', -INFO ) 00315 RETURN 00316 END IF 00317 * 00318 * Quick return if nothing to do 00319 * 00320 DO 10 I = 1, 11 00321 RESULT( I ) = -ONE 00322 10 CONTINUE 00323 * 00324 IF( N.EQ.0 ) 00325 $ RETURN 00326 * 00327 * More Important constants 00328 * 00329 ULP = DLAMCH( 'Precision' ) 00330 SMLNUM = DLAMCH( 'S' ) 00331 ULPINV = ONE / ULP 00332 * 00333 * Compute eigenvalues and eigenvectors, and test them 00334 * 00335 IF( LWORK.GE.6*N+N*N ) THEN 00336 SENSE = 'B' 00337 ISENSM = 2 00338 ELSE 00339 SENSE = 'E' 00340 ISENSM = 1 00341 END IF 00342 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 00343 CALL DGEEVX( BALANC, 'V', 'V', SENSE, N, H, LDA, WR, WI, VL, LDVL, 00344 $ VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, 00345 $ WORK, LWORK, IWORK, IINFO ) 00346 IF( IINFO.NE.0 ) THEN 00347 RESULT( 1 ) = ULPINV 00348 IF( JTYPE.NE.22 ) THEN 00349 WRITE( NOUNIT, FMT = 9998 )'DGEEVX1', IINFO, N, JTYPE, 00350 $ BALANC, ISEED 00351 ELSE 00352 WRITE( NOUNIT, FMT = 9999 )'DGEEVX1', IINFO, N, ISEED( 1 ) 00353 END IF 00354 INFO = ABS( IINFO ) 00355 RETURN 00356 END IF 00357 * 00358 * Do Test (1) 00359 * 00360 CALL DGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, WR, WI, WORK, 00361 $ RES ) 00362 RESULT( 1 ) = RES( 1 ) 00363 * 00364 * Do Test (2) 00365 * 00366 CALL DGET22( 'T', 'N', 'T', N, A, LDA, VL, LDVL, WR, WI, WORK, 00367 $ RES ) 00368 RESULT( 2 ) = RES( 1 ) 00369 * 00370 * Do Test (3) 00371 * 00372 DO 30 J = 1, N 00373 TNRM = ONE 00374 IF( WI( J ).EQ.ZERO ) THEN 00375 TNRM = DNRM2( N, VR( 1, J ), 1 ) 00376 ELSE IF( WI( J ).GT.ZERO ) THEN 00377 TNRM = DLAPY2( DNRM2( N, VR( 1, J ), 1 ), 00378 $ DNRM2( N, VR( 1, J+1 ), 1 ) ) 00379 END IF 00380 RESULT( 3 ) = MAX( RESULT( 3 ), 00381 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 00382 IF( WI( J ).GT.ZERO ) THEN 00383 VMX = ZERO 00384 VRMX = ZERO 00385 DO 20 JJ = 1, N 00386 VTST = DLAPY2( VR( JJ, J ), VR( JJ, J+1 ) ) 00387 IF( VTST.GT.VMX ) 00388 $ VMX = VTST 00389 IF( VR( JJ, J+1 ).EQ.ZERO .AND. ABS( VR( JJ, J ) ).GT. 00390 $ VRMX )VRMX = ABS( VR( JJ, J ) ) 00391 20 CONTINUE 00392 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 00393 $ RESULT( 3 ) = ULPINV 00394 END IF 00395 30 CONTINUE 00396 * 00397 * Do Test (4) 00398 * 00399 DO 50 J = 1, N 00400 TNRM = ONE 00401 IF( WI( J ).EQ.ZERO ) THEN 00402 TNRM = DNRM2( N, VL( 1, J ), 1 ) 00403 ELSE IF( WI( J ).GT.ZERO ) THEN 00404 TNRM = DLAPY2( DNRM2( N, VL( 1, J ), 1 ), 00405 $ DNRM2( N, VL( 1, J+1 ), 1 ) ) 00406 END IF 00407 RESULT( 4 ) = MAX( RESULT( 4 ), 00408 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 00409 IF( WI( J ).GT.ZERO ) THEN 00410 VMX = ZERO 00411 VRMX = ZERO 00412 DO 40 JJ = 1, N 00413 VTST = DLAPY2( VL( JJ, J ), VL( JJ, J+1 ) ) 00414 IF( VTST.GT.VMX ) 00415 $ VMX = VTST 00416 IF( VL( JJ, J+1 ).EQ.ZERO .AND. ABS( VL( JJ, J ) ).GT. 00417 $ VRMX )VRMX = ABS( VL( JJ, J ) ) 00418 40 CONTINUE 00419 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 00420 $ RESULT( 4 ) = ULPINV 00421 END IF 00422 50 CONTINUE 00423 * 00424 * Test for all options of computing condition numbers 00425 * 00426 DO 200 ISENS = 1, ISENSM 00427 * 00428 SENSE = SENS( ISENS ) 00429 * 00430 * Compute eigenvalues only, and test them 00431 * 00432 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 00433 CALL DGEEVX( BALANC, 'N', 'N', SENSE, N, H, LDA, WR1, WI1, DUM, 00434 $ 1, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1, 00435 $ RCNDV1, WORK, LWORK, IWORK, IINFO ) 00436 IF( IINFO.NE.0 ) THEN 00437 RESULT( 1 ) = ULPINV 00438 IF( JTYPE.NE.22 ) THEN 00439 WRITE( NOUNIT, FMT = 9998 )'DGEEVX2', IINFO, N, JTYPE, 00440 $ BALANC, ISEED 00441 ELSE 00442 WRITE( NOUNIT, FMT = 9999 )'DGEEVX2', IINFO, N, 00443 $ ISEED( 1 ) 00444 END IF 00445 INFO = ABS( IINFO ) 00446 GO TO 190 00447 END IF 00448 * 00449 * Do Test (5) 00450 * 00451 DO 60 J = 1, N 00452 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 00453 $ RESULT( 5 ) = ULPINV 00454 60 CONTINUE 00455 * 00456 * Do Test (8) 00457 * 00458 IF( .NOT.NOBAL ) THEN 00459 DO 70 J = 1, N 00460 IF( SCALE( J ).NE.SCALE1( J ) ) 00461 $ RESULT( 8 ) = ULPINV 00462 70 CONTINUE 00463 IF( ILO.NE.ILO1 ) 00464 $ RESULT( 8 ) = ULPINV 00465 IF( IHI.NE.IHI1 ) 00466 $ RESULT( 8 ) = ULPINV 00467 IF( ABNRM.NE.ABNRM1 ) 00468 $ RESULT( 8 ) = ULPINV 00469 END IF 00470 * 00471 * Do Test (9) 00472 * 00473 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN 00474 DO 80 J = 1, N 00475 IF( RCONDV( J ).NE.RCNDV1( J ) ) 00476 $ RESULT( 9 ) = ULPINV 00477 80 CONTINUE 00478 END IF 00479 * 00480 * Compute eigenvalues and right eigenvectors, and test them 00481 * 00482 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 00483 CALL DGEEVX( BALANC, 'N', 'V', SENSE, N, H, LDA, WR1, WI1, DUM, 00484 $ 1, LRE, LDLRE, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1, 00485 $ RCNDV1, WORK, LWORK, IWORK, IINFO ) 00486 IF( IINFO.NE.0 ) THEN 00487 RESULT( 1 ) = ULPINV 00488 IF( JTYPE.NE.22 ) THEN 00489 WRITE( NOUNIT, FMT = 9998 )'DGEEVX3', IINFO, N, JTYPE, 00490 $ BALANC, ISEED 00491 ELSE 00492 WRITE( NOUNIT, FMT = 9999 )'DGEEVX3', IINFO, N, 00493 $ ISEED( 1 ) 00494 END IF 00495 INFO = ABS( IINFO ) 00496 GO TO 190 00497 END IF 00498 * 00499 * Do Test (5) again 00500 * 00501 DO 90 J = 1, N 00502 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 00503 $ RESULT( 5 ) = ULPINV 00504 90 CONTINUE 00505 * 00506 * Do Test (6) 00507 * 00508 DO 110 J = 1, N 00509 DO 100 JJ = 1, N 00510 IF( VR( J, JJ ).NE.LRE( J, JJ ) ) 00511 $ RESULT( 6 ) = ULPINV 00512 100 CONTINUE 00513 110 CONTINUE 00514 * 00515 * Do Test (8) again 00516 * 00517 IF( .NOT.NOBAL ) THEN 00518 DO 120 J = 1, N 00519 IF( SCALE( J ).NE.SCALE1( J ) ) 00520 $ RESULT( 8 ) = ULPINV 00521 120 CONTINUE 00522 IF( ILO.NE.ILO1 ) 00523 $ RESULT( 8 ) = ULPINV 00524 IF( IHI.NE.IHI1 ) 00525 $ RESULT( 8 ) = ULPINV 00526 IF( ABNRM.NE.ABNRM1 ) 00527 $ RESULT( 8 ) = ULPINV 00528 END IF 00529 * 00530 * Do Test (9) again 00531 * 00532 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN 00533 DO 130 J = 1, N 00534 IF( RCONDV( J ).NE.RCNDV1( J ) ) 00535 $ RESULT( 9 ) = ULPINV 00536 130 CONTINUE 00537 END IF 00538 * 00539 * Compute eigenvalues and left eigenvectors, and test them 00540 * 00541 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 00542 CALL DGEEVX( BALANC, 'V', 'N', SENSE, N, H, LDA, WR1, WI1, LRE, 00543 $ LDLRE, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1, 00544 $ RCNDV1, WORK, LWORK, IWORK, IINFO ) 00545 IF( IINFO.NE.0 ) THEN 00546 RESULT( 1 ) = ULPINV 00547 IF( JTYPE.NE.22 ) THEN 00548 WRITE( NOUNIT, FMT = 9998 )'DGEEVX4', IINFO, N, JTYPE, 00549 $ BALANC, ISEED 00550 ELSE 00551 WRITE( NOUNIT, FMT = 9999 )'DGEEVX4', IINFO, N, 00552 $ ISEED( 1 ) 00553 END IF 00554 INFO = ABS( IINFO ) 00555 GO TO 190 00556 END IF 00557 * 00558 * Do Test (5) again 00559 * 00560 DO 140 J = 1, N 00561 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 00562 $ RESULT( 5 ) = ULPINV 00563 140 CONTINUE 00564 * 00565 * Do Test (7) 00566 * 00567 DO 160 J = 1, N 00568 DO 150 JJ = 1, N 00569 IF( VL( J, JJ ).NE.LRE( J, JJ ) ) 00570 $ RESULT( 7 ) = ULPINV 00571 150 CONTINUE 00572 160 CONTINUE 00573 * 00574 * Do Test (8) again 00575 * 00576 IF( .NOT.NOBAL ) THEN 00577 DO 170 J = 1, N 00578 IF( SCALE( J ).NE.SCALE1( J ) ) 00579 $ RESULT( 8 ) = ULPINV 00580 170 CONTINUE 00581 IF( ILO.NE.ILO1 ) 00582 $ RESULT( 8 ) = ULPINV 00583 IF( IHI.NE.IHI1 ) 00584 $ RESULT( 8 ) = ULPINV 00585 IF( ABNRM.NE.ABNRM1 ) 00586 $ RESULT( 8 ) = ULPINV 00587 END IF 00588 * 00589 * Do Test (9) again 00590 * 00591 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN 00592 DO 180 J = 1, N 00593 IF( RCONDV( J ).NE.RCNDV1( J ) ) 00594 $ RESULT( 9 ) = ULPINV 00595 180 CONTINUE 00596 END IF 00597 * 00598 190 CONTINUE 00599 * 00600 200 CONTINUE 00601 * 00602 * If COMP, compare condition numbers to precomputed ones 00603 * 00604 IF( COMP ) THEN 00605 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 00606 CALL DGEEVX( 'N', 'V', 'V', 'B', N, H, LDA, WR, WI, VL, LDVL, 00607 $ VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, 00608 $ WORK, LWORK, IWORK, IINFO ) 00609 IF( IINFO.NE.0 ) THEN 00610 RESULT( 1 ) = ULPINV 00611 WRITE( NOUNIT, FMT = 9999 )'DGEEVX5', IINFO, N, ISEED( 1 ) 00612 INFO = ABS( IINFO ) 00613 GO TO 250 00614 END IF 00615 * 00616 * Sort eigenvalues and condition numbers lexicographically 00617 * to compare with inputs 00618 * 00619 DO 220 I = 1, N - 1 00620 KMIN = I 00621 VRMIN = WR( I ) 00622 VIMIN = WI( I ) 00623 DO 210 J = I + 1, N 00624 IF( WR( J ).LT.VRMIN ) THEN 00625 KMIN = J 00626 VRMIN = WR( J ) 00627 VIMIN = WI( J ) 00628 END IF 00629 210 CONTINUE 00630 WR( KMIN ) = WR( I ) 00631 WI( KMIN ) = WI( I ) 00632 WR( I ) = VRMIN 00633 WI( I ) = VIMIN 00634 VRMIN = RCONDE( KMIN ) 00635 RCONDE( KMIN ) = RCONDE( I ) 00636 RCONDE( I ) = VRMIN 00637 VRMIN = RCONDV( KMIN ) 00638 RCONDV( KMIN ) = RCONDV( I ) 00639 RCONDV( I ) = VRMIN 00640 220 CONTINUE 00641 * 00642 * Compare condition numbers for eigenvectors 00643 * taking their condition numbers into account 00644 * 00645 RESULT( 10 ) = ZERO 00646 EPS = MAX( EPSIN, ULP ) 00647 V = MAX( DBLE( N )*EPS*ABNRM, SMLNUM ) 00648 IF( ABNRM.EQ.ZERO ) 00649 $ V = ONE 00650 DO 230 I = 1, N 00651 IF( V.GT.RCONDV( I )*RCONDE( I ) ) THEN 00652 TOL = RCONDV( I ) 00653 ELSE 00654 TOL = V / RCONDE( I ) 00655 END IF 00656 IF( V.GT.RCDVIN( I )*RCDEIN( I ) ) THEN 00657 TOLIN = RCDVIN( I ) 00658 ELSE 00659 TOLIN = V / RCDEIN( I ) 00660 END IF 00661 TOL = MAX( TOL, SMLNUM / EPS ) 00662 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00663 IF( EPS*( RCDVIN( I )-TOLIN ).GT.RCONDV( I )+TOL ) THEN 00664 VMAX = ONE / EPS 00665 ELSE IF( RCDVIN( I )-TOLIN.GT.RCONDV( I )+TOL ) THEN 00666 VMAX = ( RCDVIN( I )-TOLIN ) / ( RCONDV( I )+TOL ) 00667 ELSE IF( RCDVIN( I )+TOLIN.LT.EPS*( RCONDV( I )-TOL ) ) THEN 00668 VMAX = ONE / EPS 00669 ELSE IF( RCDVIN( I )+TOLIN.LT.RCONDV( I )-TOL ) THEN 00670 VMAX = ( RCONDV( I )-TOL ) / ( RCDVIN( I )+TOLIN ) 00671 ELSE 00672 VMAX = ONE 00673 END IF 00674 RESULT( 10 ) = MAX( RESULT( 10 ), VMAX ) 00675 230 CONTINUE 00676 * 00677 * Compare condition numbers for eigenvalues 00678 * taking their condition numbers into account 00679 * 00680 RESULT( 11 ) = ZERO 00681 DO 240 I = 1, N 00682 IF( V.GT.RCONDV( I ) ) THEN 00683 TOL = ONE 00684 ELSE 00685 TOL = V / RCONDV( I ) 00686 END IF 00687 IF( V.GT.RCDVIN( I ) ) THEN 00688 TOLIN = ONE 00689 ELSE 00690 TOLIN = V / RCDVIN( I ) 00691 END IF 00692 TOL = MAX( TOL, SMLNUM / EPS ) 00693 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00694 IF( EPS*( RCDEIN( I )-TOLIN ).GT.RCONDE( I )+TOL ) THEN 00695 VMAX = ONE / EPS 00696 ELSE IF( RCDEIN( I )-TOLIN.GT.RCONDE( I )+TOL ) THEN 00697 VMAX = ( RCDEIN( I )-TOLIN ) / ( RCONDE( I )+TOL ) 00698 ELSE IF( RCDEIN( I )+TOLIN.LT.EPS*( RCONDE( I )-TOL ) ) THEN 00699 VMAX = ONE / EPS 00700 ELSE IF( RCDEIN( I )+TOLIN.LT.RCONDE( I )-TOL ) THEN 00701 VMAX = ( RCONDE( I )-TOL ) / ( RCDEIN( I )+TOLIN ) 00702 ELSE 00703 VMAX = ONE 00704 END IF 00705 RESULT( 11 ) = MAX( RESULT( 11 ), VMAX ) 00706 240 CONTINUE 00707 250 CONTINUE 00708 * 00709 END IF 00710 * 00711 9999 FORMAT( ' DGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00712 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 ) 00713 9998 FORMAT( ' DGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00714 $ I6, ', JTYPE=', I6, ', BALANC = ', A, ', ISEED=(', 00715 $ 3( I5, ',' ), I5, ')' ) 00716 * 00717 RETURN 00718 * 00719 * End of DGET23 00720 * 00721 END
1.7.3 
