 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT, 00003 $ WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK, 00004 $ LWORK, IWORK, BWORK, 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 DOUBLE PRECISION THRESH 00014 * .. 00015 * .. Array Arguments .. 00016 LOGICAL BWORK( * ), DOTYPE( * ) 00017 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00018 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ), 00019 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ), 00020 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ), 00021 $ WR( * ), WRT( * ), WRTMP( * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * DDRVSX checks the nonsymmetric eigenvalue (Schur form) problem 00028 * expert driver DGEESX. 00029 * 00030 * DDRVSX 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 DDRVSX 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 WR+sqrt(-1)*WI 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 WR+sqrt(-1)*WI are eigenvalues of T 00074 * 1/ulp otherwise 00075 * If workspace sufficient, also compare WR, WI 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 signs. 00112 * (ULP = (first number larger than 1) - 1 ) 00113 * (5) A diagonal matrix with geometrically spaced entries 00114 * 1, ..., ULP and random signs. 00115 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00116 * and random signs. 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 orthogonal and 00124 * T has evenly spaced entries 1, ..., ULP with random signs 00125 * on the diagonal and random O(1) entries in the upper 00126 * triangle. 00127 * 00128 * (10) A matrix of the form U' T U, where U is orthogonal and 00129 * T has geometrically spaced entries 1, ..., ULP with random 00130 * signs on the diagonal and random O(1) entries in the upper 00131 * 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 * signs on the diagonal and random O(1) entries in the upper 00136 * triangle. 00137 * 00138 * (12) A matrix of the form U' T U, where U is orthogonal and 00139 * T has real or complex conjugate paired eigenvalues randomly 00140 * chosen from ( ULP, 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 signs on the diagonal and random O(1) entries 00146 * 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 signs on the diagonal and random 00151 * 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 signs on the diagonal and random O(1) entries 00156 * 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 real or complex conjugate paired 00160 * eigenvalues randomly chosen from ( ULP, 1 ) and random 00161 * O(1) entries in the upper 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 DGEESX 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 DGEESX 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 DDRVSX to continue the same random number 00242 * sequence. 00243 * 00244 * THRESH (input) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDA, max(NN)) 00269 * Another copy of the test matrix A, modified by DGEESX. 00270 * 00271 * HT (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN)) 00272 * Yet another copy of the test matrix A, modified by DGEESX. 00273 * 00274 * WR (workspace) DOUBLE PRECISION array, dimension (max(NN)) 00275 * WI (workspace) DOUBLE PRECISION array, dimension (max(NN)) 00276 * The real and imaginary parts of the eigenvalues of A. 00277 * On exit, WR + WI*i are the eigenvalues of the matrix in A. 00278 * 00279 * WRT (workspace) DOUBLE PRECISION array, dimension (max(NN)) 00280 * WIT (workspace) DOUBLE PRECISION array, dimension (max(NN)) 00281 * Like WR, WI, these arrays contain the eigenvalues of A, 00282 * but those computed when DGEESX only computes a partial 00283 * eigendecomposition, i.e. not Schur vectors 00284 * 00285 * WRTMP (workspace) DOUBLE PRECISION array, dimension (max(NN)) 00286 * WITMP (workspace) DOUBLE PRECISION array, dimension (max(NN)) 00287 * More temporary storage for eigenvalues. 00288 * 00289 * VS (workspace) DOUBLE PRECISION array, dimension (LDVS, max(NN)) 00290 * VS holds the computed Schur vectors. 00291 * 00292 * LDVS (input) INTEGER 00293 * Leading dimension of VS. Must be at least max(1,max(NN)). 00294 * 00295 * VS1 (workspace) DOUBLE PRECISION array, dimension (LDVS, max(NN)) 00296 * VS1 holds another copy of the computed Schur vectors. 00297 * 00298 * RESULT (output) DOUBLE PRECISION array, dimension (17) 00299 * The values computed by the 17 tests described above. 00300 * The values are currently limited to 1/ulp, to avoid overflow. 00301 * 00302 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00303 * 00304 * LWORK (input) INTEGER 00305 * The number of entries in WORK. This must be at least 00306 * max(3*NN(j),2*NN(j)**2) for all j. 00307 * 00308 * IWORK (workspace) INTEGER array, dimension (max(NN)*max(NN)) 00309 * 00310 * INFO (output) INTEGER 00311 * If 0, successful exit. 00312 * <0, input parameter -INFO is incorrect 00313 * >0, DLATMR, SLATMS, SLATME or DGET24 returned an error 00314 * code and INFO is its absolute value 00315 * 00316 *----------------------------------------------------------------------- 00317 * 00318 * Some Local Variables and Parameters: 00319 * ---- ----- --------- --- ---------- 00320 * ZERO, ONE Real 0 and 1. 00321 * MAXTYP The number of types defined. 00322 * NMAX Largest value in NN. 00323 * NERRS The number of tests which have exceeded THRESH 00324 * COND, CONDS, 00325 * IMODE Values to be passed to the matrix generators. 00326 * ANORM Norm of A; passed to matrix generators. 00327 * 00328 * OVFL, UNFL Overflow and underflow thresholds. 00329 * ULP, ULPINV Finest relative precision and its inverse. 00330 * RTULP, RTULPI Square roots of the previous 4 values. 00331 * The following four arrays decode JTYPE: 00332 * KTYPE(j) The general type (1-10) for type "j". 00333 * KMODE(j) The MODE value to be passed to the matrix 00334 * generator for type "j". 00335 * KMAGN(j) The order of magnitude ( O(1), 00336 * O(overflow^(1/2) ), O(underflow^(1/2) ) 00337 * KCONDS(j) Selectw whether CONDS is to be 1 or 00338 * 1/sqrt(ulp). (0 means irrelevant.) 00339 * 00340 * ===================================================================== 00341 * 00342 * .. Parameters .. 00343 DOUBLE PRECISION ZERO, ONE 00344 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00345 INTEGER MAXTYP 00346 PARAMETER ( MAXTYP = 21 ) 00347 * .. 00348 * .. Local Scalars .. 00349 LOGICAL BADNN 00350 CHARACTER*3 PATH 00351 INTEGER I, IINFO, IMODE, ITYPE, IWK, J, JCOL, JSIZE, 00352 $ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX, NNWORK, 00353 $ NSLCT, NTEST, NTESTF, NTESTT 00354 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN, 00355 $ RTULP, RTULPI, ULP, ULPINV, UNFL 00356 * .. 00357 * .. Local Arrays .. 00358 CHARACTER ADUMMA( 1 ) 00359 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ), 00360 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ), 00361 $ KMODE( MAXTYP ), KTYPE( MAXTYP ) 00362 * .. 00363 * .. Arrays in Common .. 00364 LOGICAL SELVAL( 20 ) 00365 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 ) 00366 * .. 00367 * .. Scalars in Common .. 00368 INTEGER SELDIM, SELOPT 00369 * .. 00370 * .. Common blocks .. 00371 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI 00372 * .. 00373 * .. External Functions .. 00374 DOUBLE PRECISION DLAMCH 00375 EXTERNAL DLAMCH 00376 * .. 00377 * .. External Subroutines .. 00378 EXTERNAL DGET24, DLABAD, DLASET, DLASUM, DLATME, DLATMR, 00379 $ DLATMS, XERBLA 00380 * .. 00381 * .. Intrinsic Functions .. 00382 INTRINSIC ABS, MAX, MIN, SQRT 00383 * .. 00384 * .. Data statements .. 00385 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 00386 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 00387 $ 3, 1, 2, 3 / 00388 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 00389 $ 1, 5, 5, 5, 4, 3, 1 / 00390 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 00391 * .. 00392 * .. Executable Statements .. 00393 * 00394 PATH( 1: 1 ) = 'Double precision' 00395 PATH( 2: 3 ) = 'SX' 00396 * 00397 * Check for errors 00398 * 00399 NTESTT = 0 00400 NTESTF = 0 00401 INFO = 0 00402 * 00403 * Important constants 00404 * 00405 BADNN = .FALSE. 00406 * 00407 * 12 is the largest dimension in the input file of precomputed 00408 * problems 00409 * 00410 NMAX = 12 00411 DO 10 J = 1, NSIZES 00412 NMAX = MAX( NMAX, NN( J ) ) 00413 IF( NN( J ).LT.0 ) 00414 $ BADNN = .TRUE. 00415 10 CONTINUE 00416 * 00417 * Check for errors 00418 * 00419 IF( NSIZES.LT.0 ) THEN 00420 INFO = -1 00421 ELSE IF( BADNN ) THEN 00422 INFO = -2 00423 ELSE IF( NTYPES.LT.0 ) THEN 00424 INFO = -3 00425 ELSE IF( THRESH.LT.ZERO ) THEN 00426 INFO = -6 00427 ELSE IF( NIUNIT.LE.0 ) THEN 00428 INFO = -7 00429 ELSE IF( NOUNIT.LE.0 ) THEN 00430 INFO = -8 00431 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN 00432 INFO = -10 00433 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN 00434 INFO = -20 00435 ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN 00436 INFO = -24 00437 END IF 00438 * 00439 IF( INFO.NE.0 ) THEN 00440 CALL XERBLA( 'DDRVSX', -INFO ) 00441 RETURN 00442 END IF 00443 * 00444 * If nothing to do check on NIUNIT 00445 * 00446 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00447 $ GO TO 150 00448 * 00449 * More Important constants 00450 * 00451 UNFL = DLAMCH( 'Safe minimum' ) 00452 OVFL = ONE / UNFL 00453 CALL DLABAD( UNFL, OVFL ) 00454 ULP = DLAMCH( 'Precision' ) 00455 ULPINV = ONE / ULP 00456 RTULP = SQRT( ULP ) 00457 RTULPI = ONE / RTULP 00458 * 00459 * Loop over sizes, types 00460 * 00461 NERRS = 0 00462 * 00463 DO 140 JSIZE = 1, NSIZES 00464 N = NN( JSIZE ) 00465 IF( NSIZES.NE.1 ) THEN 00466 MTYPES = MIN( MAXTYP, NTYPES ) 00467 ELSE 00468 MTYPES = MIN( MAXTYP+1, NTYPES ) 00469 END IF 00470 * 00471 DO 130 JTYPE = 1, MTYPES 00472 IF( .NOT.DOTYPE( JTYPE ) ) 00473 $ GO TO 130 00474 * 00475 * Save ISEED in case of an error. 00476 * 00477 DO 20 J = 1, 4 00478 IOLDSD( J ) = ISEED( J ) 00479 20 CONTINUE 00480 * 00481 * Compute "A" 00482 * 00483 * Control parameters: 00484 * 00485 * KMAGN KCONDS KMODE KTYPE 00486 * =1 O(1) 1 clustered 1 zero 00487 * =2 large large clustered 2 identity 00488 * =3 small exponential Jordan 00489 * =4 arithmetic diagonal, (w/ eigenvalues) 00490 * =5 random log symmetric, w/ eigenvalues 00491 * =6 random general, w/ eigenvalues 00492 * =7 random diagonal 00493 * =8 random symmetric 00494 * =9 random general 00495 * =10 random triangular 00496 * 00497 IF( MTYPES.GT.MAXTYP ) 00498 $ GO TO 90 00499 * 00500 ITYPE = KTYPE( JTYPE ) 00501 IMODE = KMODE( JTYPE ) 00502 * 00503 * Compute norm 00504 * 00505 GO TO ( 30, 40, 50 )KMAGN( JTYPE ) 00506 * 00507 30 CONTINUE 00508 ANORM = ONE 00509 GO TO 60 00510 * 00511 40 CONTINUE 00512 ANORM = OVFL*ULP 00513 GO TO 60 00514 * 00515 50 CONTINUE 00516 ANORM = UNFL*ULPINV 00517 GO TO 60 00518 * 00519 60 CONTINUE 00520 * 00521 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 00522 IINFO = 0 00523 COND = ULPINV 00524 * 00525 * Special Matrices -- Identity & Jordan block 00526 * 00527 * Zero 00528 * 00529 IF( ITYPE.EQ.1 ) THEN 00530 IINFO = 0 00531 * 00532 ELSE IF( ITYPE.EQ.2 ) THEN 00533 * 00534 * Identity 00535 * 00536 DO 70 JCOL = 1, N 00537 A( JCOL, JCOL ) = ANORM 00538 70 CONTINUE 00539 * 00540 ELSE IF( ITYPE.EQ.3 ) THEN 00541 * 00542 * Jordan Block 00543 * 00544 DO 80 JCOL = 1, N 00545 A( JCOL, JCOL ) = ANORM 00546 IF( JCOL.GT.1 ) 00547 $ A( JCOL, JCOL-1 ) = ONE 00548 80 CONTINUE 00549 * 00550 ELSE IF( ITYPE.EQ.4 ) THEN 00551 * 00552 * Diagonal Matrix, [Eigen]values Specified 00553 * 00554 CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00555 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 00556 $ IINFO ) 00557 * 00558 ELSE IF( ITYPE.EQ.5 ) THEN 00559 * 00560 * Symmetric, eigenvalues specified 00561 * 00562 CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00563 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 00564 $ IINFO ) 00565 * 00566 ELSE IF( ITYPE.EQ.6 ) THEN 00567 * 00568 * General, eigenvalues specified 00569 * 00570 IF( KCONDS( JTYPE ).EQ.1 ) THEN 00571 CONDS = ONE 00572 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 00573 CONDS = RTULPI 00574 ELSE 00575 CONDS = ZERO 00576 END IF 00577 * 00578 ADUMMA( 1 ) = ' ' 00579 CALL DLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 00580 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 00581 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 00582 $ IINFO ) 00583 * 00584 ELSE IF( ITYPE.EQ.7 ) THEN 00585 * 00586 * Diagonal, random eigenvalues 00587 * 00588 CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 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, IWORK, IINFO ) 00592 * 00593 ELSE IF( ITYPE.EQ.8 ) THEN 00594 * 00595 * Symmetric, random eigenvalues 00596 * 00597 CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 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, IWORK, IINFO ) 00601 * 00602 ELSE IF( ITYPE.EQ.9 ) THEN 00603 * 00604 * General, random eigenvalues 00605 * 00606 CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 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, IWORK, IINFO ) 00610 IF( N.GE.4 ) THEN 00611 CALL DLASET( 'Full', 2, N, ZERO, ZERO, A, LDA ) 00612 CALL DLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ), 00613 $ LDA ) 00614 CALL DLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ), 00615 $ LDA ) 00616 CALL DLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ), 00617 $ LDA ) 00618 END IF 00619 * 00620 ELSE IF( ITYPE.EQ.10 ) THEN 00621 * 00622 * Triangular, random eigenvalues 00623 * 00624 CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 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, IWORK, 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 = 3*N 00648 ELSE 00649 NNWORK = MAX( 3*N, 2*N*N ) 00650 END IF 00651 NNWORK = MAX( NNWORK, 1 ) 00652 * 00653 CALL DGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N, 00654 $ A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP, 00655 $ WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, 00656 $ ISLCT, RESULT, WORK, NNWORK, IWORK, BWORK, 00657 $ INFO ) 00658 * 00659 * Check for RESULT(j) > THRESH 00660 * 00661 NTEST = 0 00662 NFAIL = 0 00663 DO 100 J = 1, 15 00664 IF( RESULT( J ).GE.ZERO ) 00665 $ NTEST = NTEST + 1 00666 IF( RESULT( J ).GE.THRESH ) 00667 $ NFAIL = NFAIL + 1 00668 100 CONTINUE 00669 * 00670 IF( NFAIL.GT.0 ) 00671 $ NTESTF = NTESTF + 1 00672 IF( NTESTF.EQ.1 ) THEN 00673 WRITE( NOUNIT, FMT = 9999 )PATH 00674 WRITE( NOUNIT, FMT = 9998 ) 00675 WRITE( NOUNIT, FMT = 9997 ) 00676 WRITE( NOUNIT, FMT = 9996 ) 00677 WRITE( NOUNIT, FMT = 9995 )THRESH 00678 WRITE( NOUNIT, FMT = 9994 ) 00679 NTESTF = 2 00680 END IF 00681 * 00682 DO 110 J = 1, 15 00683 IF( RESULT( J ).GE.THRESH ) THEN 00684 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE, 00685 $ J, RESULT( J ) 00686 END IF 00687 110 CONTINUE 00688 * 00689 NERRS = NERRS + NFAIL 00690 NTESTT = NTESTT + NTEST 00691 * 00692 120 CONTINUE 00693 130 CONTINUE 00694 140 CONTINUE 00695 * 00696 150 CONTINUE 00697 * 00698 * Read in data from file to check accuracy of condition estimation 00699 * Read input data until N=0 00700 * 00701 JTYPE = 0 00702 160 CONTINUE 00703 READ( NIUNIT, FMT = *, END = 200 )N, NSLCT 00704 IF( N.EQ.0 ) 00705 $ GO TO 200 00706 JTYPE = JTYPE + 1 00707 ISEED( 1 ) = JTYPE 00708 IF( NSLCT.GT.0 ) 00709 $ READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT ) 00710 DO 170 I = 1, N 00711 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N ) 00712 170 CONTINUE 00713 READ( NIUNIT, FMT = * )RCDEIN, RCDVIN 00714 * 00715 CALL DGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, 00716 $ WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1, 00717 $ RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK, 00718 $ IWORK, BWORK, INFO ) 00719 * 00720 * Check for RESULT(j) > THRESH 00721 * 00722 NTEST = 0 00723 NFAIL = 0 00724 DO 180 J = 1, 17 00725 IF( RESULT( J ).GE.ZERO ) 00726 $ NTEST = NTEST + 1 00727 IF( RESULT( J ).GE.THRESH ) 00728 $ NFAIL = NFAIL + 1 00729 180 CONTINUE 00730 * 00731 IF( NFAIL.GT.0 ) 00732 $ NTESTF = NTESTF + 1 00733 IF( NTESTF.EQ.1 ) THEN 00734 WRITE( NOUNIT, FMT = 9999 )PATH 00735 WRITE( NOUNIT, FMT = 9998 ) 00736 WRITE( NOUNIT, FMT = 9997 ) 00737 WRITE( NOUNIT, FMT = 9996 ) 00738 WRITE( NOUNIT, FMT = 9995 )THRESH 00739 WRITE( NOUNIT, FMT = 9994 ) 00740 NTESTF = 2 00741 END IF 00742 DO 190 J = 1, 17 00743 IF( RESULT( J ).GE.THRESH ) THEN 00744 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J ) 00745 END IF 00746 190 CONTINUE 00747 * 00748 NERRS = NERRS + NFAIL 00749 NTESTT = NTESTT + NTEST 00750 GO TO 160 00751 200 CONTINUE 00752 * 00753 * Summary 00754 * 00755 CALL DLASUM( PATH, NOUNIT, NERRS, NTESTT ) 00756 * 00757 9999 FORMAT( / 1X, A3, ' -- Real Schur Form Decomposition Expert ', 00758 $ 'Driver', / ' Matrix types (see DDRVSX for details):' ) 00759 * 00760 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ', 00761 $ ' ', ' 5=Diagonal: geometr. spaced entries.', 00762 $ / ' 2=Identity matrix. ', ' 6=Diagona', 00763 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ', 00764 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ', 00765 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s', 00766 $ 'mall, evenly spaced.' ) 00767 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev', 00768 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e', 00769 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ', 00770 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond', 00771 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp', 00772 $ 'lex ', / ' 12=Well-cond., random complex ', ' ', 00773 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi', 00774 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.', 00775 $ ' complx ' ) 00776 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ', 00777 $ 'with small random entries.', / ' 20=Matrix with large ran', 00778 $ 'dom entries. ', / ) 00779 9995 FORMAT( ' Tests performed with test threshold =', F8.2, 00780 $ / ' ( A denotes A on input and T denotes A on output)', 00781 $ / / ' 1 = 0 if T in Schur form (no sort), ', 00782 $ ' 1/ulp otherwise', / 00783 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)', 00784 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', / 00785 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),', 00786 $ ' 1/ulp otherwise', / 00787 $ ' 5 = 0 if T same no matter if VS computed (no sort),', 00788 $ ' 1/ulp otherwise', / 00789 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)', 00790 $ ', 1/ulp otherwise' ) 00791 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise', 00792 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)', 00793 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ', 00794 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),', 00795 $ ' 1/ulp otherwise', / 00796 $ ' 11 = 0 if T same no matter what else computed (sort),', 00797 $ ' 1/ulp otherwise', / 00798 $ ' 12 = 0 if WR, WI same no matter what else computed ', 00799 $ '(sort), 1/ulp otherwise', / 00800 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise', 00801 $ / ' 14 = 0 if RCONDE same no matter what else computed,', 00802 $ ' 1/ulp otherwise', / 00803 $ ' 15 = 0 if RCONDv same no matter what else computed,', 00804 $ ' 1/ulp otherwise', / 00805 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),', 00806 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' ) 00807 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ), 00808 $ ' type ', I2, ', test(', I2, ')=', G10.3 ) 00809 9992 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=', 00810 $ G10.3 ) 00811 9991 FORMAT( ' DDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00812 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00813 * 00814 RETURN 00815 * 00816 * End of DDRVSX 00817 * 00818 END
1.7.3 
