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