 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 00003 $ ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK, 00004 $ RESULT, INFO ) 00005 * 00006 * -- LAPACK test routine (version 3.1.1) -- 00007 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00008 * February 2007 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES, 00012 $ NTYPES 00013 REAL THRESH 00014 * .. 00015 * .. Array Arguments .. 00016 LOGICAL DOTYPE( * ) 00017 INTEGER ISEED( 4 ), NN( * ) 00018 REAL RESULT( * ), RWORK( * ) 00019 COMPLEX A( LDA, * ), ALPHA( * ), ALPHA1( * ), 00020 $ B( LDA, * ), BETA( * ), BETA1( * ), 00021 $ Q( LDQ, * ), QE( LDQE, * ), S( LDA, * ), 00022 $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * CDRGEV checks the nonsymmetric generalized eigenvalue problem driver 00029 * routine CGGEV. 00030 * 00031 * CGGEV 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 CDRGEV 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 CGGEV: 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 * CDRGES 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, CDRGEV 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 CDRGES to continue the same random number 00209 * sequence. 00210 * 00211 * THRESH (input) REAL 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) COMPLEX 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, and T. 00230 * It must be at least 1 and at least max( NN ). 00231 * 00232 * B (input/workspace) COMPLEX 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 array, dimension (LDA, max(NN)) 00238 * The Schur form matrix computed from A by CGGEV. On exit, S 00239 * contains the Schur form matrix corresponding to the matrix 00240 * in A. 00241 * 00242 * T (workspace) COMPLEX array, dimension (LDA, max(NN)) 00243 * The upper triangular matrix computed from B by CGGEV. 00244 * 00245 * Q (workspace) COMPLEX array, dimension (LDQ, max(NN)) 00246 * The (left) eigenvectors matrix computed by CGGEV. 00247 * 00248 * LDQ (input) INTEGER 00249 * The leading dimension of Q and Z. It must 00250 * be at least 1 and at least max( NN ). 00251 * 00252 * Z (workspace) COMPLEX array, dimension( LDQ, max(NN) ) 00253 * The (right) orthogonal matrix computed by CGGEV. 00254 * 00255 * QE (workspace) COMPLEX array, dimension( LDQ, max(NN) ) 00256 * QE holds the computed right or left eigenvectors. 00257 * 00258 * LDQE (input) INTEGER 00259 * The leading dimension of QE. LDQE >= max(1,max(NN)). 00260 * 00261 * ALPHA (workspace) COMPLEX array, dimension (max(NN)) 00262 * BETA (workspace) COMPLEX array, dimension (max(NN)) 00263 * The generalized eigenvalues of (A,B) computed by CGGEV. 00264 * ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th 00265 * generalized eigenvalue of A and B. 00266 * 00267 * ALPHA1 (workspace) COMPLEX array, dimension (max(NN)) 00268 * BETA1 (workspace) COMPLEX array, dimension (max(NN)) 00269 * Like ALPHAR, ALPHAI, BETA, these arrays contain the 00270 * eigenvalues of A and B, but those computed when CGGEV only 00271 * computes a partial eigendecomposition, i.e. not the 00272 * eigenvalues and left and right eigenvectors. 00273 * 00274 * WORK (workspace) COMPLEX array, dimension (LWORK) 00275 * 00276 * LWORK (input) INTEGER 00277 * The number of entries in WORK. LWORK >= N*(N+1) 00278 * 00279 * RWORK (workspace) REAL array, dimension (8*N) 00280 * Real workspace. 00281 * 00282 * RESULT (output) REAL array, dimension (2) 00283 * The values computed by the tests described above. 00284 * The values are currently limited to 1/ulp, to avoid overflow. 00285 * 00286 * INFO (output) INTEGER 00287 * = 0: successful exit 00288 * < 0: if INFO = -i, the i-th argument had an illegal value. 00289 * > 0: A routine returned an error code. INFO is the 00290 * absolute value of the INFO value returned. 00291 * 00292 * ===================================================================== 00293 * 00294 * .. Parameters .. 00295 REAL ZERO, ONE 00296 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00297 COMPLEX CZERO, CONE 00298 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00299 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00300 INTEGER MAXTYP 00301 PARAMETER ( MAXTYP = 26 ) 00302 * .. 00303 * .. Local Scalars .. 00304 LOGICAL BADNN 00305 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE, 00306 $ MAXWRK, MINWRK, MTYPES, N, N1, NB, NERRS, 00307 $ NMATS, NMAX, NTESTT 00308 REAL SAFMAX, SAFMIN, ULP, ULPINV 00309 COMPLEX CTEMP 00310 * .. 00311 * .. Local Arrays .. 00312 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP ) 00313 INTEGER 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 REAL RMAGN( 0: 3 ) 00319 * .. 00320 * .. External Functions .. 00321 INTEGER ILAENV 00322 REAL SLAMCH 00323 COMPLEX CLARND 00324 EXTERNAL ILAENV, SLAMCH, CLARND 00325 * .. 00326 * .. External Subroutines .. 00327 EXTERNAL ALASVM, CGET52, CGGEV, CLACPY, CLARFG, CLASET, 00328 $ CLATM4, CUNM2R, SLABAD, XERBLA 00329 * .. 00330 * .. Intrinsic Functions .. 00331 INTRINSIC ABS, CONJG, MAX, MIN, REAL, SIGN 00332 * .. 00333 * .. Data statements .. 00334 DATA KCLASS / 15*1, 10*2, 1*3 / 00335 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 00336 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 00337 DATA KADD / 0, 0, 0, 0, 3, 2 / 00338 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 00339 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 00340 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 00341 $ 1, 1, -4, 2, -4, 8*8, 0 / 00342 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 00343 $ 4*5, 4*3, 1 / 00344 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 00345 $ 4*6, 4*4, 1 / 00346 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 00347 $ 2, 1 / 00348 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 00349 $ 2, 1 / 00350 DATA KTRIAN / 16*0, 10*1 / 00351 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE., 00352 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE., 00353 $ 3*.FALSE., 5*.TRUE., .FALSE. / 00354 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE., 00355 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE., 00356 $ 9*.FALSE. / 00357 * .. 00358 * .. Executable Statements .. 00359 * 00360 * Check for errors 00361 * 00362 INFO = 0 00363 * 00364 BADNN = .FALSE. 00365 NMAX = 1 00366 DO 10 J = 1, NSIZES 00367 NMAX = MAX( NMAX, NN( J ) ) 00368 IF( NN( J ).LT.0 ) 00369 $ BADNN = .TRUE. 00370 10 CONTINUE 00371 * 00372 IF( NSIZES.LT.0 ) THEN 00373 INFO = -1 00374 ELSE IF( BADNN ) THEN 00375 INFO = -2 00376 ELSE IF( NTYPES.LT.0 ) THEN 00377 INFO = -3 00378 ELSE IF( THRESH.LT.ZERO ) THEN 00379 INFO = -6 00380 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 00381 INFO = -9 00382 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN 00383 INFO = -14 00384 ELSE IF( LDQE.LE.1 .OR. LDQE.LT.NMAX ) THEN 00385 INFO = -17 00386 END IF 00387 * 00388 * Compute workspace 00389 * (Note: Comments in the code beginning "Workspace:" describe the 00390 * minimal amount of workspace needed at that point in the code, 00391 * as well as the preferred amount for good performance. 00392 * NB refers to the optimal block size for the immediately 00393 * following subroutine, as returned by ILAENV. 00394 * 00395 MINWRK = 1 00396 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 00397 MINWRK = NMAX*( NMAX+1 ) 00398 NB = MAX( 1, ILAENV( 1, 'CGEQRF', ' ', NMAX, NMAX, -1, -1 ), 00399 $ ILAENV( 1, 'CUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ), 00400 $ ILAENV( 1, 'CUNGQR', ' ', NMAX, NMAX, NMAX, -1 ) ) 00401 MAXWRK = MAX( 2*NMAX, NMAX*( NB+1 ), NMAX*( NMAX+1 ) ) 00402 WORK( 1 ) = MAXWRK 00403 END IF 00404 * 00405 IF( LWORK.LT.MINWRK ) 00406 $ INFO = -23 00407 * 00408 IF( INFO.NE.0 ) THEN 00409 CALL XERBLA( 'CDRGEV', -INFO ) 00410 RETURN 00411 END IF 00412 * 00413 * Quick return if possible 00414 * 00415 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00416 $ RETURN 00417 * 00418 ULP = SLAMCH( 'Precision' ) 00419 SAFMIN = SLAMCH( 'Safe minimum' ) 00420 SAFMIN = SAFMIN / ULP 00421 SAFMAX = ONE / SAFMIN 00422 CALL SLABAD( SAFMIN, SAFMAX ) 00423 ULPINV = ONE / ULP 00424 * 00425 * The values RMAGN(2:3) depend on N, see below. 00426 * 00427 RMAGN( 0 ) = ZERO 00428 RMAGN( 1 ) = ONE 00429 * 00430 * Loop over sizes, types 00431 * 00432 NTESTT = 0 00433 NERRS = 0 00434 NMATS = 0 00435 * 00436 DO 220 JSIZE = 1, NSIZES 00437 N = NN( JSIZE ) 00438 N1 = MAX( 1, N ) 00439 RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 ) 00440 RMAGN( 3 ) = SAFMIN*ULPINV*N1 00441 * 00442 IF( NSIZES.NE.1 ) THEN 00443 MTYPES = MIN( MAXTYP, NTYPES ) 00444 ELSE 00445 MTYPES = MIN( MAXTYP+1, NTYPES ) 00446 END IF 00447 * 00448 DO 210 JTYPE = 1, MTYPES 00449 IF( .NOT.DOTYPE( JTYPE ) ) 00450 $ GO TO 210 00451 NMATS = NMATS + 1 00452 * 00453 * Save ISEED in case of an error. 00454 * 00455 DO 20 J = 1, 4 00456 IOLDSD( J ) = ISEED( J ) 00457 20 CONTINUE 00458 * 00459 * Generate test matrices A and B 00460 * 00461 * Description of control parameters: 00462 * 00463 * KCLASS: =1 means w/o rotation, =2 means w/ rotation, 00464 * =3 means random. 00465 * KATYPE: the "type" to be passed to CLATM4 for computing A. 00466 * KAZERO: the pattern of zeros on the diagonal for A: 00467 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 00468 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 00469 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 00470 * non-zero entries.) 00471 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 00472 * =2: large, =3: small. 00473 * LASIGN: .TRUE. if the diagonal elements of A are to be 00474 * multiplied by a random magnitude 1 number. 00475 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. 00476 * KTRIAN: =0: don't fill in the upper triangle, =1: do. 00477 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 00478 * RMAGN: used to implement KAMAGN and KBMAGN. 00479 * 00480 IF( MTYPES.GT.MAXTYP ) 00481 $ GO TO 100 00482 IERR = 0 00483 IF( KCLASS( JTYPE ).LT.3 ) THEN 00484 * 00485 * Generate A (w/o rotation) 00486 * 00487 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 00488 IN = 2*( ( N-1 ) / 2 ) + 1 00489 IF( IN.NE.N ) 00490 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 00491 ELSE 00492 IN = N 00493 END IF 00494 CALL CLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 00495 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 00496 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 00497 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 00498 $ ISEED, A, LDA ) 00499 IADD = KADD( KAZERO( JTYPE ) ) 00500 IF( IADD.GT.0 .AND. IADD.LE.N ) 00501 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 00502 * 00503 * Generate B (w/o rotation) 00504 * 00505 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 00506 IN = 2*( ( N-1 ) / 2 ) + 1 00507 IF( IN.NE.N ) 00508 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 00509 ELSE 00510 IN = N 00511 END IF 00512 CALL CLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 00513 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 00514 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 00515 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 00516 $ ISEED, B, LDA ) 00517 IADD = KADD( KBZERO( JTYPE ) ) 00518 IF( IADD.NE.0 .AND. IADD.LE.N ) 00519 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 00520 * 00521 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 00522 * 00523 * Include rotations 00524 * 00525 * Generate Q, Z as Householder transformations times 00526 * a diagonal matrix. 00527 * 00528 DO 40 JC = 1, N - 1 00529 DO 30 JR = JC, N 00530 Q( JR, JC ) = CLARND( 3, ISEED ) 00531 Z( JR, JC ) = CLARND( 3, ISEED ) 00532 30 CONTINUE 00533 CALL CLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 00534 $ WORK( JC ) ) 00535 WORK( 2*N+JC ) = SIGN( ONE, REAL( Q( JC, JC ) ) ) 00536 Q( JC, JC ) = CONE 00537 CALL CLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 00538 $ WORK( N+JC ) ) 00539 WORK( 3*N+JC ) = SIGN( ONE, REAL( Z( JC, JC ) ) ) 00540 Z( JC, JC ) = CONE 00541 40 CONTINUE 00542 CTEMP = CLARND( 3, ISEED ) 00543 Q( N, N ) = CONE 00544 WORK( N ) = CZERO 00545 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 00546 CTEMP = CLARND( 3, ISEED ) 00547 Z( N, N ) = CONE 00548 WORK( 2*N ) = CZERO 00549 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 00550 * 00551 * Apply the diagonal matrices 00552 * 00553 DO 60 JC = 1, N 00554 DO 50 JR = 1, N 00555 A( JR, JC ) = WORK( 2*N+JR )* 00556 $ CONJG( WORK( 3*N+JC ) )* 00557 $ A( JR, JC ) 00558 B( JR, JC ) = WORK( 2*N+JR )* 00559 $ CONJG( WORK( 3*N+JC ) )* 00560 $ B( JR, JC ) 00561 50 CONTINUE 00562 60 CONTINUE 00563 CALL CUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A, 00564 $ LDA, WORK( 2*N+1 ), IERR ) 00565 IF( IERR.NE.0 ) 00566 $ GO TO 90 00567 CALL CUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 00568 $ A, LDA, WORK( 2*N+1 ), IERR ) 00569 IF( IERR.NE.0 ) 00570 $ GO TO 90 00571 CALL CUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B, 00572 $ LDA, WORK( 2*N+1 ), IERR ) 00573 IF( IERR.NE.0 ) 00574 $ GO TO 90 00575 CALL CUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 00576 $ B, LDA, WORK( 2*N+1 ), IERR ) 00577 IF( IERR.NE.0 ) 00578 $ GO TO 90 00579 END IF 00580 ELSE 00581 * 00582 * Random matrices 00583 * 00584 DO 80 JC = 1, N 00585 DO 70 JR = 1, N 00586 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 00587 $ CLARND( 4, ISEED ) 00588 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 00589 $ CLARND( 4, ISEED ) 00590 70 CONTINUE 00591 80 CONTINUE 00592 END IF 00593 * 00594 90 CONTINUE 00595 * 00596 IF( IERR.NE.0 ) THEN 00597 WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE, 00598 $ IOLDSD 00599 INFO = ABS( IERR ) 00600 RETURN 00601 END IF 00602 * 00603 100 CONTINUE 00604 * 00605 DO 110 I = 1, 7 00606 RESULT( I ) = -ONE 00607 110 CONTINUE 00608 * 00609 * Call CGGEV to compute eigenvalues and eigenvectors. 00610 * 00611 CALL CLACPY( ' ', N, N, A, LDA, S, LDA ) 00612 CALL CLACPY( ' ', N, N, B, LDA, T, LDA ) 00613 CALL CGGEV( 'V', 'V', N, S, LDA, T, LDA, ALPHA, BETA, Q, 00614 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR ) 00615 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00616 RESULT( 1 ) = ULPINV 00617 WRITE( NOUNIT, FMT = 9999 )'CGGEV1', IERR, N, JTYPE, 00618 $ IOLDSD 00619 INFO = ABS( IERR ) 00620 GO TO 190 00621 END IF 00622 * 00623 * Do the tests (1) and (2) 00624 * 00625 CALL CGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHA, BETA, 00626 $ WORK, RWORK, RESULT( 1 ) ) 00627 IF( RESULT( 2 ).GT.THRESH ) THEN 00628 WRITE( NOUNIT, FMT = 9998 )'Left', 'CGGEV1', 00629 $ RESULT( 2 ), N, JTYPE, IOLDSD 00630 END IF 00631 * 00632 * Do the tests (3) and (4) 00633 * 00634 CALL CGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHA, 00635 $ BETA, WORK, RWORK, RESULT( 3 ) ) 00636 IF( RESULT( 4 ).GT.THRESH ) THEN 00637 WRITE( NOUNIT, FMT = 9998 )'Right', 'CGGEV1', 00638 $ RESULT( 4 ), N, JTYPE, IOLDSD 00639 END IF 00640 * 00641 * Do test (5) 00642 * 00643 CALL CLACPY( ' ', N, N, A, LDA, S, LDA ) 00644 CALL CLACPY( ' ', N, N, B, LDA, T, LDA ) 00645 CALL CGGEV( 'N', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, Q, 00646 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR ) 00647 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00648 RESULT( 1 ) = ULPINV 00649 WRITE( NOUNIT, FMT = 9999 )'CGGEV2', IERR, N, JTYPE, 00650 $ IOLDSD 00651 INFO = ABS( IERR ) 00652 GO TO 190 00653 END IF 00654 * 00655 DO 120 J = 1, N 00656 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 00657 $ BETA1( J ) )RESULT( 5 ) = ULPINV 00658 120 CONTINUE 00659 * 00660 * Do test (6): Compute eigenvalues and left eigenvectors, 00661 * and test them 00662 * 00663 CALL CLACPY( ' ', N, N, A, LDA, S, LDA ) 00664 CALL CLACPY( ' ', N, N, B, LDA, T, LDA ) 00665 CALL CGGEV( 'V', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, QE, 00666 $ LDQE, Z, LDQ, WORK, LWORK, RWORK, IERR ) 00667 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00668 RESULT( 1 ) = ULPINV 00669 WRITE( NOUNIT, FMT = 9999 )'CGGEV3', IERR, N, JTYPE, 00670 $ IOLDSD 00671 INFO = ABS( IERR ) 00672 GO TO 190 00673 END IF 00674 * 00675 DO 130 J = 1, N 00676 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 00677 $ BETA1( J ) )RESULT( 6 ) = ULPINV 00678 130 CONTINUE 00679 * 00680 DO 150 J = 1, N 00681 DO 140 JC = 1, N 00682 IF( Q( J, JC ).NE.QE( J, JC ) ) 00683 $ RESULT( 6 ) = ULPINV 00684 140 CONTINUE 00685 150 CONTINUE 00686 * 00687 * Do test (7): Compute eigenvalues and right eigenvectors, 00688 * and test them 00689 * 00690 CALL CLACPY( ' ', N, N, A, LDA, S, LDA ) 00691 CALL CLACPY( ' ', N, N, B, LDA, T, LDA ) 00692 CALL CGGEV( 'N', 'V', N, S, LDA, T, LDA, ALPHA1, BETA1, Q, 00693 $ LDQ, QE, LDQE, WORK, LWORK, RWORK, IERR ) 00694 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00695 RESULT( 1 ) = ULPINV 00696 WRITE( NOUNIT, FMT = 9999 )'CGGEV4', IERR, N, JTYPE, 00697 $ IOLDSD 00698 INFO = ABS( IERR ) 00699 GO TO 190 00700 END IF 00701 * 00702 DO 160 J = 1, N 00703 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 00704 $ BETA1( J ) )RESULT( 7 ) = ULPINV 00705 160 CONTINUE 00706 * 00707 DO 180 J = 1, N 00708 DO 170 JC = 1, N 00709 IF( Z( J, JC ).NE.QE( J, JC ) ) 00710 $ RESULT( 7 ) = ULPINV 00711 170 CONTINUE 00712 180 CONTINUE 00713 * 00714 * End of Loop -- Check for RESULT(j) > THRESH 00715 * 00716 190 CONTINUE 00717 * 00718 NTESTT = NTESTT + 7 00719 * 00720 * Print out tests which fail. 00721 * 00722 DO 200 JR = 1, 7 00723 IF( RESULT( JR ).GE.THRESH ) THEN 00724 * 00725 * If this is the first test to fail, 00726 * print a header to the data file. 00727 * 00728 IF( NERRS.EQ.0 ) THEN 00729 WRITE( NOUNIT, FMT = 9997 )'CGV' 00730 * 00731 * Matrix types 00732 * 00733 WRITE( NOUNIT, FMT = 9996 ) 00734 WRITE( NOUNIT, FMT = 9995 ) 00735 WRITE( NOUNIT, FMT = 9994 )'Orthogonal' 00736 * 00737 * Tests performed 00738 * 00739 WRITE( NOUNIT, FMT = 9993 ) 00740 * 00741 END IF 00742 NERRS = NERRS + 1 00743 IF( RESULT( JR ).LT.10000.0 ) THEN 00744 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 00745 $ RESULT( JR ) 00746 ELSE 00747 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 00748 $ RESULT( JR ) 00749 END IF 00750 END IF 00751 200 CONTINUE 00752 * 00753 210 CONTINUE 00754 220 CONTINUE 00755 * 00756 * Summary 00757 * 00758 CALL ALASVM( 'CGV', NOUNIT, NERRS, NTESTT, 0 ) 00759 * 00760 WORK( 1 ) = MAXWRK 00761 * 00762 RETURN 00763 * 00764 9999 FORMAT( ' CDRGEV: ', A, ' returned INFO=', I6, '.', / 3X, 'N=', 00765 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00766 * 00767 9998 FORMAT( ' CDRGEV: ', A, ' Eigenvectors from ', A, ' incorrectly ', 00768 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 3X, 00769 $ 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 3( I4, ',' ), I5, 00770 $ ')' ) 00771 * 00772 9997 FORMAT( / 1X, A3, ' -- Complex Generalized eigenvalue problem ', 00773 $ 'driver' ) 00774 * 00775 9996 FORMAT( ' Matrix types (see CDRGEV for details): ' ) 00776 * 00777 9995 FORMAT( ' Special Matrices:', 23X, 00778 $ '(J''=transposed Jordan block)', 00779 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 00780 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 00781 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 00782 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 00783 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 00784 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 00785 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 00786 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 00787 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 00788 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 00789 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 00790 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 00791 $ '23=(small,large) 24=(small,small) 25=(large,large)', 00792 $ / ' 26=random O(1) matrices.' ) 00793 * 00794 9993 FORMAT( / ' Tests performed: ', 00795 $ / ' 1 = max | ( b A - a B )''*l | / const.,', 00796 $ / ' 2 = | |VR(i)| - 1 | / ulp,', 00797 $ / ' 3 = max | ( b A - a B )*r | / const.', 00798 $ / ' 4 = | |VL(i)| - 1 | / ulp,', 00799 $ / ' 5 = 0 if W same no matter if r or l computed,', 00800 $ / ' 6 = 0 if l same no matter if l computed,', 00801 $ / ' 7 = 0 if r same no matter if r computed,', / 1X ) 00802 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00803 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 00804 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00805 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 ) 00806 * 00807 * End of CDRGEV 00808 * 00809 END
1.7.3 
