 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZDRVGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ THRSHN, NOUNIT, A, LDA, B, S, T, S2, T2, Q, 00003 $ LDQ, Z, ALPHA1, BETA1, ALPHA2, BETA2, VL, VR, 00004 $ WORK, LWORK, RWORK, RESULT, INFO ) 00005 * 00006 * -- LAPACK test routine (version 3.1) -- 00007 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES 00012 DOUBLE PRECISION THRESH, THRSHN 00013 * .. 00014 * .. Array Arguments .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * ZDRVGG checks the nonsymmetric generalized eigenvalue driver 00020 * routines. 00021 * T T T 00022 * ZGEGS factors A and B as Q S Z and Q T Z , where means 00023 * transpose, T is upper triangular, S is in generalized Schur form 00024 * (upper triangular), and Q and Z are unitary. It also 00025 * computes the generalized eigenvalues (alpha(1),beta(1)), ..., 00026 * (alpha(n),beta(n)), where alpha(j)=S(j,j) and beta(j)=T(j,j) -- 00027 * thus, w(j) = alpha(j)/beta(j) is a root of the generalized 00028 * eigenvalue problem 00029 * 00030 * det( A - w(j) B ) = 0 00031 * 00032 * and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent 00033 * problem 00034 * 00035 * det( m(j) A - B ) = 0 00036 * 00037 * ZGEGV computes the generalized eigenvalues (alpha(1),beta(1)), ..., 00038 * (alpha(n),beta(n)), the matrix L whose columns contain the 00039 * generalized left eigenvectors l, and the matrix R whose columns 00040 * contain the generalized right eigenvectors r for the pair (A,B). 00041 * 00042 * When ZDRVGG is called, a number of matrix "sizes" ("n's") and a 00043 * number of matrix "types" are specified. For each size ("n") 00044 * and each type of matrix, one matrix will be generated and used 00045 * to test the nonsymmetric eigenroutines. For each matrix, 7 00046 * tests will be performed and compared with the threshhold THRESH: 00047 * 00048 * Results from ZGEGS: 00049 * 00050 * H 00051 * (1) | A - Q S Z | / ( |A| n ulp ) 00052 * 00053 * H 00054 * (2) | B - Q T Z | / ( |B| n ulp ) 00055 * 00056 * H 00057 * (3) | I - QQ | / ( n ulp ) 00058 * 00059 * H 00060 * (4) | I - ZZ | / ( n ulp ) 00061 * 00062 * (5) maximum over j of D(j) where: 00063 * 00064 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)| 00065 * D(j) = ------------------------ + ----------------------- 00066 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|) 00067 * 00068 * Results from ZGEGV: 00069 * 00070 * (6) max over all left eigenvalue/-vector pairs (beta/alpha,l) of 00071 * 00072 * | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) ) 00073 * 00074 * where l**H is the conjugate tranpose of l. 00075 * 00076 * (7) max over all right eigenvalue/-vector pairs (beta/alpha,r) of 00077 * 00078 * | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) ) 00079 * 00080 * Test Matrices 00081 * ---- -------- 00082 * 00083 * The sizes of the test matrices are specified by an array 00084 * NN(1:NSIZES); the value of each element NN(j) specifies one size. 00085 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if 00086 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00087 * Currently, the list of possible types is: 00088 * 00089 * (1) ( 0, 0 ) (a pair of zero matrices) 00090 * 00091 * (2) ( I, 0 ) (an identity and a zero matrix) 00092 * 00093 * (3) ( 0, I ) (an identity and a zero matrix) 00094 * 00095 * (4) ( I, I ) (a pair of identity matrices) 00096 * 00097 * t t 00098 * (5) ( J , J ) (a pair of transposed Jordan blocks) 00099 * 00100 * t ( I 0 ) 00101 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t ) 00102 * ( 0 I ) ( 0 J ) 00103 * and I is a k x k identity and J a (k+1)x(k+1) 00104 * Jordan block; k=(N-1)/2 00105 * 00106 * (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal 00107 * matrix with those diagonal entries.) 00108 * (8) ( I, D ) 00109 * 00110 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big 00111 * 00112 * (10) ( small*D, big*I ) 00113 * 00114 * (11) ( big*I, small*D ) 00115 * 00116 * (12) ( small*I, big*D ) 00117 * 00118 * (13) ( big*D, big*I ) 00119 * 00120 * (14) ( small*D, small*I ) 00121 * 00122 * (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and 00123 * D2 is diag( 0, N-3, N-4,..., 1, 0, 0 ) 00124 * t t 00125 * (16) Q ( J , J ) Z where Q and Z are random unitary matrices. 00126 * 00127 * (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices 00128 * with random O(1) entries above the diagonal 00129 * and diagonal entries diag(T1) = 00130 * ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) = 00131 * ( 0, N-3, N-4,..., 1, 0, 0 ) 00132 * 00133 * (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 ) 00134 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 ) 00135 * s = machine precision. 00136 * 00137 * (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 ) 00138 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 ) 00139 * 00140 * N-5 00141 * (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 ) 00142 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00143 * 00144 * (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 ) 00145 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00146 * where r1,..., r(N-4) are random. 00147 * 00148 * (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00149 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00150 * 00151 * (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00152 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00153 * 00154 * (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00155 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00156 * 00157 * (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00158 * diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00159 * 00160 * (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular 00161 * matrices. 00162 * 00163 * Arguments 00164 * ========= 00165 * 00166 * NSIZES (input) INTEGER 00167 * The number of sizes of matrices to use. If it is zero, 00168 * ZDRVGG does nothing. It must be at least zero. 00169 * 00170 * NN (input) INTEGER array, dimension (NSIZES) 00171 * An array containing the sizes to be used for the matrices. 00172 * Zero values will be skipped. The values must be at least 00173 * zero. 00174 * 00175 * NTYPES (input) INTEGER 00176 * The number of elements in DOTYPE. If it is zero, ZDRVGG 00177 * does nothing. It must be at least zero. If it is MAXTYP+1 00178 * and NSIZES is 1, then an additional type, MAXTYP+1 is 00179 * defined, which is to use whatever matrix is in A. This 00180 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00181 * DOTYPE(MAXTYP+1) is .TRUE. . 00182 * 00183 * DOTYPE (input) 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 * 00192 * ISEED (input/output) INTEGER array, dimension (4) 00193 * On entry ISEED specifies the seed of the random number 00194 * generator. The array elements should be between 0 and 4095; 00195 * if not they will be reduced mod 4096. Also, ISEED(4) must 00196 * be odd. The random number generator uses a linear 00197 * congruential sequence limited to small integers, and so 00198 * should produce machine independent random numbers. The 00199 * values of ISEED are changed on exit, and can be used in the 00200 * next call to ZDRVGG to continue the same random number 00201 * sequence. 00202 * 00203 * THRESH (input) DOUBLE PRECISION 00204 * A test will count as "failed" if the "error", computed as 00205 * described above, exceeds THRESH. Note that the error is 00206 * scaled to be O(1), so THRESH should be a reasonably small 00207 * multiple of 1, e.g., 10 or 100. In particular, it should 00208 * not depend on the precision (single vs. double) or the size 00209 * of the matrix. It must be at least zero. 00210 * 00211 * THRSHN (input) DOUBLE PRECISION 00212 * Threshhold for reporting eigenvector normalization error. 00213 * If the normalization of any eigenvector differs from 1 by 00214 * more than THRSHN*ulp, then a special error message will be 00215 * printed. (This is handled separately from the other tests, 00216 * since only a compiler or programming error should cause an 00217 * error message, at least if THRSHN is at least 5--10.) 00218 * 00219 * NOUNIT (input) INTEGER 00220 * The FORTRAN unit number for printing out error messages 00221 * (e.g., if a routine returns IINFO not equal to 0.) 00222 * 00223 * A (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN)) 00224 * Used to hold the original A matrix. Used as input only 00225 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00226 * DOTYPE(MAXTYP+1)=.TRUE. 00227 * 00228 * LDA (input) INTEGER 00229 * The leading dimension of A, B, S, T, S2, and T2. 00230 * It must be at least 1 and at least max( NN ). 00231 * 00232 * B (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN)) 00233 * Used to hold the original B matrix. Used as input only 00234 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00235 * DOTYPE(MAXTYP+1)=.TRUE. 00236 * 00237 * S (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) 00238 * The upper triangular matrix computed from A by ZGEGS. 00239 * 00240 * T (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) 00241 * The upper triangular matrix computed from B by ZGEGS. 00242 * 00243 * S2 (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) 00244 * The matrix computed from A by ZGEGV. This will be the 00245 * Schur (upper triangular) form of some matrix related to A, 00246 * but will not, in general, be the same as S. 00247 * 00248 * T2 (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) 00249 * The matrix computed from B by ZGEGV. This will be the 00250 * Schur form of some matrix related to B, but will not, in 00251 * general, be the same as T. 00252 * 00253 * Q (workspace) COMPLEX*16 array, dimension (LDQ, max(NN)) 00254 * The (left) unitary matrix computed by ZGEGS. 00255 * 00256 * LDQ (input) INTEGER 00257 * The leading dimension of Q, Z, VL, and VR. It must 00258 * be at least 1 and at least max( NN ). 00259 * 00260 * Z (workspace) COMPLEX*16 array, dimension (LDQ, max(NN)) 00261 * The (right) unitary matrix computed by ZGEGS. 00262 * 00263 * ALPHA1 (workspace) COMPLEX*16 array, dimension (max(NN)) 00264 * BETA1 (workspace) COMPLEX*16 array, dimension (max(NN)) 00265 * The generalized eigenvalues of (A,B) computed by ZGEGS. 00266 * ALPHA1(k) / BETA1(k) is the k-th generalized eigenvalue of 00267 * the matrices in A and B. 00268 * 00269 * ALPHA2 (workspace) COMPLEX*16 array, dimension (max(NN)) 00270 * BETA2 (workspace) COMPLEX*16 array, dimension (max(NN)) 00271 * The generalized eigenvalues of (A,B) computed by ZGEGV. 00272 * ALPHA2(k) / BETA2(k) is the k-th generalized eigenvalue of 00273 * the matrices in A and B. 00274 * 00275 * VL (workspace) COMPLEX*16 array, dimension (LDQ, max(NN)) 00276 * The (lower triangular) left eigenvector matrix for the 00277 * matrices in A and B. 00278 * 00279 * VR (workspace) COMPLEX*16 array, dimension (LDQ, max(NN)) 00280 * The (upper triangular) right eigenvector matrix for the 00281 * matrices in A and B. 00282 * 00283 * WORK (workspace) COMPLEX*16 array, dimension (LWORK) 00284 * 00285 * LWORK (input) INTEGER 00286 * The number of entries in WORK. This must be at least 00287 * MAX( 2*N, N*(NB+1), (k+1)*(2*k+N+1) ), where "k" is the 00288 * sum of the blocksize and number-of-shifts for ZHGEQZ, and 00289 * NB is the greatest of the blocksizes for ZGEQRF, ZUNMQR, 00290 * and ZUNGQR. (The blocksizes and the number-of-shifts are 00291 * retrieved through calls to ILAENV.) 00292 * 00293 * RWORK (workspace) DOUBLE PRECISION array, dimension (8*N) 00294 * 00295 * RESULT (output) DOUBLE PRECISION array, dimension (7) 00296 * The values computed by the tests described above. 00297 * The values are currently limited to 1/ulp, to avoid 00298 * overflow. 00299 * 00300 * INFO (output) INTEGER 00301 * = 0: successful exit 00302 * < 0: if INFO = -i, the i-th argument had an illegal value. 00303 * > 0: A routine returned an error code. INFO is the 00304 * absolute value of the INFO value returned. 00305 * 00306 * ===================================================================== 00307 * 00308 LOGICAL DOTYPE( * ) 00309 INTEGER ISEED( 4 ), NN( * ) 00310 DOUBLE PRECISION RESULT( * ), RWORK( * ) 00311 COMPLEX*16 A( LDA, * ), ALPHA1( * ), ALPHA2( * ), 00312 $ B( LDA, * ), BETA1( * ), BETA2( * ), 00313 $ Q( LDQ, * ), S( LDA, * ), S2( LDA, * ), 00314 $ T( LDA, * ), T2( LDA, * ), VL( LDQ, * ), 00315 $ VR( LDQ, * ), WORK( * ), Z( LDQ, * ) 00316 * .. 00317 * .. Parameters .. 00318 DOUBLE PRECISION ZERO, ONE 00319 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00320 COMPLEX*16 CZERO, CONE 00321 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00322 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00323 INTEGER MAXTYP 00324 PARAMETER ( MAXTYP = 26 ) 00325 * .. 00326 * .. Local Scalars .. 00327 LOGICAL BADNN 00328 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE, 00329 $ LWKOPT, MTYPES, N, N1, NB, NBZ, NERRS, NMATS, 00330 $ NMAX, NS, NTEST, NTESTT 00331 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV 00332 COMPLEX*16 CTEMP, X 00333 * .. 00334 * .. Local Arrays .. 00335 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP ) 00336 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 00337 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 00338 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 00339 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 00340 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 00341 DOUBLE PRECISION DUMMA( 4 ), RMAGN( 0: 3 ) 00342 * .. 00343 * .. External Functions .. 00344 INTEGER ILAENV 00345 DOUBLE PRECISION DLAMCH 00346 COMPLEX*16 ZLARND 00347 EXTERNAL ILAENV, DLAMCH, ZLARND 00348 * .. 00349 * .. External Subroutines .. 00350 EXTERNAL ALASVM, DLABAD, XERBLA, ZGEGS, ZGEGV, ZGET51, 00351 $ ZGET52, ZLACPY, ZLARFG, ZLASET, ZLATM4, ZUNM2R 00352 * .. 00353 * .. Intrinsic Functions .. 00354 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SIGN 00355 * .. 00356 * .. Statement Functions .. 00357 DOUBLE PRECISION ABS1 00358 * .. 00359 * .. Statement Function definitions .. 00360 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) 00361 * .. 00362 * .. Data statements .. 00363 DATA KCLASS / 15*1, 10*2, 1*3 / 00364 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 00365 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 00366 DATA KADD / 0, 0, 0, 0, 3, 2 / 00367 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 00368 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 00369 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 00370 $ 1, 1, -4, 2, -4, 8*8, 0 / 00371 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 00372 $ 4*5, 4*3, 1 / 00373 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 00374 $ 4*6, 4*4, 1 / 00375 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 00376 $ 2, 1 / 00377 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 00378 $ 2, 1 / 00379 DATA KTRIAN / 16*0, 10*1 / 00380 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE., 00381 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE., 00382 $ 3*.FALSE., 5*.TRUE., .FALSE. / 00383 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE., 00384 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE., 00385 $ 9*.FALSE. / 00386 * .. 00387 * .. Executable Statements .. 00388 * 00389 * Check for errors 00390 * 00391 INFO = 0 00392 * 00393 BADNN = .FALSE. 00394 NMAX = 1 00395 DO 10 J = 1, NSIZES 00396 NMAX = MAX( NMAX, NN( J ) ) 00397 IF( NN( J ).LT.0 ) 00398 $ BADNN = .TRUE. 00399 10 CONTINUE 00400 * 00401 * Maximum blocksize and shift -- we assume that blocksize and number 00402 * of shifts are monotone increasing functions of N. 00403 * 00404 NB = MAX( 1, ILAENV( 1, 'ZGEQRF', ' ', NMAX, NMAX, -1, -1 ), 00405 $ ILAENV( 1, 'ZUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ), 00406 $ ILAENV( 1, 'ZUNGQR', ' ', NMAX, NMAX, NMAX, -1 ) ) 00407 NBZ = ILAENV( 1, 'ZHGEQZ', 'SII', NMAX, 1, NMAX, 0 ) 00408 NS = ILAENV( 4, 'ZHGEQZ', 'SII', NMAX, 1, NMAX, 0 ) 00409 I1 = NBZ + NS 00410 LWKOPT = MAX( 2*NMAX, NMAX*( NB+1 ), ( 2*I1+NMAX+1 )*( I1+1 ) ) 00411 * 00412 * Check for errors 00413 * 00414 IF( NSIZES.LT.0 ) THEN 00415 INFO = -1 00416 ELSE IF( BADNN ) THEN 00417 INFO = -2 00418 ELSE IF( NTYPES.LT.0 ) THEN 00419 INFO = -3 00420 ELSE IF( THRESH.LT.ZERO ) THEN 00421 INFO = -6 00422 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 00423 INFO = -10 00424 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN 00425 INFO = -19 00426 ELSE IF( LWKOPT.GT.LWORK ) THEN 00427 INFO = -30 00428 END IF 00429 * 00430 IF( INFO.NE.0 ) THEN 00431 CALL XERBLA( 'ZDRVGG', -INFO ) 00432 RETURN 00433 END IF 00434 * 00435 * Quick return if possible 00436 * 00437 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00438 $ RETURN 00439 * 00440 ULP = DLAMCH( 'Precision' ) 00441 SAFMIN = DLAMCH( 'Safe minimum' ) 00442 SAFMIN = SAFMIN / ULP 00443 SAFMAX = ONE / SAFMIN 00444 CALL DLABAD( SAFMIN, SAFMAX ) 00445 ULPINV = ONE / ULP 00446 * 00447 * The values RMAGN(2:3) depend on N, see below. 00448 * 00449 RMAGN( 0 ) = ZERO 00450 RMAGN( 1 ) = ONE 00451 * 00452 * Loop over sizes, types 00453 * 00454 NTESTT = 0 00455 NERRS = 0 00456 NMATS = 0 00457 * 00458 DO 160 JSIZE = 1, NSIZES 00459 N = NN( JSIZE ) 00460 N1 = MAX( 1, N ) 00461 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 ) 00462 RMAGN( 3 ) = SAFMIN*ULPINV*N1 00463 * 00464 IF( NSIZES.NE.1 ) THEN 00465 MTYPES = MIN( MAXTYP, NTYPES ) 00466 ELSE 00467 MTYPES = MIN( MAXTYP+1, NTYPES ) 00468 END IF 00469 * 00470 DO 150 JTYPE = 1, MTYPES 00471 IF( .NOT.DOTYPE( JTYPE ) ) 00472 $ GO TO 150 00473 NMATS = NMATS + 1 00474 NTEST = 0 00475 * 00476 * Save ISEED in case of an error. 00477 * 00478 DO 20 J = 1, 4 00479 IOLDSD( J ) = ISEED( J ) 00480 20 CONTINUE 00481 * 00482 * Initialize RESULT 00483 * 00484 DO 30 J = 1, 7 00485 RESULT( J ) = ZERO 00486 30 CONTINUE 00487 * 00488 * Compute A and B 00489 * 00490 * Description of control parameters: 00491 * 00492 * KZLASS: =1 means w/o rotation, =2 means w/ rotation, 00493 * =3 means random. 00494 * KATYPE: the "type" to be passed to ZLATM4 for computing A. 00495 * KAZERO: the pattern of zeros on the diagonal for A: 00496 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 00497 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 00498 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 00499 * non-zero entries.) 00500 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 00501 * =2: large, =3: small. 00502 * LASIGN: .TRUE. if the diagonal elements of A are to be 00503 * multiplied by a random magnitude 1 number. 00504 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B. 00505 * KTRIAN: =0: don't fill in the upper triangle, =1: do. 00506 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 00507 * RMAGN: used to implement KAMAGN and KBMAGN. 00508 * 00509 IF( MTYPES.GT.MAXTYP ) 00510 $ GO TO 110 00511 IINFO = 0 00512 IF( KCLASS( JTYPE ).LT.3 ) THEN 00513 * 00514 * Generate A (w/o rotation) 00515 * 00516 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 00517 IN = 2*( ( N-1 ) / 2 ) + 1 00518 IF( IN.NE.N ) 00519 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 00520 ELSE 00521 IN = N 00522 END IF 00523 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 00524 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 00525 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 00526 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 00527 $ ISEED, A, LDA ) 00528 IADD = KADD( KAZERO( JTYPE ) ) 00529 IF( IADD.GT.0 .AND. IADD.LE.N ) 00530 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 00531 * 00532 * Generate B (w/o rotation) 00533 * 00534 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 00535 IN = 2*( ( N-1 ) / 2 ) + 1 00536 IF( IN.NE.N ) 00537 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 00538 ELSE 00539 IN = N 00540 END IF 00541 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 00542 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 00543 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 00544 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 00545 $ ISEED, B, LDA ) 00546 IADD = KADD( KBZERO( JTYPE ) ) 00547 IF( IADD.NE.0 .AND. IADD.LE.N ) 00548 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 00549 * 00550 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 00551 * 00552 * Include rotations 00553 * 00554 * Generate Q, Z as Householder transformations times 00555 * a diagonal matrix. 00556 * 00557 DO 50 JC = 1, N - 1 00558 DO 40 JR = JC, N 00559 Q( JR, JC ) = ZLARND( 3, ISEED ) 00560 Z( JR, JC ) = ZLARND( 3, ISEED ) 00561 40 CONTINUE 00562 CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 00563 $ WORK( JC ) ) 00564 WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) ) 00565 Q( JC, JC ) = CONE 00566 CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 00567 $ WORK( N+JC ) ) 00568 WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) ) 00569 Z( JC, JC ) = CONE 00570 50 CONTINUE 00571 CTEMP = ZLARND( 3, ISEED ) 00572 Q( N, N ) = CONE 00573 WORK( N ) = CZERO 00574 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 00575 CTEMP = ZLARND( 3, ISEED ) 00576 Z( N, N ) = CONE 00577 WORK( 2*N ) = CZERO 00578 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 00579 * 00580 * Apply the diagonal matrices 00581 * 00582 DO 70 JC = 1, N 00583 DO 60 JR = 1, N 00584 A( JR, JC ) = WORK( 2*N+JR )* 00585 $ DCONJG( WORK( 3*N+JC ) )* 00586 $ A( JR, JC ) 00587 B( JR, JC ) = WORK( 2*N+JR )* 00588 $ DCONJG( WORK( 3*N+JC ) )* 00589 $ B( JR, JC ) 00590 60 CONTINUE 00591 70 CONTINUE 00592 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A, 00593 $ LDA, WORK( 2*N+1 ), IINFO ) 00594 IF( IINFO.NE.0 ) 00595 $ GO TO 100 00596 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 00597 $ A, LDA, WORK( 2*N+1 ), IINFO ) 00598 IF( IINFO.NE.0 ) 00599 $ GO TO 100 00600 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B, 00601 $ LDA, WORK( 2*N+1 ), IINFO ) 00602 IF( IINFO.NE.0 ) 00603 $ GO TO 100 00604 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 00605 $ B, LDA, WORK( 2*N+1 ), IINFO ) 00606 IF( IINFO.NE.0 ) 00607 $ GO TO 100 00608 END IF 00609 ELSE 00610 * 00611 * Random matrices 00612 * 00613 DO 90 JC = 1, N 00614 DO 80 JR = 1, N 00615 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 00616 $ ZLARND( 4, ISEED ) 00617 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 00618 $ ZLARND( 4, ISEED ) 00619 80 CONTINUE 00620 90 CONTINUE 00621 END IF 00622 * 00623 100 CONTINUE 00624 * 00625 IF( IINFO.NE.0 ) THEN 00626 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 00627 $ IOLDSD 00628 INFO = ABS( IINFO ) 00629 RETURN 00630 END IF 00631 * 00632 110 CONTINUE 00633 * 00634 * Call ZGEGS to compute H, T, Q, Z, alpha, and beta. 00635 * 00636 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 00637 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 00638 NTEST = 1 00639 RESULT( 1 ) = ULPINV 00640 * 00641 CALL ZGEGS( 'V', 'V', N, S, LDA, T, LDA, ALPHA1, BETA1, Q, 00642 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IINFO ) 00643 IF( IINFO.NE.0 ) THEN 00644 WRITE( NOUNIT, FMT = 9999 )'ZGEGS', IINFO, N, JTYPE, 00645 $ IOLDSD 00646 INFO = ABS( IINFO ) 00647 GO TO 130 00648 END IF 00649 * 00650 NTEST = 4 00651 * 00652 * Do tests 1--4 00653 * 00654 CALL ZGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ, WORK, 00655 $ RWORK, RESULT( 1 ) ) 00656 CALL ZGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ, WORK, 00657 $ RWORK, RESULT( 2 ) ) 00658 CALL ZGET51( 3, N, B, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK, 00659 $ RWORK, RESULT( 3 ) ) 00660 CALL ZGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK, 00661 $ RWORK, RESULT( 4 ) ) 00662 * 00663 * Do test 5: compare eigenvalues with diagonals. 00664 * 00665 TEMP1 = ZERO 00666 * 00667 DO 120 J = 1, N 00668 TEMP2 = ( ABS1( ALPHA1( J )-S( J, J ) ) / 00669 $ MAX( SAFMIN, ABS1( ALPHA1( J ) ), ABS1( S( J, 00670 $ J ) ) )+ABS1( BETA1( J )-T( J, J ) ) / 00671 $ MAX( SAFMIN, ABS1( BETA1( J ) ), ABS1( T( J, 00672 $ J ) ) ) ) / ULP 00673 TEMP1 = MAX( TEMP1, TEMP2 ) 00674 120 CONTINUE 00675 RESULT( 5 ) = TEMP1 00676 * 00677 * Call ZGEGV to compute S2, T2, VL, and VR, do tests. 00678 * 00679 * Eigenvalues and Eigenvectors 00680 * 00681 CALL ZLACPY( ' ', N, N, A, LDA, S2, LDA ) 00682 CALL ZLACPY( ' ', N, N, B, LDA, T2, LDA ) 00683 NTEST = 6 00684 RESULT( 6 ) = ULPINV 00685 * 00686 CALL ZGEGV( 'V', 'V', N, S2, LDA, T2, LDA, ALPHA2, BETA2, 00687 $ VL, LDQ, VR, LDQ, WORK, LWORK, RWORK, IINFO ) 00688 IF( IINFO.NE.0 ) THEN 00689 WRITE( NOUNIT, FMT = 9999 )'ZGEGV', IINFO, N, JTYPE, 00690 $ IOLDSD 00691 INFO = ABS( IINFO ) 00692 GO TO 130 00693 END IF 00694 * 00695 NTEST = 7 00696 * 00697 * Do Tests 6 and 7 00698 * 00699 CALL ZGET52( .TRUE., N, A, LDA, B, LDA, VL, LDQ, ALPHA2, 00700 $ BETA2, WORK, RWORK, DUMMA( 1 ) ) 00701 RESULT( 6 ) = DUMMA( 1 ) 00702 IF( DUMMA( 2 ).GT.THRSHN ) THEN 00703 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZGEGV', DUMMA( 2 ), 00704 $ N, JTYPE, IOLDSD 00705 END IF 00706 * 00707 CALL ZGET52( .FALSE., N, A, LDA, B, LDA, VR, LDQ, ALPHA2, 00708 $ BETA2, WORK, RWORK, DUMMA( 1 ) ) 00709 RESULT( 7 ) = DUMMA( 1 ) 00710 IF( DUMMA( 2 ).GT.THRESH ) THEN 00711 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZGEGV', DUMMA( 2 ), 00712 $ N, JTYPE, IOLDSD 00713 END IF 00714 * 00715 * End of Loop -- Check for RESULT(j) > THRESH 00716 * 00717 130 CONTINUE 00718 * 00719 NTESTT = NTESTT + NTEST 00720 * 00721 * Print out tests which fail. 00722 * 00723 DO 140 JR = 1, NTEST 00724 IF( RESULT( JR ).GE.THRESH ) THEN 00725 * 00726 * If this is the first test to fail, 00727 * print a header to the data file. 00728 * 00729 IF( NERRS.EQ.0 ) THEN 00730 WRITE( NOUNIT, FMT = 9997 )'ZGG' 00731 * 00732 * Matrix types 00733 * 00734 WRITE( NOUNIT, FMT = 9996 ) 00735 WRITE( NOUNIT, FMT = 9995 ) 00736 WRITE( NOUNIT, FMT = 9994 )'Unitary' 00737 * 00738 * Tests performed 00739 * 00740 WRITE( NOUNIT, FMT = 9993 )'unitary', '*', 00741 $ 'conjugate transpose', ( '*', J = 1, 5 ) 00742 * 00743 END IF 00744 NERRS = NERRS + 1 00745 IF( RESULT( JR ).LT.10000.0D0 ) THEN 00746 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 00747 $ RESULT( JR ) 00748 ELSE 00749 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 00750 $ RESULT( JR ) 00751 END IF 00752 END IF 00753 140 CONTINUE 00754 * 00755 150 CONTINUE 00756 160 CONTINUE 00757 * 00758 * Summary 00759 * 00760 CALL ALASVM( 'ZGG', NOUNIT, NERRS, NTESTT, 0 ) 00761 RETURN 00762 * 00763 9999 FORMAT( ' ZDRVGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00764 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00765 * 00766 9998 FORMAT( ' ZDRVGG: ', A, ' Eigenvectors from ', A, ' incorrectly ', 00767 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 00768 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 00769 $ ')' ) 00770 * 00771 9997 FORMAT( / 1X, A3, 00772 $ ' -- Complex Generalized eigenvalue problem driver' ) 00773 * 00774 9996 FORMAT( ' Matrix types (see ZDRVGG for details): ' ) 00775 * 00776 9995 FORMAT( ' Special Matrices:', 23X, 00777 $ '(J''=transposed Jordan block)', 00778 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 00779 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 00780 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 00781 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 00782 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 00783 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 00784 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 00785 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 00786 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 00787 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 00788 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 00789 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 00790 $ '23=(small,large) 24=(small,small) 25=(large,large)', 00791 $ / ' 26=random O(1) matrices.' ) 00792 * 00793 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ', 00794 $ 'Q and Z are ', A, ',', / 20X, 00795 $ 'l and r are the appropriate left and right', / 19X, 00796 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A, 00797 $ ' means ', A, '.)', / ' 1 = | A - Q S Z', A, 00798 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A, 00799 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A, 00800 $ ' | / ( n ulp ) 4 = | I - ZZ', A, 00801 $ ' | / ( n ulp )', / 00802 $ ' 5 = difference between (alpha,beta) and diagonals of', 00803 $ ' (S,T)', / ' 6 = max | ( b A - a B )', A, 00804 $ ' l | / const. 7 = max | ( b A - a B ) r | / const.', 00805 $ / 1X ) 00806 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00807 $ 4( I4, ',' ), ' result ', I3, ' is', 0P, F8.2 ) 00808 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00809 $ 4( I4, ',' ), ' result ', I3, ' is', 1P, D10.3 ) 00810 * 00811 * End of ZDRVGG 00812 * 00813 END
1.7.3 
