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