|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1, 00003 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1, 00004 $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK, 00005 $ RWORK, LLWORK, RESULT, INFO ) 00006 * 00007 * -- LAPACK test routine (version 3.1) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * November 2006 00010 * 00011 * .. Scalar Arguments .. 00012 LOGICAL TSTDIF 00013 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES 00014 REAL THRESH, THRSHN 00015 * .. 00016 * .. Array Arguments .. 00017 LOGICAL DOTYPE( * ), LLWORK( * ) 00018 INTEGER ISEED( 4 ), NN( * ) 00019 REAL RESULT( 15 ), RWORK( * ) 00020 COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ), 00021 $ B( LDA, * ), BETA1( * ), BETA3( * ), 00022 $ EVECTL( LDU, * ), EVECTR( LDU, * ), 00023 $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ), 00024 $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ), 00025 $ T( LDA, * ), U( LDU, * ), V( LDU, * ), 00026 $ WORK( * ), Z( LDU, * ) 00027 * .. 00028 * 00029 * Purpose 00030 * ======= 00031 * 00032 * CCHKGG checks the nonsymmetric generalized eigenvalue problem 00033 * routines. 00034 * H H H 00035 * CGGHRD factors A and B as U H V and U T V , where means conjugate 00036 * transpose, H is hessenberg, T is triangular and U and V are unitary. 00037 * 00038 * H H 00039 * CHGEQZ factors H and T as Q S Z and Q P Z , where P and S are upper 00040 * triangular and Q and Z are unitary. It also computes the generalized 00041 * eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)), where 00042 * alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus, w(j) = alpha(j)/beta(j) 00043 * is a root of the generalized eigenvalue problem 00044 * 00045 * det( A - w(j) B ) = 0 00046 * 00047 * and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent 00048 * problem 00049 * 00050 * det( m(j) A - B ) = 0 00051 * 00052 * CTGEVC computes the matrix L of left eigenvectors and the matrix R 00053 * of right eigenvectors for the matrix pair ( S, P ). In the 00054 * description below, l and r are left and right eigenvectors 00055 * corresponding to the generalized eigenvalues (alpha,beta). 00056 * 00057 * When CCHKGG is called, a number of matrix "sizes" ("n's") and a 00058 * number of matrix "types" are specified. For each size ("n") 00059 * and each type of matrix, one matrix will be generated and used 00060 * to test the nonsymmetric eigenroutines. For each matrix, 13 00061 * tests will be performed. The first twelve "test ratios" should be 00062 * small -- O(1). They will be compared with the threshhold THRESH: 00063 * 00064 * H 00065 * (1) | A - U H V | / ( |A| n ulp ) 00066 * 00067 * H 00068 * (2) | B - U T V | / ( |B| n ulp ) 00069 * 00070 * H 00071 * (3) | I - UU | / ( n ulp ) 00072 * 00073 * H 00074 * (4) | I - VV | / ( n ulp ) 00075 * 00076 * H 00077 * (5) | H - Q S Z | / ( |H| n ulp ) 00078 * 00079 * H 00080 * (6) | T - Q P Z | / ( |T| n ulp ) 00081 * 00082 * H 00083 * (7) | I - QQ | / ( n ulp ) 00084 * 00085 * H 00086 * (8) | I - ZZ | / ( n ulp ) 00087 * 00088 * (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of 00089 * H 00090 * | (beta A - alpha B) l | / ( ulp max( |beta A|, |alpha B| ) ) 00091 * 00092 * (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of 00093 * H 00094 * | (beta H - alpha T) l' | / ( ulp max( |beta H|, |alpha T| ) ) 00095 * 00096 * where the eigenvectors l' are the result of passing Q to 00097 * STGEVC and back transforming (JOB='B'). 00098 * 00099 * (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of 00100 * 00101 * | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) ) 00102 * 00103 * (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of 00104 * 00105 * | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) ) 00106 * 00107 * where the eigenvectors r' are the result of passing Z to 00108 * STGEVC and back transforming (JOB='B'). 00109 * 00110 * The last three test ratios will usually be small, but there is no 00111 * mathematical requirement that they be so. They are therefore 00112 * compared with THRESH only if TSTDIF is .TRUE. 00113 * 00114 * (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp ) 00115 * 00116 * (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp ) 00117 * 00118 * (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| , 00119 * |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp 00120 * 00121 * In addition, the normalization of L and R are checked, and compared 00122 * with the threshhold THRSHN. 00123 * 00124 * Test Matrices 00125 * ---- -------- 00126 * 00127 * The sizes of the test matrices are specified by an array 00128 * NN(1:NSIZES); the value of each element NN(j) specifies one size. 00129 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if 00130 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00131 * Currently, the list of possible types is: 00132 * 00133 * (1) ( 0, 0 ) (a pair of zero matrices) 00134 * 00135 * (2) ( I, 0 ) (an identity and a zero matrix) 00136 * 00137 * (3) ( 0, I ) (an identity and a zero matrix) 00138 * 00139 * (4) ( I, I ) (a pair of identity matrices) 00140 * 00141 * t t 00142 * (5) ( J , J ) (a pair of transposed Jordan blocks) 00143 * 00144 * t ( I 0 ) 00145 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t ) 00146 * ( 0 I ) ( 0 J ) 00147 * and I is a k x k identity and J a (k+1)x(k+1) 00148 * Jordan block; k=(N-1)/2 00149 * 00150 * (7) ( D, I ) where D is P*D1, P is a random unitary diagonal 00151 * matrix (i.e., with random magnitude 1 entries 00152 * on the diagonal), and D1=diag( 0, 1,..., N-1 ) 00153 * (i.e., a diagonal matrix with D1(1,1)=0, 00154 * D1(2,2)=1, ..., D1(N,N)=N-1.) 00155 * (8) ( I, D ) 00156 * 00157 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big 00158 * 00159 * (10) ( small*D, big*I ) 00160 * 00161 * (11) ( big*I, small*D ) 00162 * 00163 * (12) ( small*I, big*D ) 00164 * 00165 * (13) ( big*D, big*I ) 00166 * 00167 * (14) ( small*D, small*I ) 00168 * 00169 * (15) ( D1, D2 ) where D1=P*diag( 0, 0, 1, ..., N-3, 0 ) and 00170 * D2=Q*diag( 0, N-3, N-4,..., 1, 0, 0 ), and 00171 * P and Q are random unitary diagonal matrices. 00172 * t t 00173 * (16) U ( J , J ) V where U and V are random unitary matrices. 00174 * 00175 * (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices 00176 * with random O(1) entries above the diagonal 00177 * and diagonal entries diag(T1) = 00178 * P*( 0, 0, 1, ..., N-3, 0 ) and diag(T2) = 00179 * Q*( 0, N-3, N-4,..., 1, 0, 0 ) 00180 * 00181 * (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 ) 00182 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 ) 00183 * s = machine precision. 00184 * 00185 * (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 ) 00186 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 ) 00187 * 00188 * N-5 00189 * (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 ) 00190 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00191 * 00192 * (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 ) 00193 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00194 * where r1,..., r(N-4) are random. 00195 * 00196 * (22) U ( big*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00197 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00198 * 00199 * (23) U ( small*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00200 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00201 * 00202 * (24) U ( small*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00203 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00204 * 00205 * (25) U ( big*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00206 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00207 * 00208 * (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular 00209 * matrices. 00210 * 00211 * Arguments 00212 * ========= 00213 * 00214 * NSIZES (input) INTEGER 00215 * The number of sizes of matrices to use. If it is zero, 00216 * CCHKGG does nothing. It must be at least zero. 00217 * 00218 * NN (input) INTEGER array, dimension (NSIZES) 00219 * An array containing the sizes to be used for the matrices. 00220 * Zero values will be skipped. The values must be at least 00221 * zero. 00222 * 00223 * NTYPES (input) INTEGER 00224 * The number of elements in DOTYPE. If it is zero, CCHKGG 00225 * does nothing. It must be at least zero. If it is MAXTYP+1 00226 * and NSIZES is 1, then an additional type, MAXTYP+1 is 00227 * defined, which is to use whatever matrix is in A. This 00228 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00229 * DOTYPE(MAXTYP+1) is .TRUE. . 00230 * 00231 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00232 * If DOTYPE(j) is .TRUE., then for each size in NN a 00233 * matrix of that size and of type j will be generated. 00234 * If NTYPES is smaller than the maximum number of types 00235 * defined (PARAMETER MAXTYP), then types NTYPES+1 through 00236 * MAXTYP will not be generated. If NTYPES is larger 00237 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 00238 * will be ignored. 00239 * 00240 * ISEED (input/output) INTEGER array, dimension (4) 00241 * On entry ISEED specifies the seed of the random number 00242 * generator. The array elements should be between 0 and 4095; 00243 * if not they will be reduced mod 4096. Also, ISEED(4) must 00244 * be odd. The random number generator uses a linear 00245 * congruential sequence limited to small integers, and so 00246 * should produce machine independent random numbers. The 00247 * values of ISEED are changed on exit, and can be used in the 00248 * next call to CCHKGG to continue the same random number 00249 * sequence. 00250 * 00251 * THRESH (input) REAL 00252 * A test will count as "failed" if the "error", computed as 00253 * described above, exceeds THRESH. Note that the error 00254 * is scaled to be O(1), so THRESH should be a reasonably 00255 * small multiple of 1, e.g., 10 or 100. In particular, 00256 * it should not depend on the precision (single vs. double) 00257 * or the size of the matrix. It must be at least zero. 00258 * 00259 * TSTDIF (input) LOGICAL 00260 * Specifies whether test ratios 13-15 will be computed and 00261 * compared with THRESH. 00262 * = .FALSE.: Only test ratios 1-12 will be computed and tested. 00263 * Ratios 13-15 will be set to zero. 00264 * = .TRUE.: All the test ratios 1-15 will be computed and 00265 * tested. 00266 * 00267 * THRSHN (input) REAL 00268 * Threshhold for reporting eigenvector normalization error. 00269 * If the normalization of any eigenvector differs from 1 by 00270 * more than THRSHN*ulp, then a special error message will be 00271 * printed. (This is handled separately from the other tests, 00272 * since only a compiler or programming error should cause an 00273 * error message, at least if THRSHN is at least 5--10.) 00274 * 00275 * NOUNIT (input) INTEGER 00276 * The FORTRAN unit number for printing out error messages 00277 * (e.g., if a routine returns IINFO not equal to 0.) 00278 * 00279 * A (input/workspace) COMPLEX array, dimension (LDA, max(NN)) 00280 * Used to hold the original A matrix. Used as input only 00281 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00282 * DOTYPE(MAXTYP+1)=.TRUE. 00283 * 00284 * LDA (input) INTEGER 00285 * The leading dimension of A, B, H, T, S1, P1, S2, and P2. 00286 * It must be at least 1 and at least max( NN ). 00287 * 00288 * B (input/workspace) COMPLEX array, dimension (LDA, max(NN)) 00289 * Used to hold the original B matrix. Used as input only 00290 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00291 * DOTYPE(MAXTYP+1)=.TRUE. 00292 * 00293 * H (workspace) COMPLEX array, dimension (LDA, max(NN)) 00294 * The upper Hessenberg matrix computed from A by CGGHRD. 00295 * 00296 * T (workspace) COMPLEX array, dimension (LDA, max(NN)) 00297 * The upper triangular matrix computed from B by CGGHRD. 00298 * 00299 * S1 (workspace) COMPLEX array, dimension (LDA, max(NN)) 00300 * The Schur (upper triangular) matrix computed from H by CHGEQZ 00301 * when Q and Z are also computed. 00302 * 00303 * S2 (workspace) COMPLEX array, dimension (LDA, max(NN)) 00304 * The Schur (upper triangular) matrix computed from H by CHGEQZ 00305 * when Q and Z are not computed. 00306 * 00307 * P1 (workspace) COMPLEX array, dimension (LDA, max(NN)) 00308 * The upper triangular matrix computed from T by CHGEQZ 00309 * when Q and Z are also computed. 00310 * 00311 * P2 (workspace) COMPLEX array, dimension (LDA, max(NN)) 00312 * The upper triangular matrix computed from T by CHGEQZ 00313 * when Q and Z are not computed. 00314 * 00315 * U (workspace) COMPLEX array, dimension (LDU, max(NN)) 00316 * The (left) unitary matrix computed by CGGHRD. 00317 * 00318 * LDU (input) INTEGER 00319 * The leading dimension of U, V, Q, Z, EVECTL, and EVECTR. It 00320 * must be at least 1 and at least max( NN ). 00321 * 00322 * V (workspace) COMPLEX array, dimension (LDU, max(NN)) 00323 * The (right) unitary matrix computed by CGGHRD. 00324 * 00325 * Q (workspace) COMPLEX array, dimension (LDU, max(NN)) 00326 * The (left) unitary matrix computed by CHGEQZ. 00327 * 00328 * Z (workspace) COMPLEX array, dimension (LDU, max(NN)) 00329 * The (left) unitary matrix computed by CHGEQZ. 00330 * 00331 * ALPHA1 (workspace) COMPLEX array, dimension (max(NN)) 00332 * BETA1 (workspace) COMPLEX array, dimension (max(NN)) 00333 * The generalized eigenvalues of (A,B) computed by CHGEQZ 00334 * when Q, Z, and the full Schur matrices are computed. 00335 * 00336 * ALPHA3 (workspace) COMPLEX array, dimension (max(NN)) 00337 * BETA3 (workspace) COMPLEX array, dimension (max(NN)) 00338 * The generalized eigenvalues of (A,B) computed by CHGEQZ 00339 * when neither Q, Z, nor the Schur matrices are computed. 00340 * 00341 * EVECTL (workspace) COMPLEX array, dimension (LDU, max(NN)) 00342 * The (lower triangular) left eigenvector matrix for the 00343 * matrices in S1 and P1. 00344 * 00345 * EVECTR (workspace) COMPLEX array, dimension (LDU, max(NN)) 00346 * The (upper triangular) right eigenvector matrix for the 00347 * matrices in S1 and P1. 00348 * 00349 * WORK (workspace) COMPLEX array, dimension (LWORK) 00350 * 00351 * LWORK (input) INTEGER 00352 * The number of entries in WORK. This must be at least 00353 * max( 4*N, 2 * N**2, 1 ), for all N=NN(j). 00354 * 00355 * RWORK (workspace) REAL array, dimension (2*max(NN)) 00356 * 00357 * LLWORK (workspace) LOGICAL array, dimension (max(NN)) 00358 * 00359 * RESULT (output) REAL array, dimension (15) 00360 * The values computed by the tests described above. 00361 * The values are currently limited to 1/ulp, to avoid 00362 * overflow. 00363 * 00364 * INFO (output) INTEGER 00365 * = 0: successful exit. 00366 * < 0: if INFO = -i, the i-th argument had an illegal value. 00367 * > 0: A routine returned an error code. INFO is the 00368 * absolute value of the INFO value returned. 00369 * 00370 * ===================================================================== 00371 * 00372 * .. Parameters .. 00373 REAL ZERO, ONE 00374 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00375 COMPLEX CZERO, CONE 00376 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00377 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00378 INTEGER MAXTYP 00379 PARAMETER ( MAXTYP = 26 ) 00380 * .. 00381 * .. Local Scalars .. 00382 LOGICAL BADNN 00383 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE, 00384 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX, 00385 $ NTEST, NTESTT 00386 REAL ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2, 00387 $ ULP, ULPINV 00388 COMPLEX CTEMP 00389 * .. 00390 * .. Local Arrays .. 00391 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP ) 00392 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 00393 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 00394 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 00395 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 00396 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 00397 REAL DUMMA( 4 ), RMAGN( 0: 3 ) 00398 COMPLEX CDUMMA( 4 ) 00399 * .. 00400 * .. External Functions .. 00401 REAL CLANGE, SLAMCH 00402 COMPLEX CLARND 00403 EXTERNAL CLANGE, SLAMCH, CLARND 00404 * .. 00405 * .. External Subroutines .. 00406 EXTERNAL CGEQR2, CGET51, CGET52, CGGHRD, CHGEQZ, CLACPY, 00407 $ CLARFG, CLASET, CLATM4, CTGEVC, CUNM2R, SLABAD, 00408 $ SLASUM, XERBLA 00409 * .. 00410 * .. Intrinsic Functions .. 00411 INTRINSIC ABS, CONJG, MAX, MIN, REAL, SIGN 00412 * .. 00413 * .. Data statements .. 00414 DATA KCLASS / 15*1, 10*2, 1*3 / 00415 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 00416 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 00417 DATA KADD / 0, 0, 0, 0, 3, 2 / 00418 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 00419 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 00420 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 00421 $ 1, 1, -4, 2, -4, 8*8, 0 / 00422 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 00423 $ 4*5, 4*3, 1 / 00424 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 00425 $ 4*6, 4*4, 1 / 00426 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 00427 $ 2, 1 / 00428 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 00429 $ 2, 1 / 00430 DATA KTRIAN / 16*0, 10*1 / 00431 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE., 00432 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE., 00433 $ 3*.FALSE., 5*.TRUE., .FALSE. / 00434 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE., 00435 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE., 00436 $ 9*.FALSE. / 00437 * .. 00438 * .. Executable Statements .. 00439 * 00440 * Check for errors 00441 * 00442 INFO = 0 00443 * 00444 BADNN = .FALSE. 00445 NMAX = 1 00446 DO 10 J = 1, NSIZES 00447 NMAX = MAX( NMAX, NN( J ) ) 00448 IF( NN( J ).LT.0 ) 00449 $ BADNN = .TRUE. 00450 10 CONTINUE 00451 * 00452 LWKOPT = MAX( 2*NMAX*NMAX, 4*NMAX, 1 ) 00453 * 00454 * Check for errors 00455 * 00456 IF( NSIZES.LT.0 ) THEN 00457 INFO = -1 00458 ELSE IF( BADNN ) THEN 00459 INFO = -2 00460 ELSE IF( NTYPES.LT.0 ) THEN 00461 INFO = -3 00462 ELSE IF( THRESH.LT.ZERO ) THEN 00463 INFO = -6 00464 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 00465 INFO = -10 00466 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN 00467 INFO = -19 00468 ELSE IF( LWKOPT.GT.LWORK ) THEN 00469 INFO = -30 00470 END IF 00471 * 00472 IF( INFO.NE.0 ) THEN 00473 CALL XERBLA( 'CCHKGG', -INFO ) 00474 RETURN 00475 END IF 00476 * 00477 * Quick return if possible 00478 * 00479 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00480 $ RETURN 00481 * 00482 SAFMIN = SLAMCH( 'Safe minimum' ) 00483 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00484 SAFMIN = SAFMIN / ULP 00485 SAFMAX = ONE / SAFMIN 00486 CALL SLABAD( SAFMIN, SAFMAX ) 00487 ULPINV = ONE / ULP 00488 * 00489 * The values RMAGN(2:3) depend on N, see below. 00490 * 00491 RMAGN( 0 ) = ZERO 00492 RMAGN( 1 ) = ONE 00493 * 00494 * Loop over sizes, types 00495 * 00496 NTESTT = 0 00497 NERRS = 0 00498 NMATS = 0 00499 * 00500 DO 240 JSIZE = 1, NSIZES 00501 N = NN( JSIZE ) 00502 N1 = MAX( 1, N ) 00503 RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 ) 00504 RMAGN( 3 ) = SAFMIN*ULPINV*N1 00505 * 00506 IF( NSIZES.NE.1 ) THEN 00507 MTYPES = MIN( MAXTYP, NTYPES ) 00508 ELSE 00509 MTYPES = MIN( MAXTYP+1, NTYPES ) 00510 END IF 00511 * 00512 DO 230 JTYPE = 1, MTYPES 00513 IF( .NOT.DOTYPE( JTYPE ) ) 00514 $ GO TO 230 00515 NMATS = NMATS + 1 00516 NTEST = 0 00517 * 00518 * Save ISEED in case of an error. 00519 * 00520 DO 20 J = 1, 4 00521 IOLDSD( J ) = ISEED( J ) 00522 20 CONTINUE 00523 * 00524 * Initialize RESULT 00525 * 00526 DO 30 J = 1, 15 00527 RESULT( J ) = ZERO 00528 30 CONTINUE 00529 * 00530 * Compute A and B 00531 * 00532 * Description of control parameters: 00533 * 00534 * KCLASS: =1 means w/o rotation, =2 means w/ rotation, 00535 * =3 means random. 00536 * KATYPE: the "type" to be passed to CLATM4 for computing A. 00537 * KAZERO: the pattern of zeros on the diagonal for A: 00538 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 00539 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 00540 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 00541 * non-zero entries.) 00542 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 00543 * =2: large, =3: small. 00544 * LASIGN: .TRUE. if the diagonal elements of A are to be 00545 * multiplied by a random magnitude 1 number. 00546 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. 00547 * KTRIAN: =0: don't fill in the upper triangle, =1: do. 00548 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 00549 * RMAGN: used to implement KAMAGN and KBMAGN. 00550 * 00551 IF( MTYPES.GT.MAXTYP ) 00552 $ GO TO 110 00553 IINFO = 0 00554 IF( KCLASS( JTYPE ).LT.3 ) THEN 00555 * 00556 * Generate A (w/o rotation) 00557 * 00558 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 00559 IN = 2*( ( N-1 ) / 2 ) + 1 00560 IF( IN.NE.N ) 00561 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 00562 ELSE 00563 IN = N 00564 END IF 00565 CALL CLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 00566 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 00567 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 00568 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4, 00569 $ ISEED, A, LDA ) 00570 IADD = KADD( KAZERO( JTYPE ) ) 00571 IF( IADD.GT.0 .AND. IADD.LE.N ) 00572 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 00573 * 00574 * Generate B (w/o rotation) 00575 * 00576 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 00577 IN = 2*( ( N-1 ) / 2 ) + 1 00578 IF( IN.NE.N ) 00579 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 00580 ELSE 00581 IN = N 00582 END IF 00583 CALL CLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 00584 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 00585 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 00586 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4, 00587 $ ISEED, B, LDA ) 00588 IADD = KADD( KBZERO( JTYPE ) ) 00589 IF( IADD.NE.0 ) 00590 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 00591 * 00592 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 00593 * 00594 * Include rotations 00595 * 00596 * Generate U, V as Householder transformations times a 00597 * diagonal matrix. (Note that CLARFG makes U(j,j) and 00598 * V(j,j) real.) 00599 * 00600 DO 50 JC = 1, N - 1 00601 DO 40 JR = JC, N 00602 U( JR, JC ) = CLARND( 3, ISEED ) 00603 V( JR, JC ) = CLARND( 3, ISEED ) 00604 40 CONTINUE 00605 CALL CLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1, 00606 $ WORK( JC ) ) 00607 WORK( 2*N+JC ) = SIGN( ONE, REAL( U( JC, JC ) ) ) 00608 U( JC, JC ) = CONE 00609 CALL CLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1, 00610 $ WORK( N+JC ) ) 00611 WORK( 3*N+JC ) = SIGN( ONE, REAL( V( JC, JC ) ) ) 00612 V( JC, JC ) = CONE 00613 50 CONTINUE 00614 CTEMP = CLARND( 3, ISEED ) 00615 U( N, N ) = CONE 00616 WORK( N ) = CZERO 00617 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 00618 CTEMP = CLARND( 3, ISEED ) 00619 V( N, N ) = CONE 00620 WORK( 2*N ) = CZERO 00621 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 00622 * 00623 * Apply the diagonal matrices 00624 * 00625 DO 70 JC = 1, N 00626 DO 60 JR = 1, N 00627 A( JR, JC ) = WORK( 2*N+JR )* 00628 $ CONJG( WORK( 3*N+JC ) )* 00629 $ A( JR, JC ) 00630 B( JR, JC ) = WORK( 2*N+JR )* 00631 $ CONJG( WORK( 3*N+JC ) )* 00632 $ B( JR, JC ) 00633 60 CONTINUE 00634 70 CONTINUE 00635 CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A, 00636 $ LDA, WORK( 2*N+1 ), IINFO ) 00637 IF( IINFO.NE.0 ) 00638 $ GO TO 100 00639 CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ), 00640 $ A, LDA, WORK( 2*N+1 ), IINFO ) 00641 IF( IINFO.NE.0 ) 00642 $ GO TO 100 00643 CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B, 00644 $ LDA, WORK( 2*N+1 ), IINFO ) 00645 IF( IINFO.NE.0 ) 00646 $ GO TO 100 00647 CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ), 00648 $ B, LDA, WORK( 2*N+1 ), IINFO ) 00649 IF( IINFO.NE.0 ) 00650 $ GO TO 100 00651 END IF 00652 ELSE 00653 * 00654 * Random matrices 00655 * 00656 DO 90 JC = 1, N 00657 DO 80 JR = 1, N 00658 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 00659 $ CLARND( 4, ISEED ) 00660 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 00661 $ CLARND( 4, ISEED ) 00662 80 CONTINUE 00663 90 CONTINUE 00664 END IF 00665 * 00666 ANORM = CLANGE( '1', N, N, A, LDA, RWORK ) 00667 BNORM = CLANGE( '1', N, N, B, LDA, RWORK ) 00668 * 00669 100 CONTINUE 00670 * 00671 IF( IINFO.NE.0 ) THEN 00672 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 00673 $ IOLDSD 00674 INFO = ABS( IINFO ) 00675 RETURN 00676 END IF 00677 * 00678 110 CONTINUE 00679 * 00680 * Call CGEQR2, CUNM2R, and CGGHRD to compute H, T, U, and V 00681 * 00682 CALL CLACPY( ' ', N, N, A, LDA, H, LDA ) 00683 CALL CLACPY( ' ', N, N, B, LDA, T, LDA ) 00684 NTEST = 1 00685 RESULT( 1 ) = ULPINV 00686 * 00687 CALL CGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO ) 00688 IF( IINFO.NE.0 ) THEN 00689 WRITE( NOUNIT, FMT = 9999 )'CGEQR2', IINFO, N, JTYPE, 00690 $ IOLDSD 00691 INFO = ABS( IINFO ) 00692 GO TO 210 00693 END IF 00694 * 00695 CALL CUNM2R( 'L', 'C', N, N, N, T, LDA, WORK, H, LDA, 00696 $ WORK( N+1 ), IINFO ) 00697 IF( IINFO.NE.0 ) THEN 00698 WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE, 00699 $ IOLDSD 00700 INFO = ABS( IINFO ) 00701 GO TO 210 00702 END IF 00703 * 00704 CALL CLASET( 'Full', N, N, CZERO, CONE, U, LDU ) 00705 CALL CUNM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU, 00706 $ WORK( N+1 ), IINFO ) 00707 IF( IINFO.NE.0 ) THEN 00708 WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE, 00709 $ IOLDSD 00710 INFO = ABS( IINFO ) 00711 GO TO 210 00712 END IF 00713 * 00714 CALL CGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V, 00715 $ LDU, IINFO ) 00716 IF( IINFO.NE.0 ) THEN 00717 WRITE( NOUNIT, FMT = 9999 )'CGGHRD', IINFO, N, JTYPE, 00718 $ IOLDSD 00719 INFO = ABS( IINFO ) 00720 GO TO 210 00721 END IF 00722 NTEST = 4 00723 * 00724 * Do tests 1--4 00725 * 00726 CALL CGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK, 00727 $ RWORK, RESULT( 1 ) ) 00728 CALL CGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK, 00729 $ RWORK, RESULT( 2 ) ) 00730 CALL CGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK, 00731 $ RWORK, RESULT( 3 ) ) 00732 CALL CGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK, 00733 $ RWORK, RESULT( 4 ) ) 00734 * 00735 * Call CHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests. 00736 * 00737 * Compute T1 and UZ 00738 * 00739 * Eigenvalues only 00740 * 00741 CALL CLACPY( ' ', N, N, H, LDA, S2, LDA ) 00742 CALL CLACPY( ' ', N, N, T, LDA, P2, LDA ) 00743 NTEST = 5 00744 RESULT( 5 ) = ULPINV 00745 * 00746 CALL CHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA, 00747 $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK, 00748 $ RWORK, IINFO ) 00749 IF( IINFO.NE.0 ) THEN 00750 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(E)', IINFO, N, JTYPE, 00751 $ IOLDSD 00752 INFO = ABS( IINFO ) 00753 GO TO 210 00754 END IF 00755 * 00756 * Eigenvalues and Full Schur Form 00757 * 00758 CALL CLACPY( ' ', N, N, H, LDA, S2, LDA ) 00759 CALL CLACPY( ' ', N, N, T, LDA, P2, LDA ) 00760 * 00761 CALL CHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA, 00762 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK, 00763 $ RWORK, IINFO ) 00764 IF( IINFO.NE.0 ) THEN 00765 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(S)', IINFO, N, JTYPE, 00766 $ IOLDSD 00767 INFO = ABS( IINFO ) 00768 GO TO 210 00769 END IF 00770 * 00771 * Eigenvalues, Schur Form, and Schur Vectors 00772 * 00773 CALL CLACPY( ' ', N, N, H, LDA, S1, LDA ) 00774 CALL CLACPY( ' ', N, N, T, LDA, P1, LDA ) 00775 * 00776 CALL CHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA, 00777 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK, 00778 $ RWORK, IINFO ) 00779 IF( IINFO.NE.0 ) THEN 00780 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(V)', IINFO, N, JTYPE, 00781 $ IOLDSD 00782 INFO = ABS( IINFO ) 00783 GO TO 210 00784 END IF 00785 * 00786 NTEST = 8 00787 * 00788 * Do Tests 5--8 00789 * 00790 CALL CGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK, 00791 $ RWORK, RESULT( 5 ) ) 00792 CALL CGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK, 00793 $ RWORK, RESULT( 6 ) ) 00794 CALL CGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK, 00795 $ RWORK, RESULT( 7 ) ) 00796 CALL CGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK, 00797 $ RWORK, RESULT( 8 ) ) 00798 * 00799 * Compute the Left and Right Eigenvectors of (S1,P1) 00800 * 00801 * 9: Compute the left eigenvector Matrix without 00802 * back transforming: 00803 * 00804 NTEST = 9 00805 RESULT( 9 ) = ULPINV 00806 * 00807 * To test "SELECT" option, compute half of the eigenvectors 00808 * in one call, and half in another 00809 * 00810 I1 = N / 2 00811 DO 120 J = 1, I1 00812 LLWORK( J ) = .TRUE. 00813 120 CONTINUE 00814 DO 130 J = I1 + 1, N 00815 LLWORK( J ) = .FALSE. 00816 130 CONTINUE 00817 * 00818 CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL, 00819 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 00820 IF( IINFO.NE.0 ) THEN 00821 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S1)', IINFO, N, 00822 $ JTYPE, IOLDSD 00823 INFO = ABS( IINFO ) 00824 GO TO 210 00825 END IF 00826 * 00827 I1 = IN 00828 DO 140 J = 1, I1 00829 LLWORK( J ) = .FALSE. 00830 140 CONTINUE 00831 DO 150 J = I1 + 1, N 00832 LLWORK( J ) = .TRUE. 00833 150 CONTINUE 00834 * 00835 CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, 00836 $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN, 00837 $ WORK, RWORK, IINFO ) 00838 IF( IINFO.NE.0 ) THEN 00839 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S2)', IINFO, N, 00840 $ JTYPE, IOLDSD 00841 INFO = ABS( IINFO ) 00842 GO TO 210 00843 END IF 00844 * 00845 CALL CGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU, 00846 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 00847 RESULT( 9 ) = DUMMA( 1 ) 00848 IF( DUMMA( 2 ).GT.THRSHN ) THEN 00849 WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=S)', 00850 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00851 END IF 00852 * 00853 * 10: Compute the left eigenvector Matrix with 00854 * back transforming: 00855 * 00856 NTEST = 10 00857 RESULT( 10 ) = ULPINV 00858 CALL CLACPY( 'F', N, N, Q, LDU, EVECTL, LDU ) 00859 CALL CTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL, 00860 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 00861 IF( IINFO.NE.0 ) THEN 00862 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,B)', IINFO, N, 00863 $ JTYPE, IOLDSD 00864 INFO = ABS( IINFO ) 00865 GO TO 210 00866 END IF 00867 * 00868 CALL CGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1, 00869 $ BETA1, WORK, RWORK, DUMMA( 1 ) ) 00870 RESULT( 10 ) = DUMMA( 1 ) 00871 IF( DUMMA( 2 ).GT.THRSHN ) THEN 00872 WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=B)', 00873 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00874 END IF 00875 * 00876 * 11: Compute the right eigenvector Matrix without 00877 * back transforming: 00878 * 00879 NTEST = 11 00880 RESULT( 11 ) = ULPINV 00881 * 00882 * To test "SELECT" option, compute half of the eigenvectors 00883 * in one call, and half in another 00884 * 00885 I1 = N / 2 00886 DO 160 J = 1, I1 00887 LLWORK( J ) = .TRUE. 00888 160 CONTINUE 00889 DO 170 J = I1 + 1, N 00890 LLWORK( J ) = .FALSE. 00891 170 CONTINUE 00892 * 00893 CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 00894 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 00895 IF( IINFO.NE.0 ) THEN 00896 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S1)', IINFO, N, 00897 $ JTYPE, IOLDSD 00898 INFO = ABS( IINFO ) 00899 GO TO 210 00900 END IF 00901 * 00902 I1 = IN 00903 DO 180 J = 1, I1 00904 LLWORK( J ) = .FALSE. 00905 180 CONTINUE 00906 DO 190 J = I1 + 1, N 00907 LLWORK( J ) = .TRUE. 00908 190 CONTINUE 00909 * 00910 CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 00911 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK, 00912 $ RWORK, IINFO ) 00913 IF( IINFO.NE.0 ) THEN 00914 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S2)', IINFO, N, 00915 $ JTYPE, IOLDSD 00916 INFO = ABS( IINFO ) 00917 GO TO 210 00918 END IF 00919 * 00920 CALL CGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU, 00921 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 00922 RESULT( 11 ) = DUMMA( 1 ) 00923 IF( DUMMA( 2 ).GT.THRESH ) THEN 00924 WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=S)', 00925 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00926 END IF 00927 * 00928 * 12: Compute the right eigenvector Matrix with 00929 * back transforming: 00930 * 00931 NTEST = 12 00932 RESULT( 12 ) = ULPINV 00933 CALL CLACPY( 'F', N, N, Z, LDU, EVECTR, LDU ) 00934 CALL CTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 00935 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 00936 IF( IINFO.NE.0 ) THEN 00937 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,B)', IINFO, N, 00938 $ JTYPE, IOLDSD 00939 INFO = ABS( IINFO ) 00940 GO TO 210 00941 END IF 00942 * 00943 CALL CGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU, 00944 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 00945 RESULT( 12 ) = DUMMA( 1 ) 00946 IF( DUMMA( 2 ).GT.THRESH ) THEN 00947 WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=B)', 00948 $ DUMMA( 2 ), N, JTYPE, IOLDSD 00949 END IF 00950 * 00951 * Tests 13--15 are done only on request 00952 * 00953 IF( TSTDIF ) THEN 00954 * 00955 * Do Tests 13--14 00956 * 00957 CALL CGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU, 00958 $ WORK, RWORK, RESULT( 13 ) ) 00959 CALL CGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU, 00960 $ WORK, RWORK, RESULT( 14 ) ) 00961 * 00962 * Do Test 15 00963 * 00964 TEMP1 = ZERO 00965 TEMP2 = ZERO 00966 DO 200 J = 1, N 00967 TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) ) 00968 TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) ) 00969 200 CONTINUE 00970 * 00971 TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) ) 00972 TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) ) 00973 RESULT( 15 ) = MAX( TEMP1, TEMP2 ) 00974 NTEST = 15 00975 ELSE 00976 RESULT( 13 ) = ZERO 00977 RESULT( 14 ) = ZERO 00978 RESULT( 15 ) = ZERO 00979 NTEST = 12 00980 END IF 00981 * 00982 * End of Loop -- Check for RESULT(j) > THRESH 00983 * 00984 210 CONTINUE 00985 * 00986 NTESTT = NTESTT + NTEST 00987 * 00988 * Print out tests which fail. 00989 * 00990 DO 220 JR = 1, NTEST 00991 IF( RESULT( JR ).GE.THRESH ) THEN 00992 * 00993 * If this is the first test to fail, 00994 * print a header to the data file. 00995 * 00996 IF( NERRS.EQ.0 ) THEN 00997 WRITE( NOUNIT, FMT = 9997 )'CGG' 00998 * 00999 * Matrix types 01000 * 01001 WRITE( NOUNIT, FMT = 9996 ) 01002 WRITE( NOUNIT, FMT = 9995 ) 01003 WRITE( NOUNIT, FMT = 9994 )'Unitary' 01004 * 01005 * Tests performed 01006 * 01007 WRITE( NOUNIT, FMT = 9993 )'unitary', '*', 01008 $ 'conjugate transpose', ( '*', J = 1, 10 ) 01009 * 01010 END IF 01011 NERRS = NERRS + 1 01012 IF( RESULT( JR ).LT.10000.0 ) THEN 01013 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 01014 $ RESULT( JR ) 01015 ELSE 01016 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 01017 $ RESULT( JR ) 01018 END IF 01019 END IF 01020 220 CONTINUE 01021 * 01022 230 CONTINUE 01023 240 CONTINUE 01024 * 01025 * Summary 01026 * 01027 CALL SLASUM( 'CGG', NOUNIT, NERRS, NTESTT ) 01028 RETURN 01029 * 01030 9999 FORMAT( ' CCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 01031 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 01032 * 01033 9998 FORMAT( ' CCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ', 01034 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 01035 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 01036 $ ')' ) 01037 * 01038 9997 FORMAT( 1X, A3, ' -- Complex Generalized eigenvalue problem' ) 01039 * 01040 9996 FORMAT( ' Matrix types (see CCHKGG for details): ' ) 01041 * 01042 9995 FORMAT( ' Special Matrices:', 23X, 01043 $ '(J''=transposed Jordan block)', 01044 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 01045 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 01046 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 01047 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 01048 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 01049 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 01050 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 01051 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 01052 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 01053 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 01054 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 01055 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 01056 $ '23=(small,large) 24=(small,small) 25=(large,large)', 01057 $ / ' 26=random O(1) matrices.' ) 01058 * 01059 9993 FORMAT( / ' Tests performed: (H is Hessenberg, S is Schur, B, ', 01060 $ 'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A, 01061 $ ', l and r are the', / 20X, 01062 $ 'appropriate left and right eigenvectors, resp., a is', 01063 $ / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)', 01064 $ / ' 1 = | A - U H V', A, 01065 $ ' | / ( |A| n ulp ) 2 = | B - U T V', A, 01066 $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', A, 01067 $ ' | / ( n ulp ) 4 = | I - VV', A, 01068 $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', A, 01069 $ ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A, 01070 $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A, 01071 $ ' | / ( n ulp ) 8 = | I - ZZ', A, 01072 $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A, 01073 $ ' l | / const. 10 = max | ( b H - a T )', A, 01074 $ ' l | / const.', / 01075 $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H', 01076 $ ' - a T ) r | / const.', / 1X ) 01077 * 01078 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 01079 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 01080 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 01081 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 ) 01082 * 01083 * End of CCHKGG 01084 * 01085 END
1.7.3 
