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