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