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