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