 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CDRVVX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ NIUNIT, NOUNIT, A, LDA, H, W, W1, VL, LDVL, VR, 00003 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN, 00004 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT, 00005 $ WORK, NWORK, RWORK, 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 INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NIUNIT, NOUNIT, 00013 $ NSIZES, NTYPES, NWORK 00014 REAL THRESH 00015 * .. 00016 * .. Array Arguments .. 00017 LOGICAL DOTYPE( * ) 00018 INTEGER ISEED( 4 ), NN( * ) 00019 REAL RCDEIN( * ), RCDVIN( * ), RCNDE1( * ), 00020 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ), 00021 $ RESULT( 11 ), RWORK( * ), SCALE( * ), 00022 $ SCALE1( * ) 00023 COMPLEX A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ), 00024 $ VL( LDVL, * ), VR( LDVR, * ), W( * ), W1( * ), 00025 $ WORK( * ) 00026 * .. 00027 * 00028 * Purpose 00029 * ======= 00030 * 00031 * CDRVVX checks the nonsymmetric eigenvalue problem expert driver 00032 * CGEEVX. 00033 * 00034 * CDRVVX uses both test matrices generated randomly depending on 00035 * data supplied in the calling sequence, as well as on data 00036 * read from an input file and including precomputed condition 00037 * numbers to which it compares the ones it computes. 00038 * 00039 * When CDRVVX is called, a number of matrix "sizes" ("n's") and a 00040 * number of matrix "types" are specified in the calling sequence. 00041 * For each size ("n") and each type of matrix, one matrix will be 00042 * generated and used to test the nonsymmetric eigenroutines. For 00043 * each matrix, 9 tests will be performed: 00044 * 00045 * (1) | A * VR - VR * W | / ( n |A| ulp ) 00046 * 00047 * Here VR is the matrix of unit right eigenvectors. 00048 * W is a diagonal matrix with diagonal entries W(j). 00049 * 00050 * (2) | A**H * VL - VL * W**H | / ( n |A| ulp ) 00051 * 00052 * Here VL is the matrix of unit left eigenvectors, A**H is the 00053 * conjugate transpose of A, and W is as above. 00054 * 00055 * (3) | |VR(i)| - 1 | / ulp and largest component real 00056 * 00057 * VR(i) denotes the i-th column of VR. 00058 * 00059 * (4) | |VL(i)| - 1 | / ulp and largest component real 00060 * 00061 * VL(i) denotes the i-th column of VL. 00062 * 00063 * (5) W(full) = W(partial) 00064 * 00065 * W(full) denotes the eigenvalues computed when VR, VL, RCONDV 00066 * and RCONDE are also computed, and W(partial) denotes the 00067 * eigenvalues computed when only some of VR, VL, RCONDV, and 00068 * RCONDE are computed. 00069 * 00070 * (6) VR(full) = VR(partial) 00071 * 00072 * VR(full) denotes the right eigenvectors computed when VL, RCONDV 00073 * and RCONDE are computed, and VR(partial) denotes the result 00074 * when only some of VL and RCONDV are computed. 00075 * 00076 * (7) VL(full) = VL(partial) 00077 * 00078 * VL(full) denotes the left eigenvectors computed when VR, RCONDV 00079 * and RCONDE are computed, and VL(partial) denotes the result 00080 * when only some of VR and RCONDV are computed. 00081 * 00082 * (8) 0 if SCALE, ILO, IHI, ABNRM (full) = 00083 * SCALE, ILO, IHI, ABNRM (partial) 00084 * 1/ulp otherwise 00085 * 00086 * SCALE, ILO, IHI and ABNRM describe how the matrix is balanced. 00087 * (full) is when VR, VL, RCONDE and RCONDV are also computed, and 00088 * (partial) is when some are not computed. 00089 * 00090 * (9) RCONDV(full) = RCONDV(partial) 00091 * 00092 * RCONDV(full) denotes the reciprocal condition numbers of the 00093 * right eigenvectors computed when VR, VL and RCONDE are also 00094 * computed. RCONDV(partial) denotes the reciprocal condition 00095 * numbers when only some of VR, VL and RCONDE are computed. 00096 * 00097 * The "sizes" are specified by an array NN(1:NSIZES); the value of 00098 * each element NN(j) specifies one size. 00099 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 00100 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00101 * Currently, the list of possible types is: 00102 * 00103 * (1) The zero matrix. 00104 * (2) The identity matrix. 00105 * (3) A (transposed) Jordan block, with 1's on the diagonal. 00106 * 00107 * (4) A diagonal matrix with evenly spaced entries 00108 * 1, ..., ULP and random complex angles. 00109 * (ULP = (first number larger than 1) - 1 ) 00110 * (5) A diagonal matrix with geometrically spaced entries 00111 * 1, ..., ULP and random complex angles. 00112 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00113 * and random complex angles. 00114 * 00115 * (7) Same as (4), but multiplied by a constant near 00116 * the overflow threshold 00117 * (8) Same as (4), but multiplied by a constant near 00118 * the underflow threshold 00119 * 00120 * (9) A matrix of the form U' T U, where U is unitary and 00121 * T has evenly spaced entries 1, ..., ULP with random complex 00122 * angles on the diagonal and random O(1) entries in the upper 00123 * triangle. 00124 * 00125 * (10) A matrix of the form U' T U, where U is unitary and 00126 * T has geometrically spaced entries 1, ..., ULP with random 00127 * complex angles on the diagonal and random O(1) entries in 00128 * the upper triangle. 00129 * 00130 * (11) A matrix of the form U' T U, where U is unitary and 00131 * T has "clustered" entries 1, ULP,..., ULP with random 00132 * complex angles on the diagonal and random O(1) entries in 00133 * the upper triangle. 00134 * 00135 * (12) A matrix of the form U' T U, where U is unitary and 00136 * T has complex eigenvalues randomly chosen from 00137 * ULP < |z| < 1 and random O(1) entries in the upper 00138 * triangle. 00139 * 00140 * (13) A matrix of the form X' T X, where X has condition 00141 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP 00142 * with random complex angles on the diagonal and random O(1) 00143 * entries in the upper triangle. 00144 * 00145 * (14) A matrix of the form X' T X, where X has condition 00146 * SQRT( ULP ) and T has geometrically spaced entries 00147 * 1, ..., ULP with random complex angles on the diagonal 00148 * and random O(1) entries in the upper triangle. 00149 * 00150 * (15) A matrix of the form X' T X, where X has condition 00151 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP 00152 * with random complex angles on the diagonal and random O(1) 00153 * entries in the upper triangle. 00154 * 00155 * (16) A matrix of the form X' T X, where X has condition 00156 * SQRT( ULP ) and T has complex eigenvalues randomly chosen 00157 * from ULP < |z| < 1 and random O(1) entries in the upper 00158 * triangle. 00159 * 00160 * (17) Same as (16), but multiplied by a constant 00161 * near the overflow threshold 00162 * (18) Same as (16), but multiplied by a constant 00163 * near the underflow threshold 00164 * 00165 * (19) Nonsymmetric matrix with random entries chosen from |z| < 1 00166 * If N is at least 4, all entries in first two rows and last 00167 * row, and first column and last two columns are zero. 00168 * (20) Same as (19), but multiplied by a constant 00169 * near the overflow threshold 00170 * (21) Same as (19), but multiplied by a constant 00171 * near the underflow threshold 00172 * 00173 * In addition, an input file will be read from logical unit number 00174 * NIUNIT. The file contains matrices along with precomputed 00175 * eigenvalues and reciprocal condition numbers for the eigenvalues 00176 * and right eigenvectors. For these matrices, in addition to tests 00177 * (1) to (9) we will compute the following two tests: 00178 * 00179 * (10) |RCONDV - RCDVIN| / cond(RCONDV) 00180 * 00181 * RCONDV is the reciprocal right eigenvector condition number 00182 * computed by CGEEVX and RCDVIN (the precomputed true value) 00183 * is supplied as input. cond(RCONDV) is the condition number of 00184 * RCONDV, and takes errors in computing RCONDV into account, so 00185 * that the resulting quantity should be O(ULP). cond(RCONDV) is 00186 * essentially given by norm(A)/RCONDE. 00187 * 00188 * (11) |RCONDE - RCDEIN| / cond(RCONDE) 00189 * 00190 * RCONDE is the reciprocal eigenvalue condition number 00191 * computed by CGEEVX and RCDEIN (the precomputed true value) 00192 * is supplied as input. cond(RCONDE) is the condition number 00193 * of RCONDE, and takes errors in computing RCONDE into account, 00194 * so that the resulting quantity should be O(ULP). cond(RCONDE) 00195 * is essentially given by norm(A)/RCONDV. 00196 * 00197 * Arguments 00198 * ========== 00199 * 00200 * NSIZES (input) INTEGER 00201 * The number of sizes of matrices to use. NSIZES must be at 00202 * least zero. If it is zero, no randomly generated matrices 00203 * are tested, but any test matrices read from NIUNIT will be 00204 * tested. 00205 * 00206 * NN (input) INTEGER array, dimension (NSIZES) 00207 * An array containing the sizes to be used for the matrices. 00208 * Zero values will be skipped. The values must be at least 00209 * zero. 00210 * 00211 * NTYPES (input) INTEGER 00212 * The number of elements in DOTYPE. NTYPES must be at least 00213 * zero. If it is zero, no randomly generated test matrices 00214 * are tested, but and test matrices read from NIUNIT will be 00215 * tested. If it is MAXTYP+1 and NSIZES is 1, then an 00216 * additional type, MAXTYP+1 is defined, which is to use 00217 * whatever matrix is in A. This is only useful if 00218 * DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. . 00219 * 00220 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00221 * If DOTYPE(j) is .TRUE., then for each size in NN a 00222 * matrix of that size and of type j will be generated. 00223 * If NTYPES is smaller than the maximum number of types 00224 * defined (PARAMETER MAXTYP), then types NTYPES+1 through 00225 * MAXTYP will not be generated. If NTYPES is larger 00226 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 00227 * will be ignored. 00228 * 00229 * ISEED (input/output) INTEGER array, dimension (4) 00230 * On entry ISEED specifies the seed of the random number 00231 * generator. The array elements should be between 0 and 4095; 00232 * if not they will be reduced mod 4096. Also, ISEED(4) must 00233 * be odd. The random number generator uses a linear 00234 * congruential sequence limited to small integers, and so 00235 * should produce machine independent random numbers. The 00236 * values of ISEED are changed on exit, and can be used in the 00237 * next call to CDRVVX to continue the same random number 00238 * sequence. 00239 * 00240 * THRESH (input) REAL 00241 * A test will count as "failed" if the "error", computed as 00242 * described above, exceeds THRESH. Note that the error 00243 * is scaled to be O(1), so THRESH should be a reasonably 00244 * small multiple of 1, e.g., 10 or 100. In particular, 00245 * it should not depend on the precision (single vs. double) 00246 * or the size of the matrix. It must be at least zero. 00247 * 00248 * NIUNIT (input) INTEGER 00249 * The FORTRAN unit number for reading in the data file of 00250 * problems to solve. 00251 * 00252 * NOUNIT (input) INTEGER 00253 * The FORTRAN unit number for printing out error messages 00254 * (e.g., if a routine returns INFO not equal to 0.) 00255 * 00256 * A (workspace) COMPLEX array, dimension (LDA, max(NN,12)) 00257 * Used to hold the matrix whose eigenvalues are to be 00258 * computed. On exit, A contains the last matrix actually used. 00259 * 00260 * LDA (input) INTEGER 00261 * The leading dimension of A, and H. LDA must be at 00262 * least 1 and at least max( NN, 12 ). (12 is the 00263 * dimension of the largest matrix on the precomputed 00264 * input file.) 00265 * 00266 * H (workspace) COMPLEX array, dimension (LDA, max(NN,12)) 00267 * Another copy of the test matrix A, modified by CGEEVX. 00268 * 00269 * W (workspace) COMPLEX array, dimension (max(NN,12)) 00270 * Contains the eigenvalues of A. 00271 * 00272 * W1 (workspace) COMPLEX array, dimension (max(NN,12)) 00273 * Like W, this array contains the eigenvalues of A, 00274 * but those computed when CGEEVX only computes a partial 00275 * eigendecomposition, i.e. not the eigenvalues and left 00276 * and right eigenvectors. 00277 * 00278 * VL (workspace) COMPLEX array, dimension (LDVL, max(NN,12)) 00279 * VL holds the computed left eigenvectors. 00280 * 00281 * LDVL (input) INTEGER 00282 * Leading dimension of VL. Must be at least max(1,max(NN,12)). 00283 * 00284 * VR (workspace) COMPLEX array, dimension (LDVR, max(NN,12)) 00285 * VR holds the computed right eigenvectors. 00286 * 00287 * LDVR (input) INTEGER 00288 * Leading dimension of VR. Must be at least max(1,max(NN,12)). 00289 * 00290 * LRE (workspace) COMPLEX array, dimension (LDLRE, max(NN,12)) 00291 * LRE holds the computed right or left eigenvectors. 00292 * 00293 * LDLRE (input) INTEGER 00294 * Leading dimension of LRE. Must be at least max(1,max(NN,12)) 00295 * 00296 * RESULT (output) REAL array, dimension (11) 00297 * The values computed by the seven tests described above. 00298 * The values are currently limited to 1/ulp, to avoid 00299 * overflow. 00300 * 00301 * WORK (workspace) COMPLEX array, dimension (NWORK) 00302 * 00303 * NWORK (input) INTEGER 00304 * The number of entries in WORK. This must be at least 00305 * max(6*12+2*12**2,6*NN(j)+2*NN(j)**2) = 00306 * max( 360 ,6*NN(j)+2*NN(j)**2) for all j. 00307 * 00308 * RWORK (workspace) REAL array, dimension (2*max(NN,12)) 00309 * 00310 * INFO (output) INTEGER 00311 * If 0, then successful exit. 00312 * If <0, then input paramter -INFO is incorrect. 00313 * If >0, CLATMR, CLATMS, CLATME or CGET23 returned an error 00314 * code, and INFO is its absolute value. 00315 * 00316 *----------------------------------------------------------------------- 00317 * 00318 * Some Local Variables and Parameters: 00319 * ---- ----- --------- --- ---------- 00320 * 00321 * ZERO, ONE Real 0 and 1. 00322 * MAXTYP The number of types defined. 00323 * NMAX Largest value in NN or 12. 00324 * NERRS The number of tests which have exceeded THRESH 00325 * COND, CONDS, 00326 * IMODE Values to be passed to the matrix generators. 00327 * ANORM Norm of A; passed to matrix generators. 00328 * 00329 * OVFL, UNFL Overflow and underflow thresholds. 00330 * ULP, ULPINV Finest relative precision and its inverse. 00331 * RTULP, RTULPI Square roots of the previous 4 values. 00332 * 00333 * The following four arrays decode JTYPE: 00334 * KTYPE(j) The general type (1-10) for type "j". 00335 * KMODE(j) The MODE value to be passed to the matrix 00336 * generator for type "j". 00337 * KMAGN(j) The order of magnitude ( O(1), 00338 * O(overflow^(1/2) ), O(underflow^(1/2) ) 00339 * KCONDS(j) Selectw whether CONDS is to be 1 or 00340 * 1/sqrt(ulp). (0 means irrelevant.) 00341 * 00342 * ===================================================================== 00343 * 00344 * .. Parameters .. 00345 COMPLEX CZERO 00346 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00347 COMPLEX CONE 00348 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00349 REAL ZERO, ONE 00350 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00351 INTEGER MAXTYP 00352 PARAMETER ( MAXTYP = 21 ) 00353 * .. 00354 * .. Local Scalars .. 00355 LOGICAL BADNN 00356 CHARACTER BALANC 00357 CHARACTER*3 PATH 00358 INTEGER I, IBAL, IINFO, IMODE, ISRT, ITYPE, IWK, J, 00359 $ JCOL, JSIZE, JTYPE, MTYPES, N, NERRS, 00360 $ NFAIL, NMAX, NNWORK, NTEST, NTESTF, NTESTT 00361 REAL ANORM, COND, CONDS, OVFL, RTULP, RTULPI, ULP, 00362 $ ULPINV, UNFL, WI, WR 00363 * .. 00364 * .. Local Arrays .. 00365 CHARACTER BAL( 4 ) 00366 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), 00367 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 00368 $ KTYPE( MAXTYP ) 00369 * .. 00370 * .. External Functions .. 00371 REAL SLAMCH 00372 EXTERNAL SLAMCH 00373 * .. 00374 * .. External Subroutines .. 00375 EXTERNAL CGET23, CLATME, CLATMR, CLATMS, CLASET, SLABAD, 00376 $ SLASUM, XERBLA 00377 * .. 00378 * .. Intrinsic Functions .. 00379 INTRINSIC ABS, CMPLX, MAX, MIN, SQRT 00380 * .. 00381 * .. Data statements .. 00382 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 00383 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 00384 $ 3, 1, 2, 3 / 00385 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 00386 $ 1, 5, 5, 5, 4, 3, 1 / 00387 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 00388 DATA BAL / 'N', 'P', 'S', 'B' / 00389 * .. 00390 * .. Executable Statements .. 00391 * 00392 PATH( 1: 1 ) = 'Complex precision' 00393 PATH( 2: 3 ) = 'VX' 00394 * 00395 * Check for errors 00396 * 00397 NTESTT = 0 00398 NTESTF = 0 00399 INFO = 0 00400 * 00401 * Important constants 00402 * 00403 BADNN = .FALSE. 00404 * 00405 * 7 is the largest dimension in the input file of precomputed 00406 * problems 00407 * 00408 NMAX = 7 00409 DO 10 J = 1, NSIZES 00410 NMAX = MAX( NMAX, NN( J ) ) 00411 IF( NN( J ).LT.0 ) 00412 $ BADNN = .TRUE. 00413 10 CONTINUE 00414 * 00415 * Check for errors 00416 * 00417 IF( NSIZES.LT.0 ) THEN 00418 INFO = -1 00419 ELSE IF( BADNN ) THEN 00420 INFO = -2 00421 ELSE IF( NTYPES.LT.0 ) THEN 00422 INFO = -3 00423 ELSE IF( THRESH.LT.ZERO ) THEN 00424 INFO = -6 00425 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN 00426 INFO = -10 00427 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.NMAX ) THEN 00428 INFO = -15 00429 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.NMAX ) THEN 00430 INFO = -17 00431 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.NMAX ) THEN 00432 INFO = -19 00433 ELSE IF( 6*NMAX+2*NMAX**2.GT.NWORK ) THEN 00434 INFO = -30 00435 END IF 00436 * 00437 IF( INFO.NE.0 ) THEN 00438 CALL XERBLA( 'CDRVVX', -INFO ) 00439 RETURN 00440 END IF 00441 * 00442 * If nothing to do check on NIUNIT 00443 * 00444 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00445 $ GO TO 160 00446 * 00447 * More Important constants 00448 * 00449 UNFL = SLAMCH( 'Safe minimum' ) 00450 OVFL = ONE / UNFL 00451 CALL SLABAD( UNFL, OVFL ) 00452 ULP = SLAMCH( 'Precision' ) 00453 ULPINV = ONE / ULP 00454 RTULP = SQRT( ULP ) 00455 RTULPI = ONE / RTULP 00456 * 00457 * Loop over sizes, types 00458 * 00459 NERRS = 0 00460 * 00461 DO 150 JSIZE = 1, NSIZES 00462 N = NN( JSIZE ) 00463 IF( NSIZES.NE.1 ) THEN 00464 MTYPES = MIN( MAXTYP, NTYPES ) 00465 ELSE 00466 MTYPES = MIN( MAXTYP+1, NTYPES ) 00467 END IF 00468 * 00469 DO 140 JTYPE = 1, MTYPES 00470 IF( .NOT.DOTYPE( JTYPE ) ) 00471 $ GO TO 140 00472 * 00473 * Save ISEED in case of an error. 00474 * 00475 DO 20 J = 1, 4 00476 IOLDSD( J ) = ISEED( J ) 00477 20 CONTINUE 00478 * 00479 * Compute "A" 00480 * 00481 * Control parameters: 00482 * 00483 * KMAGN KCONDS KMODE KTYPE 00484 * =1 O(1) 1 clustered 1 zero 00485 * =2 large large clustered 2 identity 00486 * =3 small exponential Jordan 00487 * =4 arithmetic diagonal, (w/ eigenvalues) 00488 * =5 random log symmetric, w/ eigenvalues 00489 * =6 random general, w/ eigenvalues 00490 * =7 random diagonal 00491 * =8 random symmetric 00492 * =9 random general 00493 * =10 random triangular 00494 * 00495 IF( MTYPES.GT.MAXTYP ) 00496 $ GO TO 90 00497 * 00498 ITYPE = KTYPE( JTYPE ) 00499 IMODE = KMODE( JTYPE ) 00500 * 00501 * Compute norm 00502 * 00503 GO TO ( 30, 40, 50 )KMAGN( JTYPE ) 00504 * 00505 30 CONTINUE 00506 ANORM = ONE 00507 GO TO 60 00508 * 00509 40 CONTINUE 00510 ANORM = OVFL*ULP 00511 GO TO 60 00512 * 00513 50 CONTINUE 00514 ANORM = UNFL*ULPINV 00515 GO TO 60 00516 * 00517 60 CONTINUE 00518 * 00519 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA ) 00520 IINFO = 0 00521 COND = ULPINV 00522 * 00523 * Special Matrices -- Identity & Jordan block 00524 * 00525 * Zero 00526 * 00527 IF( ITYPE.EQ.1 ) THEN 00528 IINFO = 0 00529 * 00530 ELSE IF( ITYPE.EQ.2 ) THEN 00531 * 00532 * Identity 00533 * 00534 DO 70 JCOL = 1, N 00535 A( JCOL, JCOL ) = ANORM 00536 70 CONTINUE 00537 * 00538 ELSE IF( ITYPE.EQ.3 ) THEN 00539 * 00540 * Jordan Block 00541 * 00542 DO 80 JCOL = 1, N 00543 A( JCOL, JCOL ) = ANORM 00544 IF( JCOL.GT.1 ) 00545 $ A( JCOL, JCOL-1 ) = ONE 00546 80 CONTINUE 00547 * 00548 ELSE IF( ITYPE.EQ.4 ) THEN 00549 * 00550 * Diagonal Matrix, [Eigen]values Specified 00551 * 00552 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, 00553 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 00554 $ IINFO ) 00555 * 00556 ELSE IF( ITYPE.EQ.5 ) THEN 00557 * 00558 * Symmetric, eigenvalues specified 00559 * 00560 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, 00561 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 00562 $ IINFO ) 00563 * 00564 ELSE IF( ITYPE.EQ.6 ) THEN 00565 * 00566 * General, eigenvalues specified 00567 * 00568 IF( KCONDS( JTYPE ).EQ.1 ) THEN 00569 CONDS = ONE 00570 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 00571 CONDS = RTULPI 00572 ELSE 00573 CONDS = ZERO 00574 END IF 00575 * 00576 CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ', 00577 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM, 00578 $ A, LDA, WORK( 2*N+1 ), IINFO ) 00579 * 00580 ELSE IF( ITYPE.EQ.7 ) THEN 00581 * 00582 * Diagonal, random eigenvalues 00583 * 00584 CALL CLATMR( N, N, 'D', ISEED, 'S', WORK, 6, ONE, CONE, 00585 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00586 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00587 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO ) 00588 * 00589 ELSE IF( ITYPE.EQ.8 ) THEN 00590 * 00591 * Symmetric, random eigenvalues 00592 * 00593 CALL CLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE, 00594 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00595 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00596 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO ) 00597 * 00598 ELSE IF( ITYPE.EQ.9 ) THEN 00599 * 00600 * General, random eigenvalues 00601 * 00602 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00603 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00604 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00605 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO ) 00606 IF( N.GE.4 ) THEN 00607 CALL CLASET( 'Full', 2, N, CZERO, CZERO, A, LDA ) 00608 CALL CLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ), 00609 $ LDA ) 00610 CALL CLASET( 'Full', N-3, 2, CZERO, CZERO, 00611 $ A( 3, N-1 ), LDA ) 00612 CALL CLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ), 00613 $ LDA ) 00614 END IF 00615 * 00616 ELSE IF( ITYPE.EQ.10 ) THEN 00617 * 00618 * Triangular, random eigenvalues 00619 * 00620 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00621 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00622 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 00623 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO ) 00624 * 00625 ELSE 00626 * 00627 IINFO = 1 00628 END IF 00629 * 00630 IF( IINFO.NE.0 ) THEN 00631 WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE, 00632 $ IOLDSD 00633 INFO = ABS( IINFO ) 00634 RETURN 00635 END IF 00636 * 00637 90 CONTINUE 00638 * 00639 * Test for minimal and generous workspace 00640 * 00641 DO 130 IWK = 1, 3 00642 IF( IWK.EQ.1 ) THEN 00643 NNWORK = 2*N 00644 ELSE IF( IWK.EQ.2 ) THEN 00645 NNWORK = 2*N + N**2 00646 ELSE 00647 NNWORK = 6*N + 2*N**2 00648 END IF 00649 NNWORK = MAX( NNWORK, 1 ) 00650 * 00651 * Test for all balancing options 00652 * 00653 DO 120 IBAL = 1, 4 00654 BALANC = BAL( IBAL ) 00655 * 00656 * Perform tests 00657 * 00658 CALL CGET23( .FALSE., 0, BALANC, JTYPE, THRESH, 00659 $ IOLDSD, NOUNIT, N, A, LDA, H, W, W1, VL, 00660 $ LDVL, VR, LDVR, LRE, LDLRE, RCONDV, 00661 $ RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN, 00662 $ SCALE, SCALE1, RESULT, WORK, NNWORK, 00663 $ RWORK, INFO ) 00664 * 00665 * Check for RESULT(j) > THRESH 00666 * 00667 NTEST = 0 00668 NFAIL = 0 00669 DO 100 J = 1, 9 00670 IF( RESULT( J ).GE.ZERO ) 00671 $ NTEST = NTEST + 1 00672 IF( RESULT( J ).GE.THRESH ) 00673 $ NFAIL = NFAIL + 1 00674 100 CONTINUE 00675 * 00676 IF( NFAIL.GT.0 ) 00677 $ NTESTF = NTESTF + 1 00678 IF( NTESTF.EQ.1 ) THEN 00679 WRITE( NOUNIT, FMT = 9999 )PATH 00680 WRITE( NOUNIT, FMT = 9998 ) 00681 WRITE( NOUNIT, FMT = 9997 ) 00682 WRITE( NOUNIT, FMT = 9996 ) 00683 WRITE( NOUNIT, FMT = 9995 )THRESH 00684 NTESTF = 2 00685 END IF 00686 * 00687 DO 110 J = 1, 9 00688 IF( RESULT( J ).GE.THRESH ) THEN 00689 WRITE( NOUNIT, FMT = 9994 )BALANC, N, IWK, 00690 $ IOLDSD, JTYPE, J, RESULT( J ) 00691 END IF 00692 110 CONTINUE 00693 * 00694 NERRS = NERRS + NFAIL 00695 NTESTT = NTESTT + NTEST 00696 * 00697 120 CONTINUE 00698 130 CONTINUE 00699 140 CONTINUE 00700 150 CONTINUE 00701 * 00702 160 CONTINUE 00703 * 00704 * Read in data from file to check accuracy of condition estimation. 00705 * Assume input eigenvalues are sorted lexicographically (increasing 00706 * by real part, then decreasing by imaginary part) 00707 * 00708 JTYPE = 0 00709 170 CONTINUE 00710 READ( NIUNIT, FMT = *, END = 220 )N, ISRT 00711 * 00712 * Read input data until N=0 00713 * 00714 IF( N.EQ.0 ) 00715 $ GO TO 220 00716 JTYPE = JTYPE + 1 00717 ISEED( 1 ) = JTYPE 00718 DO 180 I = 1, N 00719 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N ) 00720 180 CONTINUE 00721 DO 190 I = 1, N 00722 READ( NIUNIT, FMT = * )WR, WI, RCDEIN( I ), RCDVIN( I ) 00723 W1( I ) = CMPLX( WR, WI ) 00724 190 CONTINUE 00725 CALL CGET23( .TRUE., ISRT, 'N', 22, THRESH, ISEED, NOUNIT, N, A, 00726 $ LDA, H, W, W1, VL, LDVL, VR, LDVR, LRE, LDLRE, 00727 $ RCONDV, RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN, 00728 $ SCALE, SCALE1, RESULT, WORK, 6*N+2*N**2, RWORK, 00729 $ INFO ) 00730 * 00731 * Check for RESULT(j) > THRESH 00732 * 00733 NTEST = 0 00734 NFAIL = 0 00735 DO 200 J = 1, 11 00736 IF( RESULT( J ).GE.ZERO ) 00737 $ NTEST = NTEST + 1 00738 IF( RESULT( J ).GE.THRESH ) 00739 $ NFAIL = NFAIL + 1 00740 200 CONTINUE 00741 * 00742 IF( NFAIL.GT.0 ) 00743 $ NTESTF = NTESTF + 1 00744 IF( NTESTF.EQ.1 ) THEN 00745 WRITE( NOUNIT, FMT = 9999 )PATH 00746 WRITE( NOUNIT, FMT = 9998 ) 00747 WRITE( NOUNIT, FMT = 9997 ) 00748 WRITE( NOUNIT, FMT = 9996 ) 00749 WRITE( NOUNIT, FMT = 9995 )THRESH 00750 NTESTF = 2 00751 END IF 00752 * 00753 DO 210 J = 1, 11 00754 IF( RESULT( J ).GE.THRESH ) THEN 00755 WRITE( NOUNIT, FMT = 9993 )N, JTYPE, J, RESULT( J ) 00756 END IF 00757 210 CONTINUE 00758 * 00759 NERRS = NERRS + NFAIL 00760 NTESTT = NTESTT + NTEST 00761 GO TO 170 00762 220 CONTINUE 00763 * 00764 * Summary 00765 * 00766 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT ) 00767 * 00768 9999 FORMAT( / 1X, A3, ' -- Complex Eigenvalue-Eigenvector ', 00769 $ 'Decomposition Expert Driver', 00770 $ / ' Matrix types (see CDRVVX for details): ' ) 00771 * 00772 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ', 00773 $ ' ', ' 5=Diagonal: geometr. spaced entries.', 00774 $ / ' 2=Identity matrix. ', ' 6=Diagona', 00775 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ', 00776 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ', 00777 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s', 00778 $ 'mall, evenly spaced.' ) 00779 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev', 00780 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e', 00781 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ', 00782 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond', 00783 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp', 00784 $ 'lex ', / ' 12=Well-cond., random complex ', ' ', 00785 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi', 00786 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.', 00787 $ ' complx ' ) 00788 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ', 00789 $ 'with small random entries.', / ' 20=Matrix with large ran', 00790 $ 'dom entries. ', ' 22=Matrix read from input file', / ) 00791 9995 FORMAT( ' Tests performed with test threshold =', F8.2, 00792 $ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ', 00793 $ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ', 00794 $ / ' 3 = | |VR(i)| - 1 | / ulp ', 00795 $ / ' 4 = | |VL(i)| - 1 | / ulp ', 00796 $ / ' 5 = 0 if W same no matter if VR or VL computed,', 00797 $ ' 1/ulp otherwise', / 00798 $ ' 6 = 0 if VR same no matter what else computed,', 00799 $ ' 1/ulp otherwise', / 00800 $ ' 7 = 0 if VL same no matter what else computed,', 00801 $ ' 1/ulp otherwise', / 00802 $ ' 8 = 0 if RCONDV same no matter what else computed,', 00803 $ ' 1/ulp otherwise', / 00804 $ ' 9 = 0 if SCALE, ILO, IHI, ABNRM same no matter what else', 00805 $ ' computed, 1/ulp otherwise', 00806 $ / ' 10 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),', 00807 $ / ' 11 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),' ) 00808 9994 FORMAT( ' BALANC=''', A1, ''',N=', I4, ',IWK=', I1, ', seed=', 00809 $ 4( I4, ',' ), ' type ', I2, ', test(', I2, ')=', G10.3 ) 00810 9993 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=', 00811 $ G10.3 ) 00812 9992 FORMAT( ' CDRVVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00813 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00814 * 00815 RETURN 00816 * 00817 * End of CDRVVX 00818 * 00819 END
1.7.3 
