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