 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, W1, 00003 $ W3, EVECTL, EVECTR, EVECTY, EVECTX, UU, TAU, 00004 $ WORK, NWORK, RWORK, IWORK, SELECT, RESULT, 00005 $ INFO ) 00006 * 00007 * -- LAPACK test routine (version 3.3.1) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * -- April 2011 -- 00010 * 00011 * .. Scalar Arguments .. 00012 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK 00013 REAL THRESH 00014 * .. 00015 * .. Array Arguments .. 00016 LOGICAL DOTYPE( * ), SELECT( * ) 00017 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00018 REAL RESULT( 14 ), RWORK( * ) 00019 COMPLEX A( LDA, * ), EVECTL( LDU, * ), 00020 $ EVECTR( LDU, * ), EVECTX( LDU, * ), 00021 $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ), 00022 $ T2( LDA, * ), TAU( * ), U( LDU, * ), 00023 $ UU( LDU, * ), UZ( LDU, * ), W1( * ), W3( * ), 00024 $ WORK( * ), Z( LDU, * ) 00025 * .. 00026 * 00027 * Purpose 00028 * ======= 00029 * 00030 * CCHKHS checks the nonsymmetric eigenvalue problem routines. 00031 * 00032 * CGEHRD factors A as U H U' , where ' means conjugate 00033 * transpose, H is hessenberg, and U is unitary. 00034 * 00035 * CUNGHR generates the unitary matrix U. 00036 * 00037 * CUNMHR multiplies a matrix by the unitary matrix U. 00038 * 00039 * CHSEQR factors H as Z T Z' , where Z is unitary and T 00040 * is upper triangular. It also computes the eigenvalues, 00041 * w(1), ..., w(n); we define a diagonal matrix W whose 00042 * (diagonal) entries are the eigenvalues. 00043 * 00044 * CTREVC computes the left eigenvector matrix L and the 00045 * right eigenvector matrix R for the matrix T. The 00046 * columns of L are the complex conjugates of the left 00047 * eigenvectors of T. The columns of R are the right 00048 * eigenvectors of T. L is lower triangular, and R is 00049 * upper triangular. 00050 * 00051 * CHSEIN computes the left eigenvector matrix Y and the 00052 * right eigenvector matrix X for the matrix H. The 00053 * columns of Y are the complex conjugates of the left 00054 * eigenvectors of H. The columns of X are the right 00055 * eigenvectors of H. Y is lower triangular, and X is 00056 * upper triangular. 00057 * 00058 * When CCHKHS is called, a number of matrix "sizes" ("n's") and a 00059 * number of matrix "types" are specified. For each size ("n") 00060 * and each type of matrix, one matrix will be generated and used 00061 * to test the nonsymmetric eigenroutines. For each matrix, 14 00062 * tests will be performed: 00063 * 00064 * (1) | A - U H U**H | / ( |A| n ulp ) 00065 * 00066 * (2) | I - UU**H | / ( n ulp ) 00067 * 00068 * (3) | H - Z T Z**H | / ( |H| n ulp ) 00069 * 00070 * (4) | I - ZZ**H | / ( n ulp ) 00071 * 00072 * (5) | A - UZ H (UZ)**H | / ( |A| n ulp ) 00073 * 00074 * (6) | I - UZ (UZ)**H | / ( n ulp ) 00075 * 00076 * (7) | T(Z computed) - T(Z not computed) | / ( |T| ulp ) 00077 * 00078 * (8) | W(Z computed) - W(Z not computed) | / ( |W| ulp ) 00079 * 00080 * (9) | TR - RW | / ( |T| |R| ulp ) 00081 * 00082 * (10) | L**H T - W**H L | / ( |T| |L| ulp ) 00083 * 00084 * (11) | HX - XW | / ( |H| |X| ulp ) 00085 * 00086 * (12) | Y**H H - W**H Y | / ( |H| |Y| ulp ) 00087 * 00088 * (13) | AX - XW | / ( |A| |X| ulp ) 00089 * 00090 * (14) | Y**H A - W**H Y | / ( |A| |Y| ulp ) 00091 * 00092 * The "sizes" are specified by an array NN(1:NSIZES); the value of 00093 * each element NN(j) specifies one size. 00094 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 00095 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00096 * Currently, the list of possible types is: 00097 * 00098 * (1) The zero matrix. 00099 * (2) The identity matrix. 00100 * (3) A (transposed) Jordan block, with 1's on the diagonal. 00101 * 00102 * (4) A diagonal matrix with evenly spaced entries 00103 * 1, ..., ULP and random complex angles. 00104 * (ULP = (first number larger than 1) - 1 ) 00105 * (5) A diagonal matrix with geometrically spaced entries 00106 * 1, ..., ULP and random complex angles. 00107 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00108 * and random complex angles. 00109 * 00110 * (7) Same as (4), but multiplied by SQRT( overflow threshold ) 00111 * (8) Same as (4), but multiplied by SQRT( underflow threshold ) 00112 * 00113 * (9) A matrix of the form U' T U, where U is unitary and 00114 * T has evenly spaced entries 1, ..., ULP with random complex 00115 * angles on the diagonal and random O(1) entries in the upper 00116 * triangle. 00117 * 00118 * (10) A matrix of the form U' T U, where U is unitary and 00119 * T has geometrically spaced entries 1, ..., ULP with random 00120 * complex angles on the diagonal and random O(1) entries in 00121 * the upper triangle. 00122 * 00123 * (11) A matrix of the form U' T U, where U is unitary and 00124 * T has "clustered" entries 1, ULP,..., ULP with random 00125 * complex angles on the diagonal and random O(1) entries in 00126 * the upper triangle. 00127 * 00128 * (12) A matrix of the form U' T U, where U is unitary and 00129 * T has complex eigenvalues randomly chosen from 00130 * ULP < |z| < 1 and random O(1) entries in the upper 00131 * triangle. 00132 * 00133 * (13) A matrix of the form X' T X, where X has condition 00134 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP 00135 * with random complex angles on the diagonal and random O(1) 00136 * entries in the upper triangle. 00137 * 00138 * (14) A matrix of the form X' T X, where X has condition 00139 * SQRT( ULP ) and T has geometrically spaced entries 00140 * 1, ..., ULP with random complex angles on the diagonal 00141 * and random O(1) entries in the upper triangle. 00142 * 00143 * (15) A matrix of the form X' T X, where X has condition 00144 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP 00145 * with random complex angles on the diagonal and random O(1) 00146 * entries in the upper triangle. 00147 * 00148 * (16) A matrix of the form X' T X, where X has condition 00149 * SQRT( ULP ) and T has complex eigenvalues randomly chosen 00150 * from ULP < |z| < 1 and random O(1) entries in the upper 00151 * triangle. 00152 * 00153 * (17) Same as (16), but multiplied by SQRT( overflow threshold ) 00154 * (18) Same as (16), but multiplied by SQRT( underflow threshold ) 00155 * 00156 * (19) Nonsymmetric matrix with random entries chosen from |z| < 1 00157 * (20) Same as (19), but multiplied by SQRT( overflow threshold ) 00158 * (21) Same as (19), but multiplied by SQRT( underflow threshold ) 00159 * 00160 * Arguments 00161 * ========== 00162 * 00163 * NSIZES - INTEGER 00164 * The number of sizes of matrices to use. If it is zero, 00165 * CCHKHS does nothing. It must be at least zero. 00166 * Not modified. 00167 * 00168 * NN - INTEGER array, dimension (NSIZES) 00169 * An array containing the sizes to be used for the matrices. 00170 * Zero values will be skipped. The values must be at least 00171 * zero. 00172 * Not modified. 00173 * 00174 * NTYPES - INTEGER 00175 * The number of elements in DOTYPE. If it is zero, CCHKHS 00176 * does nothing. It must be at least zero. If it is MAXTYP+1 00177 * and NSIZES is 1, then an additional type, MAXTYP+1 is 00178 * defined, which is to use whatever matrix is in A. This 00179 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00180 * DOTYPE(MAXTYP+1) is .TRUE. . 00181 * Not modified. 00182 * 00183 * DOTYPE - LOGICAL array, dimension (NTYPES) 00184 * If DOTYPE(j) is .TRUE., then for each size in NN a 00185 * matrix of that size and of type j will be generated. 00186 * If NTYPES is smaller than the maximum number of types 00187 * defined (PARAMETER MAXTYP), then types NTYPES+1 through 00188 * MAXTYP will not be generated. If NTYPES is larger 00189 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 00190 * will be ignored. 00191 * Not modified. 00192 * 00193 * ISEED - INTEGER array, dimension (4) 00194 * On entry ISEED specifies the seed of the random number 00195 * generator. The array elements should be between 0 and 4095; 00196 * if not they will be reduced mod 4096. Also, ISEED(4) must 00197 * be odd. The random number generator uses a linear 00198 * congruential sequence limited to small integers, and so 00199 * should produce machine independent random numbers. The 00200 * values of ISEED are changed on exit, and can be used in the 00201 * next call to CCHKHS to continue the same random number 00202 * sequence. 00203 * Modified. 00204 * 00205 * THRESH - REAL 00206 * A test will count as "failed" if the "error", computed as 00207 * described above, exceeds THRESH. Note that the error 00208 * is scaled to be O(1), so THRESH should be a reasonably 00209 * small multiple of 1, e.g., 10 or 100. In particular, 00210 * it should not depend on the precision (single vs. double) 00211 * or the size of the matrix. It must be at least zero. 00212 * Not modified. 00213 * 00214 * NOUNIT - INTEGER 00215 * The FORTRAN unit number for printing out error messages 00216 * (e.g., if a routine returns IINFO not equal to 0.) 00217 * Not modified. 00218 * 00219 * A - COMPLEX array, dimension (LDA,max(NN)) 00220 * Used to hold the matrix whose eigenvalues are to be 00221 * computed. On exit, A contains the last matrix actually 00222 * used. 00223 * Modified. 00224 * 00225 * LDA - INTEGER 00226 * The leading dimension of A, H, T1 and T2. It must be at 00227 * least 1 and at least max( NN ). 00228 * Not modified. 00229 * 00230 * H - COMPLEX array, dimension (LDA,max(NN)) 00231 * The upper hessenberg matrix computed by CGEHRD. On exit, 00232 * H contains the Hessenberg form of the matrix in A. 00233 * Modified. 00234 * 00235 * T1 - COMPLEX array, dimension (LDA,max(NN)) 00236 * The Schur (="quasi-triangular") matrix computed by CHSEQR 00237 * if Z is computed. On exit, T1 contains the Schur form of 00238 * the matrix in A. 00239 * Modified. 00240 * 00241 * T2 - COMPLEX array, dimension (LDA,max(NN)) 00242 * The Schur matrix computed by CHSEQR when Z is not computed. 00243 * This should be identical to T1. 00244 * Modified. 00245 * 00246 * LDU - INTEGER 00247 * The leading dimension of U, Z, UZ and UU. It must be at 00248 * least 1 and at least max( NN ). 00249 * Not modified. 00250 * 00251 * U - COMPLEX array, dimension (LDU,max(NN)) 00252 * The unitary matrix computed by CGEHRD. 00253 * Modified. 00254 * 00255 * Z - COMPLEX array, dimension (LDU,max(NN)) 00256 * The unitary matrix computed by CHSEQR. 00257 * Modified. 00258 * 00259 * UZ - COMPLEX array, dimension (LDU,max(NN)) 00260 * The product of U times Z. 00261 * Modified. 00262 * 00263 * W1 - COMPLEX array, dimension (max(NN)) 00264 * The eigenvalues of A, as computed by a full Schur 00265 * decomposition H = Z T Z'. On exit, W1 contains the 00266 * eigenvalues of the matrix in A. 00267 * Modified. 00268 * 00269 * W3 - COMPLEX array, dimension (max(NN)) 00270 * The eigenvalues of A, as computed by a partial Schur 00271 * decomposition (Z not computed, T only computed as much 00272 * as is necessary for determining eigenvalues). On exit, 00273 * W3 contains the eigenvalues of the matrix in A, possibly 00274 * perturbed by CHSEIN. 00275 * Modified. 00276 * 00277 * EVECTL - COMPLEX array, dimension (LDU,max(NN)) 00278 * The conjugate transpose of the (upper triangular) left 00279 * eigenvector matrix for the matrix in T1. 00280 * Modified. 00281 * 00282 * EVECTR - COMPLEX array, dimension (LDU,max(NN)) 00283 * The (upper triangular) right eigenvector matrix for the 00284 * matrix in T1. 00285 * Modified. 00286 * 00287 * EVECTY - COMPLEX array, dimension (LDU,max(NN)) 00288 * The conjugate transpose of the left eigenvector matrix 00289 * for the matrix in H. 00290 * Modified. 00291 * 00292 * EVECTX - COMPLEX array, dimension (LDU,max(NN)) 00293 * The right eigenvector matrix for the matrix in H. 00294 * Modified. 00295 * 00296 * UU - COMPLEX array, dimension (LDU,max(NN)) 00297 * Details of the unitary matrix computed by CGEHRD. 00298 * Modified. 00299 * 00300 * TAU - COMPLEX array, dimension (max(NN)) 00301 * Further details of the unitary matrix computed by CGEHRD. 00302 * Modified. 00303 * 00304 * WORK - COMPLEX array, dimension (NWORK) 00305 * Workspace. 00306 * Modified. 00307 * 00308 * NWORK - INTEGER 00309 * The number of entries in WORK. NWORK >= 4*NN(j)*NN(j) + 2. 00310 * 00311 * RWORK - REAL array, dimension (max(NN)) 00312 * Workspace. Could be equivalenced to IWORK, but not SELECT. 00313 * Modified. 00314 * 00315 * IWORK - INTEGER array, dimension (max(NN)) 00316 * Workspace. 00317 * Modified. 00318 * 00319 * SELECT - LOGICAL array, dimension (max(NN)) 00320 * Workspace. Could be equivalenced to IWORK, but not RWORK. 00321 * Modified. 00322 * 00323 * RESULT - REAL array, dimension (14) 00324 * The values computed by the fourteen tests described above. 00325 * The values are currently limited to 1/ulp, to avoid 00326 * overflow. 00327 * Modified. 00328 * 00329 * INFO - INTEGER 00330 * If 0, then everything ran OK. 00331 * -1: NSIZES < 0 00332 * -2: Some NN(j) < 0 00333 * -3: NTYPES < 0 00334 * -6: THRESH < 0 00335 * -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). 00336 * -14: LDU < 1 or LDU < NMAX. 00337 * -26: NWORK too small. 00338 * If CLATMR, CLATMS, or CLATME returns an error code, the 00339 * absolute value of it is returned. 00340 * If 1, then CHSEQR could not find all the shifts. 00341 * If 2, then the EISPACK code (for small blocks) failed. 00342 * If >2, then 30*N iterations were not enough to find an 00343 * eigenvalue or to decompose the problem. 00344 * Modified. 00345 * 00346 *----------------------------------------------------------------------- 00347 * 00348 * Some Local Variables and Parameters: 00349 * ---- ----- --------- --- ---------- 00350 * 00351 * ZERO, ONE Real 0 and 1. 00352 * MAXTYP The number of types defined. 00353 * MTEST The number of tests defined: care must be taken 00354 * that (1) the size of RESULT, (2) the number of 00355 * tests actually performed, and (3) MTEST agree. 00356 * NTEST The number of tests performed on this matrix 00357 * so far. This should be less than MTEST, and 00358 * equal to it by the last test. It will be less 00359 * if any of the routines being tested indicates 00360 * that it could not compute the matrices that 00361 * would be tested. 00362 * NMAX Largest value in NN. 00363 * NMATS The number of matrices generated so far. 00364 * NERRS The number of tests which have exceeded THRESH 00365 * so far (computed by SLAFTS). 00366 * COND, CONDS, 00367 * IMODE Values to be passed to the matrix generators. 00368 * ANORM Norm of A; passed to matrix generators. 00369 * 00370 * OVFL, UNFL Overflow and underflow thresholds. 00371 * ULP, ULPINV Finest relative precision and its inverse. 00372 * RTOVFL, RTUNFL, 00373 * RTULP, RTULPI Square roots of the previous 4 values. 00374 * 00375 * The following four arrays decode JTYPE: 00376 * KTYPE(j) The general type (1-10) for type "j". 00377 * KMODE(j) The MODE value to be passed to the matrix 00378 * generator for type "j". 00379 * KMAGN(j) The order of magnitude ( O(1), 00380 * O(overflow^(1/2) ), O(underflow^(1/2) ) 00381 * KCONDS(j) Selects whether CONDS is to be 1 or 00382 * 1/sqrt(ulp). (0 means irrelevant.) 00383 * 00384 * ===================================================================== 00385 * 00386 * .. Parameters .. 00387 REAL ZERO, ONE 00388 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00389 COMPLEX CZERO, CONE 00390 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00391 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00392 INTEGER MAXTYP 00393 PARAMETER ( MAXTYP = 21 ) 00394 * .. 00395 * .. Local Scalars .. 00396 LOGICAL BADNN, MATCH 00397 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL, 00398 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS, 00399 $ NMATS, NMAX, NTEST, NTESTT 00400 REAL ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP, 00401 $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL 00402 * .. 00403 * .. Local Arrays .. 00404 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), 00405 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 00406 $ KTYPE( MAXTYP ) 00407 REAL DUMMA( 4 ) 00408 COMPLEX CDUMMA( 4 ) 00409 * .. 00410 * .. External Functions .. 00411 REAL SLAMCH 00412 EXTERNAL SLAMCH 00413 * .. 00414 * .. External Subroutines .. 00415 EXTERNAL CCOPY, CGEHRD, CGEMM, CGET10, CGET22, CHSEIN, 00416 $ CHSEQR, CHST01, CLACPY, CLASET, CLATME, CLATMR, 00417 $ CLATMS, CTREVC, CUNGHR, CUNMHR, SLABAD, SLAFTS, 00418 $ SLASUM, XERBLA 00419 * .. 00420 * .. Intrinsic Functions .. 00421 INTRINSIC ABS, MAX, MIN, REAL, SQRT 00422 * .. 00423 * .. Data statements .. 00424 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 00425 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 00426 $ 3, 1, 2, 3 / 00427 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 00428 $ 1, 5, 5, 5, 4, 3, 1 / 00429 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 00430 * .. 00431 * .. Executable Statements .. 00432 * 00433 * Check for errors 00434 * 00435 NTESTT = 0 00436 INFO = 0 00437 * 00438 BADNN = .FALSE. 00439 NMAX = 0 00440 DO 10 J = 1, NSIZES 00441 NMAX = MAX( NMAX, NN( J ) ) 00442 IF( NN( J ).LT.0 ) 00443 $ BADNN = .TRUE. 00444 10 CONTINUE 00445 * 00446 * Check for errors 00447 * 00448 IF( NSIZES.LT.0 ) THEN 00449 INFO = -1 00450 ELSE IF( BADNN ) THEN 00451 INFO = -2 00452 ELSE IF( NTYPES.LT.0 ) THEN 00453 INFO = -3 00454 ELSE IF( THRESH.LT.ZERO ) THEN 00455 INFO = -6 00456 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 00457 INFO = -9 00458 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN 00459 INFO = -14 00460 ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN 00461 INFO = -26 00462 END IF 00463 * 00464 IF( INFO.NE.0 ) THEN 00465 CALL XERBLA( 'CCHKHS', -INFO ) 00466 RETURN 00467 END IF 00468 * 00469 * Quick return if possible 00470 * 00471 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00472 $ RETURN 00473 * 00474 * More important constants 00475 * 00476 UNFL = SLAMCH( 'Safe minimum' ) 00477 OVFL = SLAMCH( 'Overflow' ) 00478 CALL SLABAD( UNFL, OVFL ) 00479 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00480 ULPINV = ONE / ULP 00481 RTUNFL = SQRT( UNFL ) 00482 RTOVFL = SQRT( OVFL ) 00483 RTULP = SQRT( ULP ) 00484 RTULPI = ONE / RTULP 00485 * 00486 * Loop over sizes, types 00487 * 00488 NERRS = 0 00489 NMATS = 0 00490 * 00491 DO 260 JSIZE = 1, NSIZES 00492 N = NN( JSIZE ) 00493 IF( N.EQ.0 ) 00494 $ GO TO 260 00495 N1 = MAX( 1, N ) 00496 ANINV = ONE / REAL( N1 ) 00497 * 00498 IF( NSIZES.NE.1 ) THEN 00499 MTYPES = MIN( MAXTYP, NTYPES ) 00500 ELSE 00501 MTYPES = MIN( MAXTYP+1, NTYPES ) 00502 END IF 00503 * 00504 DO 250 JTYPE = 1, MTYPES 00505 IF( .NOT.DOTYPE( JTYPE ) ) 00506 $ GO TO 250 00507 NMATS = NMATS + 1 00508 NTEST = 0 00509 * 00510 * Save ISEED in case of an error. 00511 * 00512 DO 20 J = 1, 4 00513 IOLDSD( J ) = ISEED( J ) 00514 20 CONTINUE 00515 * 00516 * Initialize RESULT 00517 * 00518 DO 30 J = 1, 14 00519 RESULT( J ) = ZERO 00520 30 CONTINUE 00521 * 00522 * Compute "A" 00523 * 00524 * Control parameters: 00525 * 00526 * KMAGN KCONDS KMODE KTYPE 00527 * =1 O(1) 1 clustered 1 zero 00528 * =2 large large clustered 2 identity 00529 * =3 small exponential Jordan 00530 * =4 arithmetic diagonal, (w/ eigenvalues) 00531 * =5 random log hermitian, w/ eigenvalues 00532 * =6 random general, w/ eigenvalues 00533 * =7 random diagonal 00534 * =8 random hermitian 00535 * =9 random general 00536 * =10 random triangular 00537 * 00538 IF( MTYPES.GT.MAXTYP ) 00539 $ GO TO 100 00540 * 00541 ITYPE = KTYPE( JTYPE ) 00542 IMODE = KMODE( JTYPE ) 00543 * 00544 * Compute norm 00545 * 00546 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 00547 * 00548 40 CONTINUE 00549 ANORM = ONE 00550 GO TO 70 00551 * 00552 50 CONTINUE 00553 ANORM = ( RTOVFL*ULP )*ANINV 00554 GO TO 70 00555 * 00556 60 CONTINUE 00557 ANORM = RTUNFL*N*ULPINV 00558 GO TO 70 00559 * 00560 70 CONTINUE 00561 * 00562 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA ) 00563 IINFO = 0 00564 COND = ULPINV 00565 * 00566 * Special Matrices 00567 * 00568 IF( ITYPE.EQ.1 ) THEN 00569 * 00570 * Zero 00571 * 00572 IINFO = 0 00573 ELSE IF( ITYPE.EQ.2 ) THEN 00574 * 00575 * Identity 00576 * 00577 DO 80 JCOL = 1, N 00578 A( JCOL, JCOL ) = ANORM 00579 80 CONTINUE 00580 * 00581 ELSE IF( ITYPE.EQ.3 ) THEN 00582 * 00583 * Jordan Block 00584 * 00585 DO 90 JCOL = 1, N 00586 A( JCOL, JCOL ) = ANORM 00587 IF( JCOL.GT.1 ) 00588 $ A( JCOL, JCOL-1 ) = ONE 00589 90 CONTINUE 00590 * 00591 ELSE IF( ITYPE.EQ.4 ) THEN 00592 * 00593 * Diagonal Matrix, [Eigen]values Specified 00594 * 00595 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, IMODE, COND, 00596 $ CONE, 'T', 'N', WORK( N+1 ), 1, ONE, 00597 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00598 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00599 * 00600 ELSE IF( ITYPE.EQ.5 ) THEN 00601 * 00602 * Hermitian, eigenvalues specified 00603 * 00604 CALL CLATMS( N, N, 'D', ISEED, 'H', RWORK, IMODE, COND, 00605 $ ANORM, N, N, 'N', A, LDA, WORK, IINFO ) 00606 * 00607 ELSE IF( ITYPE.EQ.6 ) THEN 00608 * 00609 * General, eigenvalues specified 00610 * 00611 IF( KCONDS( JTYPE ).EQ.1 ) THEN 00612 CONDS = ONE 00613 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 00614 CONDS = RTULPI 00615 ELSE 00616 CONDS = ZERO 00617 END IF 00618 * 00619 CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ', 00620 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM, 00621 $ A, LDA, WORK( N+1 ), IINFO ) 00622 * 00623 ELSE IF( ITYPE.EQ.7 ) THEN 00624 * 00625 * Diagonal, random eigenvalues 00626 * 00627 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00628 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00629 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00630 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00631 * 00632 ELSE IF( ITYPE.EQ.8 ) THEN 00633 * 00634 * Hermitian, random eigenvalues 00635 * 00636 CALL CLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE, 00637 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00638 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00639 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00640 * 00641 ELSE IF( ITYPE.EQ.9 ) THEN 00642 * 00643 * General, random eigenvalues 00644 * 00645 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00646 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00647 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00648 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00649 * 00650 ELSE IF( ITYPE.EQ.10 ) THEN 00651 * 00652 * Triangular, random eigenvalues 00653 * 00654 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00655 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00656 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 00657 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00658 * 00659 ELSE 00660 * 00661 IINFO = 1 00662 END IF 00663 * 00664 IF( IINFO.NE.0 ) THEN 00665 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 00666 $ IOLDSD 00667 INFO = ABS( IINFO ) 00668 RETURN 00669 END IF 00670 * 00671 100 CONTINUE 00672 * 00673 * Call CGEHRD to compute H and U, do tests. 00674 * 00675 CALL CLACPY( ' ', N, N, A, LDA, H, LDA ) 00676 NTEST = 1 00677 * 00678 ILO = 1 00679 IHI = N 00680 * 00681 CALL CGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ), 00682 $ NWORK-N, IINFO ) 00683 * 00684 IF( IINFO.NE.0 ) THEN 00685 RESULT( 1 ) = ULPINV 00686 WRITE( NOUNIT, FMT = 9999 )'CGEHRD', IINFO, N, JTYPE, 00687 $ IOLDSD 00688 INFO = ABS( IINFO ) 00689 GO TO 240 00690 END IF 00691 * 00692 DO 120 J = 1, N - 1 00693 UU( J+1, J ) = CZERO 00694 DO 110 I = J + 2, N 00695 U( I, J ) = H( I, J ) 00696 UU( I, J ) = H( I, J ) 00697 H( I, J ) = CZERO 00698 110 CONTINUE 00699 120 CONTINUE 00700 CALL CCOPY( N-1, WORK, 1, TAU, 1 ) 00701 CALL CUNGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ), 00702 $ NWORK-N, IINFO ) 00703 NTEST = 2 00704 * 00705 CALL CHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK, 00706 $ NWORK, RWORK, RESULT( 1 ) ) 00707 * 00708 * Call CHSEQR to compute T1, T2 and Z, do tests. 00709 * 00710 * Eigenvalues only (W3) 00711 * 00712 CALL CLACPY( ' ', N, N, H, LDA, T2, LDA ) 00713 NTEST = 3 00714 RESULT( 3 ) = ULPINV 00715 * 00716 CALL CHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, W3, UZ, LDU, 00717 $ WORK, NWORK, IINFO ) 00718 IF( IINFO.NE.0 ) THEN 00719 WRITE( NOUNIT, FMT = 9999 )'CHSEQR(E)', IINFO, N, JTYPE, 00720 $ IOLDSD 00721 IF( IINFO.LE.N+2 ) THEN 00722 INFO = ABS( IINFO ) 00723 GO TO 240 00724 END IF 00725 END IF 00726 * 00727 * Eigenvalues (W1) and Full Schur Form (T2) 00728 * 00729 CALL CLACPY( ' ', N, N, H, LDA, T2, LDA ) 00730 * 00731 CALL CHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, W1, UZ, LDU, 00732 $ WORK, NWORK, IINFO ) 00733 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 00734 WRITE( NOUNIT, FMT = 9999 )'CHSEQR(S)', IINFO, N, JTYPE, 00735 $ IOLDSD 00736 INFO = ABS( IINFO ) 00737 GO TO 240 00738 END IF 00739 * 00740 * Eigenvalues (W1), Schur Form (T1), and Schur Vectors (UZ) 00741 * 00742 CALL CLACPY( ' ', N, N, H, LDA, T1, LDA ) 00743 CALL CLACPY( ' ', N, N, U, LDU, UZ, LDU ) 00744 * 00745 CALL CHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, W1, UZ, LDU, 00746 $ WORK, NWORK, IINFO ) 00747 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 00748 WRITE( NOUNIT, FMT = 9999 )'CHSEQR(V)', IINFO, N, JTYPE, 00749 $ IOLDSD 00750 INFO = ABS( IINFO ) 00751 GO TO 240 00752 END IF 00753 * 00754 * Compute Z = U' UZ 00755 * 00756 CALL CGEMM( 'C', 'N', N, N, N, CONE, U, LDU, UZ, LDU, CZERO, 00757 $ Z, LDU ) 00758 NTEST = 8 00759 * 00760 * Do Tests 3: | H - Z T Z' | / ( |H| n ulp ) 00761 * and 4: | I - Z Z' | / ( n ulp ) 00762 * 00763 CALL CHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK, 00764 $ NWORK, RWORK, RESULT( 3 ) ) 00765 * 00766 * Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp ) 00767 * and 6: | I - UZ (UZ)' | / ( n ulp ) 00768 * 00769 CALL CHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK, 00770 $ NWORK, RWORK, RESULT( 5 ) ) 00771 * 00772 * Do Test 7: | T2 - T1 | / ( |T| n ulp ) 00773 * 00774 CALL CGET10( N, N, T2, LDA, T1, LDA, WORK, RWORK, 00775 $ RESULT( 7 ) ) 00776 * 00777 * Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp ) 00778 * 00779 TEMP1 = ZERO 00780 TEMP2 = ZERO 00781 DO 130 J = 1, N 00782 TEMP1 = MAX( TEMP1, ABS( W1( J ) ), ABS( W3( J ) ) ) 00783 TEMP2 = MAX( TEMP2, ABS( W1( J )-W3( J ) ) ) 00784 130 CONTINUE 00785 * 00786 RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 00787 * 00788 * Compute the Left and Right Eigenvectors of T 00789 * 00790 * Compute the Right eigenvector Matrix: 00791 * 00792 NTEST = 9 00793 RESULT( 9 ) = ULPINV 00794 * 00795 * Select every other eigenvector 00796 * 00797 DO 140 J = 1, N 00798 SELECT( J ) = .FALSE. 00799 140 CONTINUE 00800 DO 150 J = 1, N, 2 00801 SELECT( J ) = .TRUE. 00802 150 CONTINUE 00803 CALL CTREVC( 'Right', 'All', SELECT, N, T1, LDA, CDUMMA, 00804 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 00805 IF( IINFO.NE.0 ) THEN 00806 WRITE( NOUNIT, FMT = 9999 )'CTREVC(R,A)', IINFO, N, 00807 $ JTYPE, IOLDSD 00808 INFO = ABS( IINFO ) 00809 GO TO 240 00810 END IF 00811 * 00812 * Test 9: | TR - RW | / ( |T| |R| ulp ) 00813 * 00814 CALL CGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, W1, 00815 $ WORK, RWORK, DUMMA( 1 ) ) 00816 RESULT( 9 ) = DUMMA( 1 ) 00817 IF( DUMMA( 2 ).GT.THRESH ) THEN 00818 WRITE( NOUNIT, FMT = 9998 )'Right', 'CTREVC', 00819 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00820 END IF 00821 * 00822 * Compute selected right eigenvectors and confirm that 00823 * they agree with previous right eigenvectors 00824 * 00825 CALL CTREVC( 'Right', 'Some', SELECT, N, T1, LDA, CDUMMA, 00826 $ LDU, EVECTL, LDU, N, IN, WORK, RWORK, IINFO ) 00827 IF( IINFO.NE.0 ) THEN 00828 WRITE( NOUNIT, FMT = 9999 )'CTREVC(R,S)', IINFO, N, 00829 $ JTYPE, IOLDSD 00830 INFO = ABS( IINFO ) 00831 GO TO 240 00832 END IF 00833 * 00834 K = 1 00835 MATCH = .TRUE. 00836 DO 170 J = 1, N 00837 IF( SELECT( J ) ) THEN 00838 DO 160 JJ = 1, N 00839 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN 00840 MATCH = .FALSE. 00841 GO TO 180 00842 END IF 00843 160 CONTINUE 00844 K = K + 1 00845 END IF 00846 170 CONTINUE 00847 180 CONTINUE 00848 IF( .NOT.MATCH ) 00849 $ WRITE( NOUNIT, FMT = 9997 )'Right', 'CTREVC', N, JTYPE, 00850 $ IOLDSD 00851 * 00852 * Compute the Left eigenvector Matrix: 00853 * 00854 NTEST = 10 00855 RESULT( 10 ) = ULPINV 00856 CALL CTREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU, 00857 $ CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 00858 IF( IINFO.NE.0 ) THEN 00859 WRITE( NOUNIT, FMT = 9999 )'CTREVC(L,A)', IINFO, N, 00860 $ JTYPE, IOLDSD 00861 INFO = ABS( IINFO ) 00862 GO TO 240 00863 END IF 00864 * 00865 * Test 10: | LT - WL | / ( |T| |L| ulp ) 00866 * 00867 CALL CGET22( 'C', 'N', 'C', N, T1, LDA, EVECTL, LDU, W1, 00868 $ WORK, RWORK, DUMMA( 3 ) ) 00869 RESULT( 10 ) = DUMMA( 3 ) 00870 IF( DUMMA( 4 ).GT.THRESH ) THEN 00871 WRITE( NOUNIT, FMT = 9998 )'Left', 'CTREVC', DUMMA( 4 ), 00872 $ N, JTYPE, IOLDSD 00873 END IF 00874 * 00875 * Compute selected left eigenvectors and confirm that 00876 * they agree with previous left eigenvectors 00877 * 00878 CALL CTREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR, 00879 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 00880 IF( IINFO.NE.0 ) THEN 00881 WRITE( NOUNIT, FMT = 9999 )'CTREVC(L,S)', IINFO, N, 00882 $ JTYPE, IOLDSD 00883 INFO = ABS( IINFO ) 00884 GO TO 240 00885 END IF 00886 * 00887 K = 1 00888 MATCH = .TRUE. 00889 DO 200 J = 1, N 00890 IF( SELECT( J ) ) THEN 00891 DO 190 JJ = 1, N 00892 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN 00893 MATCH = .FALSE. 00894 GO TO 210 00895 END IF 00896 190 CONTINUE 00897 K = K + 1 00898 END IF 00899 200 CONTINUE 00900 210 CONTINUE 00901 IF( .NOT.MATCH ) 00902 $ WRITE( NOUNIT, FMT = 9997 )'Left', 'CTREVC', N, JTYPE, 00903 $ IOLDSD 00904 * 00905 * Call CHSEIN for Right eigenvectors of H, do test 11 00906 * 00907 NTEST = 11 00908 RESULT( 11 ) = ULPINV 00909 DO 220 J = 1, N 00910 SELECT( J ) = .TRUE. 00911 220 CONTINUE 00912 * 00913 CALL CHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3, 00914 $ CDUMMA, LDU, EVECTX, LDU, N1, IN, WORK, RWORK, 00915 $ IWORK, IWORK, IINFO ) 00916 IF( IINFO.NE.0 ) THEN 00917 WRITE( NOUNIT, FMT = 9999 )'CHSEIN(R)', IINFO, N, JTYPE, 00918 $ IOLDSD 00919 INFO = ABS( IINFO ) 00920 IF( IINFO.LT.0 ) 00921 $ GO TO 240 00922 ELSE 00923 * 00924 * Test 11: | HX - XW | / ( |H| |X| ulp ) 00925 * 00926 * (from inverse iteration) 00927 * 00928 CALL CGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, W3, 00929 $ WORK, RWORK, DUMMA( 1 ) ) 00930 IF( DUMMA( 1 ).LT.ULPINV ) 00931 $ RESULT( 11 ) = DUMMA( 1 )*ANINV 00932 IF( DUMMA( 2 ).GT.THRESH ) THEN 00933 WRITE( NOUNIT, FMT = 9998 )'Right', 'CHSEIN', 00934 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00935 END IF 00936 END IF 00937 * 00938 * Call CHSEIN for Left eigenvectors of H, do test 12 00939 * 00940 NTEST = 12 00941 RESULT( 12 ) = ULPINV 00942 DO 230 J = 1, N 00943 SELECT( J ) = .TRUE. 00944 230 CONTINUE 00945 * 00946 CALL CHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3, 00947 $ EVECTY, LDU, CDUMMA, LDU, N1, IN, WORK, RWORK, 00948 $ IWORK, IWORK, IINFO ) 00949 IF( IINFO.NE.0 ) THEN 00950 WRITE( NOUNIT, FMT = 9999 )'CHSEIN(L)', IINFO, N, JTYPE, 00951 $ IOLDSD 00952 INFO = ABS( IINFO ) 00953 IF( IINFO.LT.0 ) 00954 $ GO TO 240 00955 ELSE 00956 * 00957 * Test 12: | YH - WY | / ( |H| |Y| ulp ) 00958 * 00959 * (from inverse iteration) 00960 * 00961 CALL CGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, W3, 00962 $ WORK, RWORK, DUMMA( 3 ) ) 00963 IF( DUMMA( 3 ).LT.ULPINV ) 00964 $ RESULT( 12 ) = DUMMA( 3 )*ANINV 00965 IF( DUMMA( 4 ).GT.THRESH ) THEN 00966 WRITE( NOUNIT, FMT = 9998 )'Left', 'CHSEIN', 00967 $ DUMMA( 4 ), N, JTYPE, IOLDSD 00968 END IF 00969 END IF 00970 * 00971 * Call CUNMHR for Right eigenvectors of A, do test 13 00972 * 00973 NTEST = 13 00974 RESULT( 13 ) = ULPINV 00975 * 00976 CALL CUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 00977 $ LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO ) 00978 IF( IINFO.NE.0 ) THEN 00979 WRITE( NOUNIT, FMT = 9999 )'CUNMHR(L)', IINFO, N, JTYPE, 00980 $ IOLDSD 00981 INFO = ABS( IINFO ) 00982 IF( IINFO.LT.0 ) 00983 $ GO TO 240 00984 ELSE 00985 * 00986 * Test 13: | AX - XW | / ( |A| |X| ulp ) 00987 * 00988 * (from inverse iteration) 00989 * 00990 CALL CGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, W3, 00991 $ WORK, RWORK, DUMMA( 1 ) ) 00992 IF( DUMMA( 1 ).LT.ULPINV ) 00993 $ RESULT( 13 ) = DUMMA( 1 )*ANINV 00994 END IF 00995 * 00996 * Call CUNMHR for Left eigenvectors of A, do test 14 00997 * 00998 NTEST = 14 00999 RESULT( 14 ) = ULPINV 01000 * 01001 CALL CUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 01002 $ LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO ) 01003 IF( IINFO.NE.0 ) THEN 01004 WRITE( NOUNIT, FMT = 9999 )'CUNMHR(L)', IINFO, N, JTYPE, 01005 $ IOLDSD 01006 INFO = ABS( IINFO ) 01007 IF( IINFO.LT.0 ) 01008 $ GO TO 240 01009 ELSE 01010 * 01011 * Test 14: | YA - WY | / ( |A| |Y| ulp ) 01012 * 01013 * (from inverse iteration) 01014 * 01015 CALL CGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, W3, 01016 $ WORK, RWORK, DUMMA( 3 ) ) 01017 IF( DUMMA( 3 ).LT.ULPINV ) 01018 $ RESULT( 14 ) = DUMMA( 3 )*ANINV 01019 END IF 01020 * 01021 * End of Loop -- Check for RESULT(j) > THRESH 01022 * 01023 240 CONTINUE 01024 * 01025 NTESTT = NTESTT + NTEST 01026 CALL SLAFTS( 'CHS', N, N, JTYPE, NTEST, RESULT, IOLDSD, 01027 $ THRESH, NOUNIT, NERRS ) 01028 * 01029 250 CONTINUE 01030 260 CONTINUE 01031 * 01032 * Summary 01033 * 01034 CALL SLASUM( 'CHS', NOUNIT, NERRS, NTESTT ) 01035 * 01036 RETURN 01037 * 01038 9999 FORMAT( ' CCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 01039 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 01040 9998 FORMAT( ' CCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ', 01041 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 01042 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 01043 $ ')' ) 01044 9997 FORMAT( ' CCHKHS: Selected ', A, ' Eigenvectors from ', A, 01045 $ ' do not match other eigenvectors ', 9X, 'N=', I6, 01046 $ ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 01047 * 01048 * End of CCHKHS 01049 * 01050 END
1.7.3 
