|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZCHKHS( 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.1.1) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * February 2007 00010 * 00011 * .. Scalar Arguments .. 00012 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK 00013 DOUBLE PRECISION THRESH 00014 * .. 00015 * .. Array Arguments .. 00016 LOGICAL DOTYPE( * ), SELECT( * ) 00017 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00018 DOUBLE PRECISION RESULT( 14 ), RWORK( * ) 00019 COMPLEX*16 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 * ZCHKHS checks the nonsymmetric eigenvalue problem routines. 00031 * 00032 * ZGEHRD factors A as U H U' , where ' means conjugate 00033 * transpose, H is hessenberg, and U is unitary. 00034 * 00035 * ZUNGHR generates the unitary matrix U. 00036 * 00037 * ZUNMHR multiplies a matrix by the unitary matrix U. 00038 * 00039 * ZHSEQR 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 * ZTREVC 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 * ZHSEIN 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 ZCHKHS 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 * ZCHKHS 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, ZCHKHS 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 ZCHKHS to continue the same random number 00202 * sequence. 00203 * Modified. 00204 * 00205 * THRESH - DOUBLE PRECISION 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*16 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*16 array, dimension (LDA,max(NN)) 00231 * The upper hessenberg matrix computed by ZGEHRD. On exit, 00232 * H contains the Hessenberg form of the matrix in A. 00233 * Modified. 00234 * 00235 * T1 - COMPLEX*16 array, dimension (LDA,max(NN)) 00236 * The Schur (="quasi-triangular") matrix computed by ZHSEQR 00237 * if Z is computed. On exit, T1 contains the Schur form of 00238 * the matrix in A. 00239 * Modified. 00240 * 00241 * T2 - COMPLEX*16 array, dimension (LDA,max(NN)) 00242 * The Schur matrix computed by ZHSEQR 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*16 array, dimension (LDU,max(NN)) 00252 * The unitary matrix computed by ZGEHRD. 00253 * Modified. 00254 * 00255 * Z - COMPLEX*16 array, dimension (LDU,max(NN)) 00256 * The unitary matrix computed by ZHSEQR. 00257 * Modified. 00258 * 00259 * UZ - COMPLEX*16 array, dimension (LDU,max(NN)) 00260 * The product of U times Z. 00261 * Modified. 00262 * 00263 * W1 - COMPLEX*16 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*16 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 ZHSEIN. 00275 * Modified. 00276 * 00277 * EVECTL - COMPLEX*16 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 * EVEZTR - COMPLEX*16 array, dimension (LDU,max(NN)) 00283 * The (upper triangular) right eigenvector matrix for the 00284 * matrix in T1. 00285 * Modified. 00286 * 00287 * EVECTY - COMPLEX*16 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*16 array, dimension (LDU,max(NN)) 00293 * The right eigenvector matrix for the matrix in H. 00294 * Modified. 00295 * 00296 * UU - COMPLEX*16 array, dimension (LDU,max(NN)) 00297 * Details of the unitary matrix computed by ZGEHRD. 00298 * Modified. 00299 * 00300 * TAU - COMPLEX*16 array, dimension (max(NN)) 00301 * Further details of the unitary matrix computed by ZGEHRD. 00302 * Modified. 00303 * 00304 * WORK - COMPLEX*16 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 ZLATMR, CLATMS, or CLATME returns an error code, the 00339 * absolute value of it is returned. 00340 * If 1, then ZHSEQR 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 DLAFTS). 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 DOUBLE PRECISION ZERO, ONE 00388 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00389 COMPLEX*16 CZERO, CONE 00390 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00391 $ CONE = ( 1.0D+0, 0.0D+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 DOUBLE PRECISION 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 DOUBLE PRECISION DUMMA( 4 ) 00408 COMPLEX*16 CDUMMA( 4 ) 00409 * .. 00410 * .. External Functions .. 00411 DOUBLE PRECISION DLAMCH 00412 EXTERNAL DLAMCH 00413 * .. 00414 * .. External Subroutines .. 00415 EXTERNAL DLABAD, DLAFTS, DLASUM, XERBLA, ZCOPY, ZGEHRD, 00416 $ ZGEMM, ZGET10, ZGET22, ZHSEIN, ZHSEQR, ZHST01, 00417 $ ZLACPY, ZLASET, ZLATME, ZLATMR, ZLATMS, ZTREVC, 00418 $ ZUNGHR, ZUNMHR 00419 * .. 00420 * .. Intrinsic Functions .. 00421 INTRINSIC ABS, DBLE, MAX, MIN, 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( 'ZCHKHS', -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 = DLAMCH( 'Safe minimum' ) 00477 OVFL = DLAMCH( 'Overflow' ) 00478 CALL DLABAD( UNFL, OVFL ) 00479 ULP = DLAMCH( 'Epsilon' )*DLAMCH( '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 N1 = MAX( 1, N ) 00494 ANINV = ONE / DBLE( N1 ) 00495 * 00496 IF( NSIZES.NE.1 ) THEN 00497 MTYPES = MIN( MAXTYP, NTYPES ) 00498 ELSE 00499 MTYPES = MIN( MAXTYP+1, NTYPES ) 00500 END IF 00501 * 00502 DO 250 JTYPE = 1, MTYPES 00503 IF( .NOT.DOTYPE( JTYPE ) ) 00504 $ GO TO 250 00505 NMATS = NMATS + 1 00506 NTEST = 0 00507 * 00508 * Save ISEED in case of an error. 00509 * 00510 DO 20 J = 1, 4 00511 IOLDSD( J ) = ISEED( J ) 00512 20 CONTINUE 00513 * 00514 * Initialize RESULT 00515 * 00516 DO 30 J = 1, 14 00517 RESULT( J ) = ZERO 00518 30 CONTINUE 00519 * 00520 * Compute "A" 00521 * 00522 * Control parameters: 00523 * 00524 * KMAGN KCONDS KMODE KTYPE 00525 * =1 O(1) 1 clustered 1 zero 00526 * =2 large large clustered 2 identity 00527 * =3 small exponential Jordan 00528 * =4 arithmetic diagonal, (w/ eigenvalues) 00529 * =5 random log hermitian, w/ eigenvalues 00530 * =6 random general, w/ eigenvalues 00531 * =7 random diagonal 00532 * =8 random hermitian 00533 * =9 random general 00534 * =10 random triangular 00535 * 00536 IF( MTYPES.GT.MAXTYP ) 00537 $ GO TO 100 00538 * 00539 ITYPE = KTYPE( JTYPE ) 00540 IMODE = KMODE( JTYPE ) 00541 * 00542 * Compute norm 00543 * 00544 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 00545 * 00546 40 CONTINUE 00547 ANORM = ONE 00548 GO TO 70 00549 * 00550 50 CONTINUE 00551 ANORM = ( RTOVFL*ULP )*ANINV 00552 GO TO 70 00553 * 00554 60 CONTINUE 00555 ANORM = RTUNFL*N*ULPINV 00556 GO TO 70 00557 * 00558 70 CONTINUE 00559 * 00560 CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA ) 00561 IINFO = 0 00562 COND = ULPINV 00563 * 00564 * Special Matrices 00565 * 00566 IF( ITYPE.EQ.1 ) THEN 00567 * 00568 * Zero 00569 * 00570 IINFO = 0 00571 ELSE IF( ITYPE.EQ.2 ) THEN 00572 * 00573 * Identity 00574 * 00575 DO 80 JCOL = 1, N 00576 A( JCOL, JCOL ) = ANORM 00577 80 CONTINUE 00578 * 00579 ELSE IF( ITYPE.EQ.3 ) THEN 00580 * 00581 * Jordan Block 00582 * 00583 DO 90 JCOL = 1, N 00584 A( JCOL, JCOL ) = ANORM 00585 IF( JCOL.GT.1 ) 00586 $ A( JCOL, JCOL-1 ) = ONE 00587 90 CONTINUE 00588 * 00589 ELSE IF( ITYPE.EQ.4 ) THEN 00590 * 00591 * Diagonal Matrix, [Eigen]values Specified 00592 * 00593 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, IMODE, COND, 00594 $ CONE, 'T', 'N', WORK( N+1 ), 1, ONE, 00595 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00596 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00597 * 00598 ELSE IF( ITYPE.EQ.5 ) THEN 00599 * 00600 * Hermitian, eigenvalues specified 00601 * 00602 CALL ZLATMS( N, N, 'D', ISEED, 'H', RWORK, IMODE, COND, 00603 $ ANORM, N, N, 'N', A, LDA, WORK, IINFO ) 00604 * 00605 ELSE IF( ITYPE.EQ.6 ) THEN 00606 * 00607 * General, eigenvalues specified 00608 * 00609 IF( KCONDS( JTYPE ).EQ.1 ) THEN 00610 CONDS = ONE 00611 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 00612 CONDS = RTULPI 00613 ELSE 00614 CONDS = ZERO 00615 END IF 00616 * 00617 CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ', 00618 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM, 00619 $ A, LDA, WORK( N+1 ), IINFO ) 00620 * 00621 ELSE IF( ITYPE.EQ.7 ) THEN 00622 * 00623 * Diagonal, random eigenvalues 00624 * 00625 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00626 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00627 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00628 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00629 * 00630 ELSE IF( ITYPE.EQ.8 ) THEN 00631 * 00632 * Hermitian, random eigenvalues 00633 * 00634 CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE, 00635 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00636 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00637 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00638 * 00639 ELSE IF( ITYPE.EQ.9 ) THEN 00640 * 00641 * General, random eigenvalues 00642 * 00643 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00644 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00645 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00646 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00647 * 00648 ELSE IF( ITYPE.EQ.10 ) THEN 00649 * 00650 * Triangular, random eigenvalues 00651 * 00652 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE, 00653 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00654 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 00655 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00656 * 00657 ELSE 00658 * 00659 IINFO = 1 00660 END IF 00661 * 00662 IF( IINFO.NE.0 ) THEN 00663 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 00664 $ IOLDSD 00665 INFO = ABS( IINFO ) 00666 RETURN 00667 END IF 00668 * 00669 100 CONTINUE 00670 * 00671 * Call ZGEHRD to compute H and U, do tests. 00672 * 00673 CALL ZLACPY( ' ', N, N, A, LDA, H, LDA ) 00674 NTEST = 1 00675 * 00676 ILO = 1 00677 IHI = N 00678 * 00679 CALL ZGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ), 00680 $ NWORK-N, IINFO ) 00681 * 00682 IF( IINFO.NE.0 ) THEN 00683 RESULT( 1 ) = ULPINV 00684 WRITE( NOUNIT, FMT = 9999 )'ZGEHRD', IINFO, N, JTYPE, 00685 $ IOLDSD 00686 INFO = ABS( IINFO ) 00687 GO TO 240 00688 END IF 00689 * 00690 DO 120 J = 1, N - 1 00691 UU( J+1, J ) = CZERO 00692 DO 110 I = J + 2, N 00693 U( I, J ) = H( I, J ) 00694 UU( I, J ) = H( I, J ) 00695 H( I, J ) = CZERO 00696 110 CONTINUE 00697 120 CONTINUE 00698 CALL ZCOPY( N-1, WORK, 1, TAU, 1 ) 00699 CALL ZUNGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ), 00700 $ NWORK-N, IINFO ) 00701 NTEST = 2 00702 * 00703 CALL ZHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK, 00704 $ NWORK, RWORK, RESULT( 1 ) ) 00705 * 00706 * Call ZHSEQR to compute T1, T2 and Z, do tests. 00707 * 00708 * Eigenvalues only (W3) 00709 * 00710 CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA ) 00711 NTEST = 3 00712 RESULT( 3 ) = ULPINV 00713 * 00714 CALL ZHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, W3, UZ, LDU, 00715 $ WORK, NWORK, IINFO ) 00716 IF( IINFO.NE.0 ) THEN 00717 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(E)', IINFO, N, JTYPE, 00718 $ IOLDSD 00719 IF( IINFO.LE.N+2 ) THEN 00720 INFO = ABS( IINFO ) 00721 GO TO 240 00722 END IF 00723 END IF 00724 * 00725 * Eigenvalues (W1) and Full Schur Form (T2) 00726 * 00727 CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA ) 00728 * 00729 CALL ZHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, W1, UZ, LDU, 00730 $ WORK, NWORK, IINFO ) 00731 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 00732 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(S)', IINFO, N, JTYPE, 00733 $ IOLDSD 00734 INFO = ABS( IINFO ) 00735 GO TO 240 00736 END IF 00737 * 00738 * Eigenvalues (W1), Schur Form (T1), and Schur Vectors (UZ) 00739 * 00740 CALL ZLACPY( ' ', N, N, H, LDA, T1, LDA ) 00741 CALL ZLACPY( ' ', N, N, U, LDU, UZ, LDU ) 00742 * 00743 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, W1, UZ, LDU, 00744 $ WORK, NWORK, IINFO ) 00745 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 00746 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(V)', IINFO, N, JTYPE, 00747 $ IOLDSD 00748 INFO = ABS( IINFO ) 00749 GO TO 240 00750 END IF 00751 * 00752 * Compute Z = U' UZ 00753 * 00754 CALL ZGEMM( 'C', 'N', N, N, N, CONE, U, LDU, UZ, LDU, CZERO, 00755 $ Z, LDU ) 00756 NTEST = 8 00757 * 00758 * Do Tests 3: | H - Z T Z' | / ( |H| n ulp ) 00759 * and 4: | I - Z Z' | / ( n ulp ) 00760 * 00761 CALL ZHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK, 00762 $ NWORK, RWORK, RESULT( 3 ) ) 00763 * 00764 * Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp ) 00765 * and 6: | I - UZ (UZ)' | / ( n ulp ) 00766 * 00767 CALL ZHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK, 00768 $ NWORK, RWORK, RESULT( 5 ) ) 00769 * 00770 * Do Test 7: | T2 - T1 | / ( |T| n ulp ) 00771 * 00772 CALL ZGET10( N, N, T2, LDA, T1, LDA, WORK, RWORK, 00773 $ RESULT( 7 ) ) 00774 * 00775 * Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp ) 00776 * 00777 TEMP1 = ZERO 00778 TEMP2 = ZERO 00779 DO 130 J = 1, N 00780 TEMP1 = MAX( TEMP1, ABS( W1( J ) ), ABS( W3( J ) ) ) 00781 TEMP2 = MAX( TEMP2, ABS( W1( J )-W3( J ) ) ) 00782 130 CONTINUE 00783 * 00784 RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 00785 * 00786 * Compute the Left and Right Eigenvectors of T 00787 * 00788 * Compute the Right eigenvector Matrix: 00789 * 00790 NTEST = 9 00791 RESULT( 9 ) = ULPINV 00792 * 00793 * Select every other eigenvector 00794 * 00795 DO 140 J = 1, N 00796 SELECT( J ) = .FALSE. 00797 140 CONTINUE 00798 DO 150 J = 1, N, 2 00799 SELECT( J ) = .TRUE. 00800 150 CONTINUE 00801 CALL ZTREVC( 'Right', 'All', SELECT, N, T1, LDA, CDUMMA, 00802 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 00803 IF( IINFO.NE.0 ) THEN 00804 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,A)', IINFO, N, 00805 $ JTYPE, IOLDSD 00806 INFO = ABS( IINFO ) 00807 GO TO 240 00808 END IF 00809 * 00810 * Test 9: | TR - RW | / ( |T| |R| ulp ) 00811 * 00812 CALL ZGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, W1, 00813 $ WORK, RWORK, DUMMA( 1 ) ) 00814 RESULT( 9 ) = DUMMA( 1 ) 00815 IF( DUMMA( 2 ).GT.THRESH ) THEN 00816 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTREVC', 00817 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00818 END IF 00819 * 00820 * Compute selected right eigenvectors and confirm that 00821 * they agree with previous right eigenvectors 00822 * 00823 CALL ZTREVC( 'Right', 'Some', SELECT, N, T1, LDA, CDUMMA, 00824 $ LDU, EVECTL, LDU, N, IN, WORK, RWORK, IINFO ) 00825 IF( IINFO.NE.0 ) THEN 00826 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,S)', IINFO, N, 00827 $ JTYPE, IOLDSD 00828 INFO = ABS( IINFO ) 00829 GO TO 240 00830 END IF 00831 * 00832 K = 1 00833 MATCH = .TRUE. 00834 DO 170 J = 1, N 00835 IF( SELECT( J ) ) THEN 00836 DO 160 JJ = 1, N 00837 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN 00838 MATCH = .FALSE. 00839 GO TO 180 00840 END IF 00841 160 CONTINUE 00842 K = K + 1 00843 END IF 00844 170 CONTINUE 00845 180 CONTINUE 00846 IF( .NOT.MATCH ) 00847 $ WRITE( NOUNIT, FMT = 9997 )'Right', 'ZTREVC', N, JTYPE, 00848 $ IOLDSD 00849 * 00850 * Compute the Left eigenvector Matrix: 00851 * 00852 NTEST = 10 00853 RESULT( 10 ) = ULPINV 00854 CALL ZTREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU, 00855 $ CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 00856 IF( IINFO.NE.0 ) THEN 00857 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,A)', IINFO, N, 00858 $ JTYPE, IOLDSD 00859 INFO = ABS( IINFO ) 00860 GO TO 240 00861 END IF 00862 * 00863 * Test 10: | LT - WL | / ( |T| |L| ulp ) 00864 * 00865 CALL ZGET22( 'C', 'N', 'C', N, T1, LDA, EVECTL, LDU, W1, 00866 $ WORK, RWORK, DUMMA( 3 ) ) 00867 RESULT( 10 ) = DUMMA( 3 ) 00868 IF( DUMMA( 4 ).GT.THRESH ) THEN 00869 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTREVC', DUMMA( 4 ), 00870 $ N, JTYPE, IOLDSD 00871 END IF 00872 * 00873 * Compute selected left eigenvectors and confirm that 00874 * they agree with previous left eigenvectors 00875 * 00876 CALL ZTREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR, 00877 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 00878 IF( IINFO.NE.0 ) THEN 00879 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,S)', IINFO, N, 00880 $ JTYPE, IOLDSD 00881 INFO = ABS( IINFO ) 00882 GO TO 240 00883 END IF 00884 * 00885 K = 1 00886 MATCH = .TRUE. 00887 DO 200 J = 1, N 00888 IF( SELECT( J ) ) THEN 00889 DO 190 JJ = 1, N 00890 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN 00891 MATCH = .FALSE. 00892 GO TO 210 00893 END IF 00894 190 CONTINUE 00895 K = K + 1 00896 END IF 00897 200 CONTINUE 00898 210 CONTINUE 00899 IF( .NOT.MATCH ) 00900 $ WRITE( NOUNIT, FMT = 9997 )'Left', 'ZTREVC', N, JTYPE, 00901 $ IOLDSD 00902 * 00903 * Call ZHSEIN for Right eigenvectors of H, do test 11 00904 * 00905 NTEST = 11 00906 RESULT( 11 ) = ULPINV 00907 DO 220 J = 1, N 00908 SELECT( J ) = .TRUE. 00909 220 CONTINUE 00910 * 00911 CALL ZHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3, 00912 $ CDUMMA, LDU, EVECTX, LDU, N1, IN, WORK, RWORK, 00913 $ IWORK, IWORK, IINFO ) 00914 IF( IINFO.NE.0 ) THEN 00915 WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(R)', IINFO, N, JTYPE, 00916 $ IOLDSD 00917 INFO = ABS( IINFO ) 00918 IF( IINFO.LT.0 ) 00919 $ GO TO 240 00920 ELSE 00921 * 00922 * Test 11: | HX - XW | / ( |H| |X| ulp ) 00923 * 00924 * (from inverse iteration) 00925 * 00926 CALL ZGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, W3, 00927 $ WORK, RWORK, DUMMA( 1 ) ) 00928 IF( DUMMA( 1 ).LT.ULPINV ) 00929 $ RESULT( 11 ) = DUMMA( 1 )*ANINV 00930 IF( DUMMA( 2 ).GT.THRESH ) THEN 00931 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZHSEIN', 00932 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00933 END IF 00934 END IF 00935 * 00936 * Call ZHSEIN for Left eigenvectors of H, do test 12 00937 * 00938 NTEST = 12 00939 RESULT( 12 ) = ULPINV 00940 DO 230 J = 1, N 00941 SELECT( J ) = .TRUE. 00942 230 CONTINUE 00943 * 00944 CALL ZHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3, 00945 $ EVECTY, LDU, CDUMMA, LDU, N1, IN, WORK, RWORK, 00946 $ IWORK, IWORK, IINFO ) 00947 IF( IINFO.NE.0 ) THEN 00948 WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(L)', IINFO, N, JTYPE, 00949 $ IOLDSD 00950 INFO = ABS( IINFO ) 00951 IF( IINFO.LT.0 ) 00952 $ GO TO 240 00953 ELSE 00954 * 00955 * Test 12: | YH - WY | / ( |H| |Y| ulp ) 00956 * 00957 * (from inverse iteration) 00958 * 00959 CALL ZGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, W3, 00960 $ WORK, RWORK, DUMMA( 3 ) ) 00961 IF( DUMMA( 3 ).LT.ULPINV ) 00962 $ RESULT( 12 ) = DUMMA( 3 )*ANINV 00963 IF( DUMMA( 4 ).GT.THRESH ) THEN 00964 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZHSEIN', 00965 $ DUMMA( 4 ), N, JTYPE, IOLDSD 00966 END IF 00967 END IF 00968 * 00969 * Call ZUNMHR for Right eigenvectors of A, do test 13 00970 * 00971 NTEST = 13 00972 RESULT( 13 ) = ULPINV 00973 * 00974 CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 00975 $ LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO ) 00976 IF( IINFO.NE.0 ) THEN 00977 WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE, 00978 $ IOLDSD 00979 INFO = ABS( IINFO ) 00980 IF( IINFO.LT.0 ) 00981 $ GO TO 240 00982 ELSE 00983 * 00984 * Test 13: | AX - XW | / ( |A| |X| ulp ) 00985 * 00986 * (from inverse iteration) 00987 * 00988 CALL ZGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, W3, 00989 $ WORK, RWORK, DUMMA( 1 ) ) 00990 IF( DUMMA( 1 ).LT.ULPINV ) 00991 $ RESULT( 13 ) = DUMMA( 1 )*ANINV 00992 END IF 00993 * 00994 * Call ZUNMHR for Left eigenvectors of A, do test 14 00995 * 00996 NTEST = 14 00997 RESULT( 14 ) = ULPINV 00998 * 00999 CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 01000 $ LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO ) 01001 IF( IINFO.NE.0 ) THEN 01002 WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE, 01003 $ IOLDSD 01004 INFO = ABS( IINFO ) 01005 IF( IINFO.LT.0 ) 01006 $ GO TO 240 01007 ELSE 01008 * 01009 * Test 14: | YA - WY | / ( |A| |Y| ulp ) 01010 * 01011 * (from inverse iteration) 01012 * 01013 CALL ZGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, W3, 01014 $ WORK, RWORK, DUMMA( 3 ) ) 01015 IF( DUMMA( 3 ).LT.ULPINV ) 01016 $ RESULT( 14 ) = DUMMA( 3 )*ANINV 01017 END IF 01018 * 01019 * End of Loop -- Check for RESULT(j) > THRESH 01020 * 01021 240 CONTINUE 01022 * 01023 NTESTT = NTESTT + NTEST 01024 CALL DLAFTS( 'ZHS', N, N, JTYPE, NTEST, RESULT, IOLDSD, 01025 $ THRESH, NOUNIT, NERRS ) 01026 * 01027 250 CONTINUE 01028 260 CONTINUE 01029 * 01030 * Summary 01031 * 01032 CALL DLASUM( 'ZHS', NOUNIT, NERRS, NTESTT ) 01033 * 01034 RETURN 01035 * 01036 9999 FORMAT( ' ZCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 01037 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 01038 9998 FORMAT( ' ZCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ', 01039 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 01040 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 01041 $ ')' ) 01042 9997 FORMAT( ' ZCHKHS: Selected ', A, ' Eigenvectors from ', A, 01043 $ ' do not match other eigenvectors ', 9X, 'N=', I6, 01044 $ ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 01045 * 01046 * End of ZCHKHS 01047 * 01048 END
1.7.3 
