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