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