 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CDRVES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ NOUNIT, A, LDA, H, HT, W, WT, VS, LDVS, RESULT, 00003 $ WORK, NWORK, RWORK, 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 RESULT( 13 ), RWORK( * ) 00017 COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ), 00018 $ VS( LDVS, * ), W( * ), WORK( * ), WT( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CDRVES checks the nonsymmetric eigenvalue (Schur form) problem 00025 * driver CGEES. 00026 * 00027 * When CDRVES 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 W 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 W 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 complex angles. 00092 * (ULP = (first number larger than 1) - 1 ) 00093 * (5) A diagonal matrix with geometrically spaced entries 00094 * 1, ..., ULP and random complex angles. 00095 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00096 * and random complex angles. 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 unitary and 00104 * T has evenly spaced entries 1, ..., ULP with random 00105 * complex angles on the diagonal and random O(1) entries in 00106 * the upper triangle. 00107 * 00108 * (10) A matrix of the form U' T U, where U is unitary and 00109 * T has geometrically spaced entries 1, ..., ULP with random 00110 * complex angles on the diagonal and random O(1) entries in 00111 * the upper 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 * complex angles on the diagonal and random O(1) entries in 00116 * the upper triangle. 00117 * 00118 * (12) A matrix of the form U' T U, where U is unitary and 00119 * T has complex eigenvalues randomly chosen from 00120 * ULP < |z| < 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 complex angles on the diagonal and random O(1) 00126 * entries 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 complex angles on the diagonal 00131 * and random 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 complex angles on the diagonal and random O(1) 00136 * entries 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 complex eigenvalues randomly chosen 00140 * from ULP < |z| < 1 and random O(1) entries in the upper 00141 * 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 * CDRVES 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, CDRVES 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 CDRVES 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) COMPLEX 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) COMPLEX array, dimension (LDA, max(NN)) 00217 * Another copy of the test matrix A, modified by CGEES. 00218 * 00219 * HT (workspace) COMPLEX array, dimension (LDA, max(NN)) 00220 * Yet another copy of the test matrix A, modified by CGEES. 00221 * 00222 * W (workspace) COMPLEX array, dimension (max(NN)) 00223 * The computed eigenvalues of A. 00224 * 00225 * WT (workspace) COMPLEX array, dimension (max(NN)) 00226 * Like W, this array contains the eigenvalues of A, 00227 * but those computed when CGEES only computes a partial 00228 * eigendecomposition, i.e. not Schur vectors 00229 * 00230 * VS (workspace) COMPLEX array, dimension (LDVS, max(NN)) 00231 * VS holds the computed Schur vectors. 00232 * 00233 * LDVS (input) INTEGER 00234 * Leading dimension of VS. Must be at least max(1,max(NN)). 00235 * 00236 * RESULT (output) REAL array, dimension (13) 00237 * The values computed by the 13 tests described above. 00238 * The values are currently limited to 1/ulp, to avoid overflow. 00239 * 00240 * WORK (workspace) COMPLEX array, dimension (NWORK) 00241 * 00242 * NWORK (input) INTEGER 00243 * The number of entries in WORK. This must be at least 00244 * 5*NN(j)+2*NN(j)**2 for all j. 00245 * 00246 * RWORK (workspace) REAL array, dimension (max(NN)) 00247 * 00248 * IWORK (workspace) INTEGER array, dimension (max(NN)) 00249 * 00250 * INFO (output) INTEGER 00251 * If 0, then everything ran OK. 00252 * -1: NSIZES < 0 00253 * -2: Some NN(j) < 0 00254 * -3: NTYPES < 0 00255 * -6: THRESH < 0 00256 * -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). 00257 * -15: LDVS < 1 or LDVS < NMAX, where NMAX is max( NN(j) ). 00258 * -18: NWORK too small. 00259 * If CLATMR, CLATMS, CLATME or CGEES returns an error code, 00260 * the absolute value of it is returned. 00261 * 00262 *----------------------------------------------------------------------- 00263 * 00264 * Some Local Variables and Parameters: 00265 * ---- ----- --------- --- ---------- 00266 * ZERO, ONE Real 0 and 1. 00267 * MAXTYP The number of types defined. 00268 * NMAX Largest value in NN. 00269 * NERRS The number of tests which have exceeded THRESH 00270 * COND, CONDS, 00271 * IMODE Values to be passed to the matrix generators. 00272 * ANORM Norm of A; passed to matrix generators. 00273 * 00274 * OVFL, UNFL Overflow and underflow thresholds. 00275 * ULP, ULPINV Finest relative precision and its inverse. 00276 * RTULP, RTULPI Square roots of the previous 4 values. 00277 * The following four arrays decode JTYPE: 00278 * KTYPE(j) The general type (1-10) for type "j". 00279 * KMODE(j) The MODE value to be passed to the matrix 00280 * generator for type "j". 00281 * KMAGN(j) The order of magnitude ( O(1), 00282 * O(overflow^(1/2) ), O(underflow^(1/2) ) 00283 * KCONDS(j) Select whether CONDS is to be 1 or 00284 * 1/sqrt(ulp). (0 means irrelevant.) 00285 * 00286 * ===================================================================== 00287 * 00288 * .. Parameters .. 00289 COMPLEX CZERO 00290 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00291 COMPLEX CONE 00292 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00293 REAL ZERO, ONE 00294 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00295 INTEGER MAXTYP 00296 PARAMETER ( MAXTYP = 21 ) 00297 * .. 00298 * .. Local Scalars .. 00299 LOGICAL BADNN 00300 CHARACTER SORT 00301 CHARACTER*3 PATH 00302 INTEGER I, IINFO, IMODE, ISORT, ITYPE, IWK, J, JCOL, 00303 $ JSIZE, JTYPE, KNTEIG, LWORK, MTYPES, N, 00304 $ NERRS, NFAIL, NMAX, NNWORK, NTEST, NTESTF, 00305 $ NTESTT, RSUB, SDIM 00306 REAL ANORM, COND, CONDS, OVFL, RTULP, RTULPI, ULP, 00307 $ ULPINV, UNFL 00308 * .. 00309 * .. Local Arrays .. 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 CSLECT 00327 REAL SLAMCH 00328 EXTERNAL CSLECT, SLAMCH 00329 * .. 00330 * .. External Subroutines .. 00331 EXTERNAL CGEES, CHST01, CLACPY, CLATME, CLATMR, CLATMS, 00332 $ CLASET, SLABAD, SLASUM, XERBLA 00333 * .. 00334 * .. Intrinsic Functions .. 00335 INTRINSIC ABS, CMPLX, MAX, MIN, 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 ) = 'Complex 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 = -15 00383 ELSE IF( 5*NMAX+2*NMAX**2.GT.NWORK ) THEN 00384 INFO = -18 00385 END IF 00386 * 00387 IF( INFO.NE.0 ) THEN 00388 CALL XERBLA( 'CDRVES', -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 240 JSIZE = 1, NSIZES 00412 N = NN( JSIZE ) 00413 IF( NSIZES.NE.1 ) THEN 00414 MTYPES = MIN( MAXTYP, NTYPES ) 00415 ELSE 00416 MTYPES = MIN( MAXTYP+1, NTYPES ) 00417 END IF 00418 * 00419 DO 230 JTYPE = 1, MTYPES 00420 IF( .NOT.DOTYPE( JTYPE ) ) 00421 $ GO TO 230 00422 * 00423 * Save ISEED in case of an error. 00424 * 00425 DO 20 J = 1, 4 00426 IOLDSD( J ) = ISEED( J ) 00427 20 CONTINUE 00428 * 00429 * Compute "A" 00430 * 00431 * Control parameters: 00432 * 00433 * KMAGN KCONDS KMODE KTYPE 00434 * =1 O(1) 1 clustered 1 zero 00435 * =2 large large clustered 2 identity 00436 * =3 small exponential Jordan 00437 * =4 arithmetic diagonal, (w/ eigenvalues) 00438 * =5 random log symmetric, w/ eigenvalues 00439 * =6 random general, w/ eigenvalues 00440 * =7 random diagonal 00441 * =8 random symmetric 00442 * =9 random general 00443 * =10 random triangular 00444 * 00445 IF( MTYPES.GT.MAXTYP ) 00446 $ GO TO 90 00447 * 00448 ITYPE = KTYPE( JTYPE ) 00449 IMODE = KMODE( JTYPE ) 00450 * 00451 * Compute norm 00452 * 00453 GO TO ( 30, 40, 50 )KMAGN( JTYPE ) 00454 * 00455 30 CONTINUE 00456 ANORM = ONE 00457 GO TO 60 00458 * 00459 40 CONTINUE 00460 ANORM = OVFL*ULP 00461 GO TO 60 00462 * 00463 50 CONTINUE 00464 ANORM = UNFL*ULPINV 00465 GO TO 60 00466 * 00467 60 CONTINUE 00468 * 00469 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA ) 00470 IINFO = 0 00471 COND = ULPINV 00472 * 00473 * Special Matrices -- Identity & Jordan block 00474 * 00475 IF( ITYPE.EQ.1 ) THEN 00476 * 00477 * Zero 00478 * 00479 IINFO = 0 00480 * 00481 ELSE IF( ITYPE.EQ.2 ) THEN 00482 * 00483 * Identity 00484 * 00485 DO 70 JCOL = 1, N 00486 A( JCOL, JCOL ) = CMPLX( ANORM ) 00487 70 CONTINUE 00488 * 00489 ELSE IF( ITYPE.EQ.3 ) THEN 00490 * 00491 * Jordan Block 00492 * 00493 DO 80 JCOL = 1, N 00494 A( JCOL, JCOL ) = CMPLX( ANORM ) 00495 IF( JCOL.GT.1 ) 00496 $ A( JCOL, JCOL-1 ) = CONE 00497 80 CONTINUE 00498 * 00499 ELSE IF( ITYPE.EQ.4 ) THEN 00500 * 00501 * Diagonal Matrix, [Eigen]values Specified 00502 * 00503 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, 00504 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 00505 $ IINFO ) 00506 * 00507 ELSE IF( ITYPE.EQ.5 ) THEN 00508 * 00509 * Symmetric, eigenvalues specified 00510 * 00511 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, 00512 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 00513 $ IINFO ) 00514 * 00515 ELSE IF( ITYPE.EQ.6 ) THEN 00516 * 00517 * General, eigenvalues specified 00518 * 00519 IF( KCONDS( JTYPE ).EQ.1 ) THEN 00520 CONDS = ONE 00521 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 00522 CONDS = RTULPI 00523 ELSE 00524 CONDS = ZERO 00525 END IF 00526 * 00527 CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ', 00528 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM, 00529 $ A, LDA, WORK( 2*N+1 ), IINFO ) 00530 * 00531 ELSE IF( ITYPE.EQ.7 ) THEN 00532 * 00533 * Diagonal, random eigenvalues 00534 * 00535 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00536 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00537 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00538 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00539 * 00540 ELSE IF( ITYPE.EQ.8 ) THEN 00541 * 00542 * Symmetric, random eigenvalues 00543 * 00544 CALL CLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE, 00545 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00546 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00547 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00548 * 00549 ELSE IF( ITYPE.EQ.9 ) THEN 00550 * 00551 * General, random eigenvalues 00552 * 00553 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00554 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00555 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00556 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00557 IF( N.GE.4 ) THEN 00558 CALL CLASET( 'Full', 2, N, CZERO, CZERO, A, LDA ) 00559 CALL CLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ), 00560 $ LDA ) 00561 CALL CLASET( 'Full', N-3, 2, CZERO, CZERO, 00562 $ A( 3, N-1 ), LDA ) 00563 CALL CLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ), 00564 $ LDA ) 00565 END IF 00566 * 00567 ELSE IF( ITYPE.EQ.10 ) THEN 00568 * 00569 * Triangular, random eigenvalues 00570 * 00571 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00572 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00573 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 00574 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00575 * 00576 ELSE 00577 * 00578 IINFO = 1 00579 END IF 00580 * 00581 IF( IINFO.NE.0 ) THEN 00582 WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE, 00583 $ IOLDSD 00584 INFO = ABS( IINFO ) 00585 RETURN 00586 END IF 00587 * 00588 90 CONTINUE 00589 * 00590 * Test for minimal and generous workspace 00591 * 00592 DO 220 IWK = 1, 2 00593 IF( IWK.EQ.1 ) THEN 00594 NNWORK = 3*N 00595 ELSE 00596 NNWORK = 5*N + 2*N**2 00597 END IF 00598 NNWORK = MAX( NNWORK, 1 ) 00599 * 00600 * Initialize RESULT 00601 * 00602 DO 100 J = 1, 13 00603 RESULT( J ) = -ONE 00604 100 CONTINUE 00605 * 00606 * Test with and without sorting of eigenvalues 00607 * 00608 DO 180 ISORT = 0, 1 00609 IF( ISORT.EQ.0 ) THEN 00610 SORT = 'N' 00611 RSUB = 0 00612 ELSE 00613 SORT = 'S' 00614 RSUB = 6 00615 END IF 00616 * 00617 * Compute Schur form and Schur vectors, and test them 00618 * 00619 CALL CLACPY( 'F', N, N, A, LDA, H, LDA ) 00620 CALL CGEES( 'V', SORT, CSLECT, N, H, LDA, SDIM, W, VS, 00621 $ LDVS, WORK, NNWORK, RWORK, BWORK, IINFO ) 00622 IF( IINFO.NE.0 ) THEN 00623 RESULT( 1+RSUB ) = ULPINV 00624 WRITE( NOUNIT, FMT = 9992 )'CGEES1', IINFO, N, 00625 $ JTYPE, IOLDSD 00626 INFO = ABS( IINFO ) 00627 GO TO 190 00628 END IF 00629 * 00630 * Do Test (1) or Test (7) 00631 * 00632 RESULT( 1+RSUB ) = ZERO 00633 DO 120 J = 1, N - 1 00634 DO 110 I = J + 1, N 00635 IF( H( I, J ).NE.ZERO ) 00636 $ RESULT( 1+RSUB ) = ULPINV 00637 110 CONTINUE 00638 120 CONTINUE 00639 * 00640 * Do Tests (2) and (3) or Tests (8) and (9) 00641 * 00642 LWORK = MAX( 1, 2*N*N ) 00643 CALL CHST01( N, 1, N, A, LDA, H, LDA, VS, LDVS, WORK, 00644 $ LWORK, RWORK, RES ) 00645 RESULT( 2+RSUB ) = RES( 1 ) 00646 RESULT( 3+RSUB ) = RES( 2 ) 00647 * 00648 * Do Test (4) or Test (10) 00649 * 00650 RESULT( 4+RSUB ) = ZERO 00651 DO 130 I = 1, N 00652 IF( H( I, I ).NE.W( I ) ) 00653 $ RESULT( 4+RSUB ) = ULPINV 00654 130 CONTINUE 00655 * 00656 * Do Test (5) or Test (11) 00657 * 00658 CALL CLACPY( 'F', N, N, A, LDA, HT, LDA ) 00659 CALL CGEES( 'N', SORT, CSLECT, N, HT, LDA, SDIM, WT, 00660 $ VS, LDVS, WORK, NNWORK, RWORK, BWORK, 00661 $ IINFO ) 00662 IF( IINFO.NE.0 ) THEN 00663 RESULT( 5+RSUB ) = ULPINV 00664 WRITE( NOUNIT, FMT = 9992 )'CGEES2', IINFO, N, 00665 $ JTYPE, IOLDSD 00666 INFO = ABS( IINFO ) 00667 GO TO 190 00668 END IF 00669 * 00670 RESULT( 5+RSUB ) = ZERO 00671 DO 150 J = 1, N 00672 DO 140 I = 1, N 00673 IF( H( I, J ).NE.HT( I, J ) ) 00674 $ RESULT( 5+RSUB ) = ULPINV 00675 140 CONTINUE 00676 150 CONTINUE 00677 * 00678 * Do Test (6) or Test (12) 00679 * 00680 RESULT( 6+RSUB ) = ZERO 00681 DO 160 I = 1, N 00682 IF( W( I ).NE.WT( I ) ) 00683 $ RESULT( 6+RSUB ) = ULPINV 00684 160 CONTINUE 00685 * 00686 * Do Test (13) 00687 * 00688 IF( ISORT.EQ.1 ) THEN 00689 RESULT( 13 ) = ZERO 00690 KNTEIG = 0 00691 DO 170 I = 1, N 00692 IF( CSLECT( W( I ) ) ) 00693 $ KNTEIG = KNTEIG + 1 00694 IF( I.LT.N ) THEN 00695 IF( CSLECT( W( I+1 ) ) .AND. 00696 $ ( .NOT.CSLECT( W( I ) ) ) )RESULT( 13 ) 00697 $ = ULPINV 00698 END IF 00699 170 CONTINUE 00700 IF( SDIM.NE.KNTEIG ) 00701 $ RESULT( 13 ) = ULPINV 00702 END IF 00703 * 00704 180 CONTINUE 00705 * 00706 * End of Loop -- Check for RESULT(j) > THRESH 00707 * 00708 190 CONTINUE 00709 * 00710 NTEST = 0 00711 NFAIL = 0 00712 DO 200 J = 1, 13 00713 IF( RESULT( J ).GE.ZERO ) 00714 $ NTEST = NTEST + 1 00715 IF( RESULT( J ).GE.THRESH ) 00716 $ NFAIL = NFAIL + 1 00717 200 CONTINUE 00718 * 00719 IF( NFAIL.GT.0 ) 00720 $ NTESTF = NTESTF + 1 00721 IF( NTESTF.EQ.1 ) THEN 00722 WRITE( NOUNIT, FMT = 9999 )PATH 00723 WRITE( NOUNIT, FMT = 9998 ) 00724 WRITE( NOUNIT, FMT = 9997 ) 00725 WRITE( NOUNIT, FMT = 9996 ) 00726 WRITE( NOUNIT, FMT = 9995 )THRESH 00727 WRITE( NOUNIT, FMT = 9994 ) 00728 NTESTF = 2 00729 END IF 00730 * 00731 DO 210 J = 1, 13 00732 IF( RESULT( J ).GE.THRESH ) THEN 00733 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE, 00734 $ J, RESULT( J ) 00735 END IF 00736 210 CONTINUE 00737 * 00738 NERRS = NERRS + NFAIL 00739 NTESTT = NTESTT + NTEST 00740 * 00741 220 CONTINUE 00742 230 CONTINUE 00743 240 CONTINUE 00744 * 00745 * Summary 00746 * 00747 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT ) 00748 * 00749 9999 FORMAT( / 1X, A3, ' -- Complex Schur Form Decomposition Driver', 00750 $ / ' Matrix types (see CDRVES for details): ' ) 00751 * 00752 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ', 00753 $ ' ', ' 5=Diagonal: geometr. spaced entries.', 00754 $ / ' 2=Identity matrix. ', ' 6=Diagona', 00755 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ', 00756 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ', 00757 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s', 00758 $ 'mall, evenly spaced.' ) 00759 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev', 00760 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e', 00761 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ', 00762 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond', 00763 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp', 00764 $ 'lex ', A6, / ' 12=Well-cond., random complex ', A6, ' ', 00765 $ ' 17=Ill-cond., large rand. complx ', A4, / ' 13=Ill-condi', 00766 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.', 00767 $ ' complx ', A4 ) 00768 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ', 00769 $ 'with small random entries.', / ' 20=Matrix with large ran', 00770 $ 'dom entries. ', / ) 00771 9995 FORMAT( ' Tests performed with test threshold =', F8.2, 00772 $ / ' ( A denotes A on input and T denotes A on output)', 00773 $ / / ' 1 = 0 if T in Schur form (no sort), ', 00774 $ ' 1/ulp otherwise', / 00775 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)', 00776 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', 00777 $ / ' 4 = 0 if W are eigenvalues of T (no sort),', 00778 $ ' 1/ulp otherwise', / 00779 $ ' 5 = 0 if T same no matter if VS computed (no sort),', 00780 $ ' 1/ulp otherwise', / 00781 $ ' 6 = 0 if W same no matter if VS computed (no sort)', 00782 $ ', 1/ulp otherwise' ) 00783 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise', 00784 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)', 00785 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ', 00786 $ / ' 10 = 0 if W are eigenvalues of T (sort),', 00787 $ ' 1/ulp otherwise', / 00788 $ ' 11 = 0 if T same no matter if VS computed (sort),', 00789 $ ' 1/ulp otherwise', / 00790 $ ' 12 = 0 if W same no matter if VS computed (sort),', 00791 $ ' 1/ulp otherwise', / 00792 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise', / ) 00793 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ), 00794 $ ' type ', I2, ', test(', I2, ')=', G10.3 ) 00795 9992 FORMAT( ' CDRVES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00796 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00797 * 00798 RETURN 00799 * 00800 * End of CDRVES 00801 * 00802 END
1.7.3 
