*DECK DDZRO SUBROUTINE DDZRO (AE, F, H, N, NQ, IROOT, RE, T, YH, UROUND, B, C, 8 FB, FC, Y) C***BEGIN PROLOGUE DDZRO C***SUBSIDIARY C***PURPOSE DDZRO SEARCHES FOR A ZERO OF A FUNCTION F(N, T, Y, IROOT) C BETWEEN THE GIVEN VALUES B AND C UNTIL THE WIDTH OF THE C INTERVAL (B, C) HAS COLLAPSED TO WITHIN A TOLERANCE C SPECIFIED BY THE STOPPING CRITERION, C ABS(B - C) .LE. 2.*(RW*ABS(B) + AE). C***LIBRARY SLATEC (SDRIVE) C***TYPE DOUBLE PRECISION (SDZRO-S, DDZRO-D, CDZRO-C) C***AUTHOR KAHANER, D. K., (NIST) C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C SUTHERLAND, C. D., (LANL) C MAIL STOP D466 C LOS ALAMOS NATIONAL LABORATORY C LOS ALAMOS, NM 87545 C***DESCRIPTION C C THIS IS A SPECIAL PURPOSE VERSION OF ZEROIN, MODIFIED FOR USE WITH C THE DDRIV PACKAGE. C C SANDIA MATHEMATICAL PROGRAM LIBRARY C MATHEMATICAL COMPUTING SERVICES DIVISION 5422 C SANDIA LABORATORIES C P. O. BOX 5800 C ALBUQUERQUE, NEW MEXICO 87115 C CONTROL DATA 6600 VERSION 4.5, 1 NOVEMBER 1971 C C PARAMETERS C F - NAME OF THE EXTERNAL FUNCTION, WHICH RETURNS A C DOUBLE PRECISION RESULT. THIS NAME MUST BE IN AN C EXTERNAL STATEMENT IN THE CALLING PROGRAM. C B - ONE END OF THE INTERVAL (B, C). THE VALUE RETURNED FOR C B USUALLY IS THE BETTER APPROXIMATION TO A ZERO OF F. C C - THE OTHER END OF THE INTERVAL (B, C). C RE - RELATIVE ERROR USED FOR RW IN THE STOPPING CRITERION. C IF THE REQUESTED RE IS LESS THAN MACHINE PRECISION, C THEN RW IS SET TO APPROXIMATELY MACHINE PRECISION. C AE - ABSOLUTE ERROR USED IN THE STOPPING CRITERION. IF THE C GIVEN INTERVAL (B, C) CONTAINS THE ORIGIN, THEN A C NONZERO VALUE SHOULD BE CHOSEN FOR AE. C C***REFERENCES L. F. SHAMPINE AND H. A. WATTS, ZEROIN, A ROOT-SOLVING C ROUTINE, SC-TM-70-631, SEPT 1970. C T. J. DEKKER, FINDING A ZERO BY MEANS OF SUCCESSIVE C LINEAR INTERPOLATION, CONSTRUCTIVE ASPECTS OF THE C FUNDAMENTAL THEOREM OF ALGEBRA, EDITED BY B. DEJON C AND P. HENRICI, 1969. C***ROUTINES CALLED DDNTP C***REVISION HISTORY (YYMMDD) C 790601 DATE WRITTEN C 900329 INITIAL SUBMISSION TO SLATEC. C***END PROLOGUE DDZRO INTEGER IC, IROOT, KOUNT, N, NQ DOUBLE PRECISION A, ACBS, ACMB, AE, B, C, CMB, ER, F, FA, FB, FC, 8 H, P, Q, RE, RW, T, TOL, UROUND, Y(*), YH(N,*) C***FIRST EXECUTABLE STATEMENT DDZRO ER = 4.D0*UROUND RW = MAX(RE, ER) IC = 0 ACBS = ABS(B - C) A = C FA = FC KOUNT = 0 C PERFORM INTERCHANGE 10 IF (ABS(FC) .LT. ABS(FB)) THEN A = B FA = FB B = C FB = FC C = A FC = FA END IF CMB = 0.5D0*(C - B) ACMB = ABS(CMB) TOL = RW*ABS(B) + AE C TEST STOPPING CRITERION IF (ACMB .LE. TOL) RETURN IF (KOUNT .GT. 50) RETURN C CALCULATE NEW ITERATE IMPLICITLY AS C B + P/Q, WHERE WE ARRANGE P .GE. 0. C THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW. P = (B - A)*FB Q = FA - FB IF (P .LT. 0.D0) THEN P = -P Q = -Q END IF C UPDATE A AND CHECK FOR SATISFACTORY REDUCTION C IN THE SIZE OF OUR BOUNDING INTERVAL. A = B FA = FB IC = IC + 1 IF (IC .GE. 4) THEN IF (8.D0*ACMB .GE. ACBS) THEN C BISECT B = 0.5D0*(C + B) GO TO 20 END IF IC = 0 END IF ACBS = ACMB C TEST FOR TOO SMALL A CHANGE IF (P .LE. ABS(Q)*TOL) THEN C INCREMENT BY TOLERANCE B = B + SIGN(TOL, CMB) C ROOT OUGHT TO BE BETWEEN C B AND (C + B)/2. ELSE IF (P .LT. CMB*Q) THEN C INTERPOLATE B = B + P/Q ELSE C BISECT B = 0.5D0*(C + B) END IF C HAVE COMPLETED COMPUTATION C FOR NEW ITERATE B. 20 CALL DDNTP (H, 0, N, NQ, T, B, YH, Y) FB = F(N, B, Y, IROOT) IF (N .EQ. 0) RETURN IF (FB .EQ. 0.D0) RETURN KOUNT = KOUNT + 1 C C DECIDE WHETHER NEXT STEP IS INTERPOLATION OR EXTRAPOLATION C IF (SIGN(1.0D0, FB) .EQ. SIGN(1.0D0, FC)) THEN C = A FC = FA END IF GO TO 10 END *DECK DDNTP SUBROUTINE DDNTP (H, K, N, NQ, T, TOUT, YH, Y) C***BEGIN PROLOGUE DDNTP C***SUBSIDIARY C***PURPOSE SUBROUTINE DDNTP INTERPOLATES THE K-TH DERIVATIVE OF Y AT C TOUT, USING THE DATA IN THE YH ARRAY. IF K HAS A VALUE C GREATER THAN NQ, THE NQ-TH DERIVATIVE IS CALCULATED. C***LIBRARY SLATEC (SDRIVE) C***TYPE DOUBLE PRECISION (SDNTP-S, DDNTP-D, CDNTP-C) C***AUTHOR KAHANER, D. K., (NIST) C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C SUTHERLAND, C. D., (LANL) C MAIL STOP D466 C LOS ALAMOS NATIONAL LABORATORY C LOS ALAMOS, NM 87545 C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 790601 DATE WRITTEN C 900329 INITIAL SUBMISSION TO SLATEC. C***END PROLOGUE DDNTP INTEGER I, J, JJ, K, KK, KUSED, N, NQ DOUBLE PRECISION FACTOR, H, R, T, TOUT, Y(*), YH(N,*) C***FIRST EXECUTABLE STATEMENT DDNTP IF (K .EQ. 0) THEN DO 10 I = 1,N 10 Y(I) = YH(I,NQ+1) R = ((TOUT - T)/H) DO 20 JJ = 1,NQ J = NQ + 1 - JJ DO 20 I = 1,N 20 Y(I) = YH(I,J) + R*Y(I) ELSE KUSED = MIN(K, NQ) FACTOR = 1.D0 DO 40 KK = 1,KUSED 40 FACTOR = FACTOR*(NQ+1-KK) DO 50 I = 1,N 50 Y(I) = FACTOR*YH(I,NQ+1) R = ((TOUT - T)/H) DO 80 JJ = KUSED+1,NQ J = KUSED + 1 + NQ - JJ FACTOR = 1.D0 DO 60 KK = 1,KUSED 60 FACTOR = FACTOR*(J-KK) DO 70 I = 1,N 70 Y(I) = FACTOR*YH(I,J) + R*Y(I) 80 CONTINUE DO 100 I = 1,N 100 Y(I) = Y(I)*H**(-KUSED) END IF RETURN END C------------------------------------------------------------------------ *DECK DLNGAM DOUBLE PRECISION FUNCTION DLNGAM (X) C***BEGIN PROLOGUE DLNGAM C***PURPOSE Compute the logarithm of the absolute value of the Gamma C function. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C7A C***TYPE DOUBLE PRECISION (ALNGAM-S, DLNGAM-D, CLNGAM-C) C***KEYWORDS ABSOLUTE VALUE, COMPLETE GAMMA FUNCTION, FNLIB, LOGARITHM, C SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C DLNGAM(X) calculates the double precision logarithm of the C absolute value of the Gamma function for double precision C argument X. C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, D9LGMC, DGAMMA, XERMSG C***REVISION HISTORY (YYMMDD) C 770601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900727 Added EXTERNAL statement. (WRB) C***END PROLOGUE DLNGAM DOUBLE PRECISION X, DXREL, PI, SINPIY, SQPI2L, SQ2PIL, XMAX, 1 Y, DGAMMA, D9LGMC, D1MACH, TEMP LOGICAL FIRST EXTERNAL DGAMMA SAVE SQ2PIL, SQPI2L, PI, XMAX, DXREL, FIRST DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 / DATA SQPI2L / +.2257913526 4472743236 3097614947 441 D+0 / DATA PI / 3.1415926535 8979323846 2643383279 50 D0 / DATA FIRST /.TRUE./ C***FIRST EXECUTABLE STATEMENT DLNGAM IF (FIRST) THEN TEMP = 1.D0/LOG(D1MACH(2)) XMAX = TEMP*D1MACH(2) DXREL = SQRT(D1MACH(4)) ENDIF FIRST = .FALSE. C Y = ABS (X) IF (Y.GT.10.D0) GO TO 20 C C LOG (ABS (DGAMMA(X)) ) FOR ABS(X) .LE. 10.0 C DLNGAM = LOG (ABS (DGAMMA(X)) ) RETURN C C LOG ( ABS (DGAMMA(X)) ) FOR ABS(X) .GT. 10.0 C 20 IF (Y .GT. XMAX) CALL XERMSG ('SLATEC', 'DLNGAM', + 'ABS(X) SO BIG DLNGAM OVERFLOWS', 2, 2) C IF (X.GT.0.D0) DLNGAM = SQ2PIL + (X-0.5D0)*LOG(X) - X + D9LGMC(Y) IF (X.GT.0.D0) RETURN C SINPIY = ABS (SIN(PI*Y)) IF (SINPIY .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DLNGAM', + 'X IS A NEGATIVE INTEGER', 3, 2) C IF (ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) CALL XERMSG ('SLATEC', + 'DLNGAM', + 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER', + 1, 1) C DLNGAM = SQPI2L + (X-0.5D0)*LOG(Y) - X - LOG(SINPIY) - D9LGMC(Y) RETURN C END *DECK D9LGMC DOUBLE PRECISION FUNCTION D9LGMC (X) C***BEGIN PROLOGUE D9LGMC C***SUBSIDIARY C***PURPOSE Compute the log Gamma correction factor so that C LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X C + D9LGMC(X). C***LIBRARY SLATEC (FNLIB) C***CATEGORY C7E C***TYPE DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C) C***KEYWORDS COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB, C LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C Compute the log gamma correction factor for X .GE. 10. so that C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X) C C Series for ALGM on the interval 0. to 1.00000E-02 C with weighted error 1.28E-31 C log weighted error 30.89 C significant figures required 29.81 C decimal places required 31.48 C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, DCSEVL, INITDS, XERMSG C***REVISION HISTORY (YYMMDD) C 770601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900720 Routine changed from user-callable to subsidiary. (WRB) C***END PROLOGUE D9LGMC DOUBLE PRECISION X, ALGMCS(15), XBIG, XMAX, DCSEVL, D1MACH LOGICAL FIRST SAVE ALGMCS, NALGM, XBIG, XMAX, FIRST DATA ALGMCS( 1) / +.1666389480 4518632472 0572965082 2 D+0 / DATA ALGMCS( 2) / -.1384948176 0675638407 3298605913 5 D-4 / DATA ALGMCS( 3) / +.9810825646 9247294261 5717154748 7 D-8 / DATA ALGMCS( 4) / -.1809129475 5724941942 6330626671 9 D-10 / DATA ALGMCS( 5) / +.6221098041 8926052271 2601554341 6 D-13 / DATA ALGMCS( 6) / -.3399615005 4177219443 0333059966 6 D-15 / DATA ALGMCS( 7) / +.2683181998 4826987489 5753884666 6 D-17 / DATA ALGMCS( 8) / -.2868042435 3346432841 4462239999 9 D-19 / DATA ALGMCS( 9) / +.3962837061 0464348036 7930666666 6 D-21 / DATA ALGMCS( 10) / -.6831888753 9857668701 1199999999 9 D-23 / DATA ALGMCS( 11) / +.1429227355 9424981475 7333333333 3 D-24 / DATA ALGMCS( 12) / -.3547598158 1010705471 9999999999 9 D-26 / DATA ALGMCS( 13) / +.1025680058 0104709120 0000000000 0 D-27 / DATA ALGMCS( 14) / -.3401102254 3167487999 9999999999 9 D-29 / DATA ALGMCS( 15) / +.1276642195 6300629333 3333333333 3 D-30 / DATA FIRST /.TRUE./ C***FIRST EXECUTABLE STATEMENT D9LGMC IF (FIRST) THEN NALGM = INITDS (ALGMCS, 15, REAL(D1MACH(3)) ) XBIG = 1.0D0/SQRT(D1MACH(3)) XMAX = EXP (MIN(LOG(D1MACH(2)/12.D0), -LOG(12.D0*D1MACH(1)))) ENDIF FIRST = .FALSE. C IF (X .LT. 10.D0) CALL XERMSG ('SLATEC', 'D9LGMC', + 'X MUST BE GE 10', 1, 2) IF (X.GE.XMAX) GO TO 20 C D9LGMC = 1.D0/(12.D0*X) IF (X.LT.XBIG) D9LGMC = DCSEVL (2.0D0*(10.D0/X)**2-1.D0, ALGMCS, 1 NALGM) / X RETURN C 20 D9LGMC = 0.D0 CALL XERMSG ('SLATEC', 'D9LGMC', 'X SO BIG D9LGMC UNDERFLOWS', 2, + 1) RETURN C END *DECK DCSEVL DOUBLE PRECISION FUNCTION DCSEVL (X, CS, N) C***BEGIN PROLOGUE DCSEVL C***PURPOSE Evaluate a Chebyshev series. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C3A2 C***TYPE DOUBLE PRECISION (CSEVL-S, DCSEVL-D) C***KEYWORDS CHEBYSHEV SERIES, FNLIB, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C Evaluate the N-term Chebyshev series CS at X. Adapted from C a method presented in the paper by Broucke referenced below. C C Input Arguments -- C X value at which the series is to be evaluated. C CS array of N terms of a Chebyshev series. In evaluating C CS, only half the first coefficient is summed. C N number of terms in array CS. C C***REFERENCES R. Broucke, Ten subroutines for the manipulation of C Chebyshev series, Algorithm 446, Communications of C the A.C.M. 16, (1973) pp. 254-256. C L. Fox and I. B. Parker, Chebyshev Polynomials in C Numerical Analysis, Oxford University Press, 1968, C page 56. C***ROUTINES CALLED D1MACH, XERMSG C***REVISION HISTORY (YYMMDD) C 770401 DATE WRITTEN C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900329 Prologued revised extensively and code rewritten to allow C X to be slightly outside interval (-1,+1). (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DCSEVL DOUBLE PRECISION B0, B1, B2, CS(*), ONEPL, TWOX, X, D1MACH LOGICAL FIRST SAVE FIRST, ONEPL DATA FIRST /.TRUE./ C***FIRST EXECUTABLE STATEMENT DCSEVL IF (FIRST) ONEPL = 1.0D0 + D1MACH(4) FIRST = .FALSE. IF (N .LT. 1) CALL XERMSG ('SLATEC', 'DCSEVL', + 'NUMBER OF TERMS .LE. 0', 2, 2) IF (N .GT. 1000) CALL XERMSG ('SLATEC', 'DCSEVL', + 'NUMBER OF TERMS .GT. 1000', 3, 2) IF (ABS(X) .GT. ONEPL) CALL XERMSG ('SLATEC', 'DCSEVL', + 'X OUTSIDE THE INTERVAL (-1,+1)', 1, 1) C B1 = 0.0D0 B0 = 0.0D0 TWOX = 2.0D0*X DO 10 I = 1,N B2 = B1 B1 = B0 NI = N + 1 - I B0 = TWOX*B1 - B2 + CS(NI) 10 CONTINUE C DCSEVL = 0.5D0*(B0-B2) C RETURN END *DECK DGAMLM SUBROUTINE DGAMLM (XMIN, XMAX) C***BEGIN PROLOGUE DGAMLM C***PURPOSE Compute the minimum and maximum bounds for the argument in C the Gamma function. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C7A, R2 C***TYPE DOUBLE PRECISION (GAMLIM-S, DGAMLM-D) C***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C Calculate the minimum and maximum legal bounds for X in gamma(X). C XMIN and XMAX are not the only bounds, but they are the only non- C trivial ones to calculate. C C Output Arguments -- C XMIN double precision minimum legal value of X in gamma(X). Any C smaller value of X might result in underflow. C XMAX double precision maximum legal value of X in gamma(X). Any C larger value of X might cause overflow. C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, XERMSG C***REVISION HISTORY (YYMMDD) C 770601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C***END PROLOGUE DGAMLM DOUBLE PRECISION XMIN, XMAX, ALNBIG, ALNSML, XLN, XOLD, D1MACH C***FIRST EXECUTABLE STATEMENT DGAMLM ALNSML = LOG(D1MACH(1)) XMIN = -ALNSML DO 10 I=1,10 XOLD = XMIN XLN = LOG(XMIN) XMIN = XMIN - XMIN*((XMIN+0.5D0)*XLN - XMIN - 0.2258D0 + ALNSML) 1 / (XMIN*XLN+0.5D0) IF (ABS(XMIN-XOLD).LT.0.005D0) GO TO 20 10 CONTINUE CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMIN', 1, 2) C 20 XMIN = -XMIN + 0.01D0 C ALNBIG = LOG (D1MACH(2)) XMAX = ALNBIG DO 30 I=1,10 XOLD = XMAX XLN = LOG(XMAX) XMAX = XMAX - XMAX*((XMAX-0.5D0)*XLN - XMAX + 0.9189D0 - ALNBIG) 1 / (XMAX*XLN-0.5D0) IF (ABS(XMAX-XOLD).LT.0.005D0) GO TO 40 30 CONTINUE CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMAX', 2, 2) C 40 XMAX = XMAX - 0.01D0 XMIN = MAX (XMIN, -XMAX+1.D0) C RETURN END *DECK DGAMMA DOUBLE PRECISION FUNCTION DGAMMA (X) C***BEGIN PROLOGUE DGAMMA C***PURPOSE Compute the complete Gamma function. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C7A C***TYPE DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C) C***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C DGAMMA(X) calculates the double precision complete Gamma function C for double precision argument X. C C Series for GAM on the interval 0. to 1.00000E+00 C with weighted error 5.79E-32 C log weighted error 31.24 C significant figures required 30.00 C decimal places required 32.05 C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG C***REVISION HISTORY (YYMMDD) C 770601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890911 Removed unnecessary intrinsics. (WRB) C 890911 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 920618 Removed space from variable name. (RWC, WRB) C***END PROLOGUE DGAMMA DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX, 1 XMIN, Y, D9LGMC, DCSEVL, D1MACH LOGICAL FIRST C SAVE GAMCS, PI, SQ2PIL, NGAM, XMIN, XMAX, DXREL, FIRST DATA GAMCS( 1) / +.8571195590 9893314219 2006239994 2 D-2 / DATA GAMCS( 2) / +.4415381324 8410067571 9131577165 2 D-2 / DATA GAMCS( 3) / +.5685043681 5993633786 3266458878 9 D-1 / DATA GAMCS( 4) / -.4219835396 4185605010 1250018662 4 D-2 / DATA GAMCS( 5) / +.1326808181 2124602205 8400679635 2 D-2 / DATA GAMCS( 6) / -.1893024529 7988804325 2394702388 6 D-3 / DATA GAMCS( 7) / +.3606925327 4412452565 7808221722 5 D-4 / DATA GAMCS( 8) / -.6056761904 4608642184 8554829036 5 D-5 / DATA GAMCS( 9) / +.1055829546 3022833447 3182350909 3 D-5 / DATA GAMCS( 10) / -.1811967365 5423840482 9185589116 6 D-6 / DATA GAMCS( 11) / +.3117724964 7153222777 9025459316 9 D-7 / DATA GAMCS( 12) / -.5354219639 0196871408 7408102434 7 D-8 / DATA GAMCS( 13) / +.9193275519 8595889468 8778682594 0 D-9 / DATA GAMCS( 14) / -.1577941280 2883397617 6742327395 3 D-9 / DATA GAMCS( 15) / +.2707980622 9349545432 6654043308 9 D-10 / DATA GAMCS( 16) / -.4646818653 8257301440 8166105893 3 D-11 / DATA GAMCS( 17) / +.7973350192 0074196564 6076717535 9 D-12 / DATA GAMCS( 18) / -.1368078209 8309160257 9949917230 9 D-12 / DATA GAMCS( 19) / +.2347319486 5638006572 3347177168 8 D-13 / DATA GAMCS( 20) / -.4027432614 9490669327 6657053469 9 D-14 / DATA GAMCS( 21) / +.6910051747 3721009121 3833697525 7 D-15 / DATA GAMCS( 22) / -.1185584500 2219929070 5238712619 2 D-15 / DATA GAMCS( 23) / +.2034148542 4963739552 0102605193 2 D-16 / DATA GAMCS( 24) / -.3490054341 7174058492 7401294910 8 D-17 / DATA GAMCS( 25) / +.5987993856 4853055671 3505106602 6 D-18 / DATA GAMCS( 26) / -.1027378057 8722280744 9006977843 1 D-18 / DATA GAMCS( 27) / +.1762702816 0605298249 4275966074 8 D-19 / DATA GAMCS( 28) / -.3024320653 7353062609 5877211204 2 D-20 / DATA GAMCS( 29) / +.5188914660 2183978397 1783355050 6 D-21 / DATA GAMCS( 30) / -.8902770842 4565766924 4925160106 6 D-22 / DATA GAMCS( 31) / +.1527474068 4933426022 7459689130 6 D-22 / DATA GAMCS( 32) / -.2620731256 1873629002 5732833279 9 D-23 / DATA GAMCS( 33) / +.4496464047 8305386703 3104657066 6 D-24 / DATA GAMCS( 34) / -.7714712731 3368779117 0390152533 3 D-25 / DATA GAMCS( 35) / +.1323635453 1260440364 8657271466 6 D-25 / DATA GAMCS( 36) / -.2270999412 9429288167 0231381333 3 D-26 / DATA GAMCS( 37) / +.3896418998 0039914493 2081663999 9 D-27 / DATA GAMCS( 38) / -.6685198115 1259533277 9212799999 9 D-28 / DATA GAMCS( 39) / +.1146998663 1400243843 4761386666 6 D-28 / DATA GAMCS( 40) / -.1967938586 3451346772 9510399999 9 D-29 / DATA GAMCS( 41) / +.3376448816 5853380903 3489066666 6 D-30 / DATA GAMCS( 42) / -.5793070335 7821357846 2549333333 3 D-31 / DATA PI / 3.1415926535 8979323846 2643383279 50 D0 / DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 / DATA FIRST /.TRUE./ C***FIRST EXECUTABLE STATEMENT DGAMMA IF (FIRST) THEN NGAM = INITDS (GAMCS, 42, 0.1*REAL(D1MACH(3)) ) C CALL DGAMLM (XMIN, XMAX) DXREL = SQRT(D1MACH(4)) ENDIF FIRST = .FALSE. C Y = ABS(X) IF (Y.GT.10.D0) GO TO 50 C C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND. REDUCE INTERVAL AND FIND C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL. C N = X IF (X.LT.0.D0) N = N - 1 Y = X - N N = N - 1 DGAMMA = 0.9375D0 + DCSEVL (2.D0*Y-1.D0, GAMCS, NGAM) IF (N.EQ.0) RETURN C IF (N.GT.0) GO TO 30 C C COMPUTE GAMMA(X) FOR X .LT. 1.0 C N = -N IF (X .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA', 'X IS 0', 4, 2) IF (X .LT. 0.0 .AND. X+N-2 .EQ. 0.D0) CALL XERMSG ('SLATEC', + 'DGAMMA', 'X IS A NEGATIVE INTEGER', 4, 2) IF (X .LT. (-0.5D0) .AND. ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) + CALL XERMSG ('SLATEC', 'DGAMMA', + 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER', + 1, 1) C DO 20 I=1,N DGAMMA = DGAMMA/(X+I-1 ) 20 CONTINUE RETURN C C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0 C 30 DO 40 I=1,N DGAMMA = (Y+I) * DGAMMA 40 CONTINUE RETURN C C GAMMA(X) FOR ABS(X) .GT. 10.0. RECALL Y = ABS(X). C 50 IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'DGAMMA', + 'X SO BIG GAMMA OVERFLOWS', 3, 2) C DGAMMA = 0.D0 IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'DGAMMA', + 'X SO SMALL GAMMA UNDERFLOWS', 2, 1) IF (X.LT.XMIN) RETURN C DGAMMA = EXP ((Y-0.5D0)*LOG(Y) - Y + SQ2PIL + D9LGMC(Y) ) IF (X.GT.0.D0) RETURN C IF (ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) CALL XERMSG ('SLATEC', + 'DGAMMA', + 'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1) C SINPIY = SIN (PI*Y) IF (SINPIY .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA', + 'X IS A NEGATIVE INTEGER', 4, 2) C DGAMMA = -PI/(Y*SINPIY*DGAMMA) C RETURN END *DECK INITDS FUNCTION INITDS (OS, NOS, ETA) C***BEGIN PROLOGUE INITDS C***PURPOSE Determine the number of terms needed in an orthogonal C polynomial series so that it meets a specified accuracy. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C3A2 C***TYPE DOUBLE PRECISION (INITS-S, INITDS-D) C***KEYWORDS CHEBYSHEV, FNLIB, INITIALIZE, ORTHOGONAL POLYNOMIAL, C ORTHOGONAL SERIES, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C Initialize the orthogonal series, represented by the array OS, so C that INITDS is the number of terms needed to insure the error is no C larger than ETA. Ordinarily, ETA will be chosen to be one-tenth C machine precision. C C Input Arguments -- C OS double precision array of NOS coefficients in an orthogonal C series. C NOS number of coefficients in OS. C ETA single precision scalar containing requested accuracy of C series. C C***REFERENCES (NONE) C***ROUTINES CALLED XERMSG C***REVISION HISTORY (YYMMDD) C 770601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 891115 Modified error message. (WRB) C 891115 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C***END PROLOGUE INITDS DOUBLE PRECISION OS(*) C***FIRST EXECUTABLE STATEMENT INITDS IF (NOS .LT. 1) CALL XERMSG ('SLATEC', 'INITDS', + 'Number of coefficients is less than 1', 2, 1) C ERR = 0. DO 10 II = 1,NOS I = NOS + 1 - II ERR = ERR + ABS(REAL(OS(I))) IF (ERR.GT.ETA) GO TO 20 10 CONTINUE C 20 IF (I .EQ. NOS) CALL XERMSG ('SLATEC', 'INITDS', + 'Chebyshev series too short for specified accuracy', 1, 1) INITDS = I C RETURN END *DECK D1MACH DOUBLE PRECISION FUNCTION D1MACH (I) C***BEGIN PROLOGUE D1MACH C***PURPOSE Return floating point machine dependent constants. C***LIBRARY SLATEC C***CATEGORY R1 C***TYPE DOUBLE PRECISION (R1MACH-S, D1MACH-D) C***KEYWORDS MACHINE CONSTANTS C***AUTHOR Fox, P. A., (Bell Labs) C Hall, A. D., (Bell Labs) C Schryer, N. L., (Bell Labs) C***DESCRIPTION C C D1MACH can be used to obtain machine-dependent parameters for the C local machine environment. It is a function subprogram with one C (input) argument, and can be referenced as follows: C C D = D1MACH(I) C C where I=1,...,5. The (output) value of D above is determined by C the (input) value of I. The results for various values of I are C discussed below. C C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C D1MACH( 3) = B**(-T), the smallest relative spacing. C D1MACH( 4) = B**(1-T), the largest relative spacing. C D1MACH( 5) = LOG10(B) C C Assume double precision numbers are represented in the T-digit, C base-B form C C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and C EMIN .LE. E .LE. EMAX. C C The values of B, T, EMIN and EMAX are provided in I1MACH as C follows: C I1MACH(10) = B, the base. C I1MACH(14) = T, the number of base-B digits. C I1MACH(15) = EMIN, the smallest exponent E. C I1MACH(16) = EMAX, the largest exponent E. C C To alter this function for a particular environment, the desired C set of DATA statements should be activated by removing the C from C column 1. Also, the values of D1MACH(1) - D1MACH(4) should be C checked for consistency with the local operating system. C C***REFERENCES P. A. Fox, A. D. Hall and N. L. Schryer, Framework for C a portable library, ACM Transactions on Mathematical C Software 4, 2 (June 1978), pp. 177-188. C***ROUTINES CALLED XERMSG C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 890213 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900618 Added DEC RISC constants. (WRB) C 900723 Added IBM RS 6000 constants. (WRB) C 900911 Added SUN 386i constants. (WRB) C 910710 Added HP 730 constants. (SMR) C 911114 Added Convex IEEE constants. (WRB) C 920121 Added SUN -r8 compiler option constants. (WRB) C 920229 Added Touchstone Delta i860 constants. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C 920625 Added CONVEX -p8 and -pd8 compiler option constants. C (BKS, WRB) C 930201 Added DEC Alpha and SGI constants. (RWC and WRB) C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) SAVE DMACH C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FDFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE APOLLO C C DATA SMALL(1), SMALL(2) / 16#00100000, 16#00000000 / C DATA LARGE(1), LARGE(2) / 16#7FFFFFFF, 16#FFFFFFFF / C DATA RIGHT(1), RIGHT(2) / 16#3CA00000, 16#00000000 / C DATA DIVER(1), DIVER(2) / 16#3CB00000, 16#00000000 / C DATA LOG10(1), LOG10(2) / 16#3FD34413, 16#509F79FF / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE C C DATA SMALL(1) / Z"3001800000000000" / C DATA SMALL(2) / Z"3001000000000000" / C DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" / C DATA LARGE(2) / Z"4FFE000000000000" / C DATA RIGHT(1) / Z"3FD2800000000000" / C DATA RIGHT(2) / Z"3FD2000000000000" / C DATA DIVER(1) / Z"3FD3800000000000" / C DATA DIVER(2) / Z"3FD3000000000000" / C DATA LOG10(1) / Z"3FFF9A209A84FBCF" / C DATA LOG10(2) / Z"3FFFF7988F8959AC" / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777777B / C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CELERITY C1260 C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -fn OR -pd8 COMPILER OPTION C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FFFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CC0000000000000' / C DATA DMACH(4) / Z'3CD0000000000000' / C DATA DMACH(5) / Z'3FF34413509F79FF' / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -fi COMPILER OPTION C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -p8 COMPILER OPTION C C DATA DMACH(1) / Z'00010000000000000000000000000000' / C DATA DMACH(2) / Z'7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3F900000000000000000000000000000' / C DATA DMACH(4) / Z'3F910000000000000000000000000000' / C DATA DMACH(5) / Z'3FFF34413509F79FEF311F12B35816F9' / C C MACHINE CONSTANTS FOR THE CRAY C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777774B / C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL / 20K, 3*0 / C DATA LARGE / 77777K, 3*177777K / C DATA RIGHT / 31420K, 3*0 / C DATA DIVER / 32020K, 3*0 / C DATA LOG10 / 40423K, 42023K, 50237K, 74776K / C C MACHINE CONSTANTS FOR THE DEC ALPHA C USING G_FLOAT C C DATA DMACH(1) / '0000000000000010'X / C DATA DMACH(2) / 'FFFFFFFFFFFF7FFF'X / C DATA DMACH(3) / '0000000000003CC0'X / C DATA DMACH(4) / '0000000000003CD0'X / C DATA DMACH(5) / '79FF509F44133FF3'X / C C MACHINE CONSTANTS FOR THE DEC ALPHA C USING IEEE_FORMAT C C DATA DMACH(1) / '0010000000000000'X / C DATA DMACH(2) / '7FEFFFFFFFFFFFFF'X / C DATA DMACH(3) / '3CA0000000000000'X / C DATA DMACH(4) / '3CB0000000000000'X / C DATA DMACH(5) / '3FD34413509F79FF'X / C C MACHINE CONSTANTS FOR THE DEC RISC C C DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000'/ C DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF'/ C DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000'/ C DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000'/ C DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413'/ C C MACHINE CONSTANTS FOR THE DEC VAX C USING D_FLOATING C (EXPRESSED IN INTEGER AND HEXADECIMAL) C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR THE DEC VAX C USING G_FLOATING C (EXPRESSED IN INTEGER AND HEXADECIMAL) C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / C DATA LOG10(1), LOG10(2) / '3FD34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / C DATA LARGE(1), LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1), DIVER(2) / '20000000, '00000334 / C DATA LOG10(1), LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES C C DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 730 C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C MACHINE CONSTANTS FOR THE HP 9000 C C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE IBM PC C ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION C ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087. C C DATA SMALL(1) / 2.23D-308 / C DATA LARGE(1) / 1.79D+308 / C DATA RIGHT(1) / 1.11D-16 / C DATA DIVER(1) / 2.22D-16 / C DATA LOG10(1) / 0.301029995663981195D0 / C C MACHINE CONSTANTS FOR THE IBM RS 6000 C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE INTEL i860 C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR) C C DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR) C C DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 8388608, 0 / C DATA LARGE(1), LARGE(2) / 2147483647, -1 / C DATA RIGHT(1), RIGHT(2) / 612368384, 0 / C DATA DIVER(1), DIVER(2) / 620756992, 0 / C DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA SMALL(3), SMALL(4) / 0, 0 / C DATA LARGE(1), LARGE(2) / 32767, -1 / C DATA LARGE(3), LARGE(4) / -1, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA DIVER(3), DIVER(4) / 0, 0 / C DATA LOG10(1), LOG10(2) / 16282, 8346 / C DATA LOG10(3), LOG10(4) / -31493, -12296 / C C DATA SMALL(1), SMALL(2) / O000200, O000000 / C DATA SMALL(3), SMALL(4) / O000000, O000000 / C DATA LARGE(1), LARGE(2) / O077777, O177777 / C DATA LARGE(3), LARGE(4) / O177777, O177777 / C DATA RIGHT(1), RIGHT(2) / O022200, O000000 / C DATA RIGHT(3), RIGHT(4) / O000000, O000000 / C DATA DIVER(1), DIVER(2) / O022400, O000000 / C DATA DIVER(3), DIVER(4) / O000000, O000000 / C DATA LOG10(1), LOG10(2) / O037632, O020232 / C DATA LOG10(3), LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE SILICON GRAPHICS C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE SUN C DATA DMACH(1) / 2.2250738585072D-308 / DATA DMACH(2) / 1.7976931348623D+308 / DATA DMACH(3) / 1.1102230246252D-16 / DATA DMACH(4) / 2.2204460492503D-16 / DATA DMACH(5) / 0.30102999566398D0 / C C MACHINE CONSTANTS FOR THE SUN C USING THE -r8 COMPILER OPTION C C DATA DMACH(1) / Z'00010000000000000000000000000000' / C DATA DMACH(2) / Z'7FFEFFFFFFFFFFFFFFFFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3F8E0000000000000000000000000000' / C DATA DMACH(4) / Z'3F8F0000000000000000000000000000' / C DATA DMACH(5) / Z'3FFD34413509F79FEF311F12B35816F9' / C C MACHINE CONSTANTS FOR THE SUN 386i C C DATA SMALL(1), SMALL(2) / Z'FFFFFFFD', Z'000FFFFF' / C DATA LARGE(1), LARGE(2) / Z'FFFFFFB0', Z'7FEFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'000000B0', Z'3CA00000' / C DATA DIVER(1), DIVER(2) / Z'FFFFFFCB', Z'3CAFFFFF' C DATA LOG10(1), LOG10(2) / Z'509F79E9', Z'3FD34413' / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER C C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 / C C***FIRST EXECUTABLE STATEMENT D1MACH IF (I .LT. 1 .OR. I .GT. 5) CALL XERMSG ('SLATEC', 'D1MACH', + 'I OUT OF BOUNDS', 1, 2) C D1MACH = DMACH(I) RETURN C END *DECK FDUMP SUBROUTINE FDUMP C***BEGIN PROLOGUE FDUMP C***PURPOSE Symbolic dump (should be locally written). C***LIBRARY SLATEC (XERROR) C***CATEGORY R3 C***TYPE ALL (FDUMP-A) C***KEYWORDS ERROR, XERMSG C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C ***Note*** Machine Dependent Routine C FDUMP is intended to be replaced by a locally written C version which produces a symbolic dump. Failing this, C it should be replaced by a version which prints the C subprogram nesting list. Note that this dump must be C printed on each of up to five files, as indicated by the C XGETUA routine. See XSETUA and XGETUA for details. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C C***REFERENCES (NONE) C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 790801 DATE WRITTEN C 861211 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE FDUMP C***FIRST EXECUTABLE STATEMENT FDUMP RETURN END *DECK I1MACH INTEGER FUNCTION I1MACH (I) C***BEGIN PROLOGUE I1MACH C***PURPOSE Return integer machine dependent constants. C***LIBRARY SLATEC C***CATEGORY R1 C***TYPE INTEGER (I1MACH-I) C***KEYWORDS MACHINE CONSTANTS C***AUTHOR Fox, P. A., (Bell Labs) C Hall, A. D., (Bell Labs) C Schryer, N. L., (Bell Labs) C***DESCRIPTION C C I1MACH can be used to obtain machine-dependent parameters for the C local machine environment. It is a function subprogram with one C (input) argument and can be referenced as follows: C C K = I1MACH(I) C C where I=1,...,16. The (output) value of K above is determined by C the (input) value of I. The results for various values of I are C discussed below. C C I/O unit numbers: C I1MACH( 1) = the standard input unit. C I1MACH( 2) = the standard output unit. C I1MACH( 3) = the standard punch unit. C I1MACH( 4) = the standard error message unit. C C Words: C I1MACH( 5) = the number of bits per integer storage unit. C I1MACH( 6) = the number of characters per integer storage unit. C C Integers: C assume integers are represented in the S-digit, base-A form C C sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C where 0 .LE. X(I) .LT. A for I=0,...,S-1. C I1MACH( 7) = A, the base. C I1MACH( 8) = S, the number of base-A digits. C I1MACH( 9) = A**S - 1, the largest magnitude. C C Floating-Point Numbers: C Assume floating-point numbers are represented in the T-digit, C base-B form C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C where 0 .LE. X(I) .LT. B for I=1,...,T, C 0 .LT. X(1), and EMIN .LE. E .LE. EMAX. C I1MACH(10) = B, the base. C C Single-Precision: C I1MACH(11) = T, the number of base-B digits. C I1MACH(12) = EMIN, the smallest exponent E. C I1MACH(13) = EMAX, the largest exponent E. C C Double-Precision: C I1MACH(14) = T, the number of base-B digits. C I1MACH(15) = EMIN, the smallest exponent E. C I1MACH(16) = EMAX, the largest exponent E. C C To alter this function for a particular environment, the desired C set of DATA statements should be activated by removing the C from C column 1. Also, the values of I1MACH(1) - I1MACH(4) should be C checked for consistency with the local operating system. C C***REFERENCES P. A. Fox, A. D. Hall and N. L. Schryer, Framework for C a portable library, ACM Transactions on Mathematical C Software 4, 2 (June 1978), pp. 177-188. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 891012 Added VAX G-floating constants. (WRB) C 891012 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900618 Added DEC RISC constants. (WRB) C 900723 Added IBM RS 6000 constants. (WRB) C 901009 Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16. C (RWC) C 910710 Added HP 730 constants. (SMR) C 911114 Added Convex IEEE constants. (WRB) C 920121 Added SUN -r8 compiler option constants. (WRB) C 920229 Added Touchstone Delta i860 constants. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C 920625 Added Convex -p8 and -pd8 compiler option constants. C (BKS, WRB) C 930201 Added DEC Alpha and SGI constants. (RWC and WRB) C 930618 Corrected I1MACH(5) for Convex -p8 and -pd8 compiler C options. (DWL, RWC and WRB). C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT SAVE IMACH EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT COMPILER C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE APOLLO C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 129 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1025 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM C C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 / C C MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 9223372036854775807 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -4095 / C DATA IMACH(13) / 4094 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -4095 / C DATA IMACH(16) / 4094 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6LOUTPUT/ C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -929 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -929 / C DATA IMACH(16) / 1069 / C C MACHINE CONSTANTS FOR THE CELERITY C1260 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z'7FFFFFFF' / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -fn COMPILER OPTION C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -fi COMPILER OPTION C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -p8 COMPILER OPTION C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 9223372036854775807 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1023 / C DATA IMACH(13) / 1023 / C DATA IMACH(14) / 113 / C DATA IMACH(15) / -16383 / C DATA IMACH(16) / 16383 / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -pd8 COMPILER OPTION C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 9223372036854775807 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1023 / C DATA IMACH(13) / 1023 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CRAY C USING THE 46 BIT INTEGER COMPILER OPTION C C DATA IMACH( 1) / 100 / C DATA IMACH( 2) / 101 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 101 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 46 / C DATA IMACH( 9) / 1777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 / C C MACHINE CONSTANTS FOR THE CRAY C USING THE 64 BIT INTEGER COMPILER OPTION C C DATA IMACH( 1) / 100 / C DATA IMACH( 2) / 101 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 101 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE DEC ALPHA C USING G_FLOAT C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE DEC ALPHA C USING IEEE_FLOAT C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE DEC RISC C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE DEC VAX C USING D_FLOATING C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE DEC VAX C USING G_FLOATING C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 730 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 4 / C DATA IMACH( 4) / 1 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 39 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 4 / C DATA IMACH( 4) / 1 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 55 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 9000 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 7 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1015 / C DATA IMACH(16) / 1017 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE IBM PC C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE IBM RS 6000 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE INTEL i860 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR) C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR) C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE SILICON GRAPHICS C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE SUN C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE SUN C USING THE -r8 COMPILER OPTION C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1021 / C DATA IMACH(13) / 1024 / C DATA IMACH(14) / 113 / C DATA IMACH(15) / -16381 / C DATA IMACH(16) / 16384 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 1 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -1024 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR C C DATA IMACH( 1) / 1 / C DATA IMACH( 2) / 1 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 1 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C***FIRST EXECUTABLE STATEMENT I1MACH IF (I .LT. 1 .OR. I .GT. 16) GO TO 10 C I1MACH = IMACH(I) RETURN C 10 CONTINUE WRITE (UNIT = OUTPUT, FMT = 9000) 9000 FORMAT ('1ERROR 1 IN I1MACH - I OUT OF BOUNDS') C C CALL FDUMP C STOP END *DECK J4SAVE FUNCTION J4SAVE (IWHICH, IVALUE, ISET) C***BEGIN PROLOGUE J4SAVE C***SUBSIDIARY C***PURPOSE Save or recall global variables needed by error C handling routines. C***LIBRARY SLATEC (XERROR) C***TYPE INTEGER (J4SAVE-I) C***KEYWORDS ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C Abstract C J4SAVE saves and recalls several global variables needed C by the library error handling routines. C C Description of Parameters C --Input-- C IWHICH - Index of item desired. C = 1 Refers to current error number. C = 2 Refers to current error control flag. C = 3 Refers to current unit number to which error C messages are to be sent. (0 means use standard.) C = 4 Refers to the maximum number of times any C message is to be printed (as set by XERMAX). C = 5 Refers to the total number of units to which C each error message is to be written. C = 6 Refers to the 2nd unit for error messages C = 7 Refers to the 3rd unit for error messages C = 8 Refers to the 4th unit for error messages C = 9 Refers to the 5th unit for error messages C IVALUE - The value to be set for the IWHICH-th parameter, C if ISET is .TRUE. . C ISET - If ISET=.TRUE., the IWHICH-th parameter will BE C given the value, IVALUE. If ISET=.FALSE., the C IWHICH-th parameter will be unchanged, and IVALUE C is a dummy parameter. C --Output-- C The (old) value of the IWHICH-th parameter will be returned C in the function value, J4SAVE. C C***SEE ALSO XERMSG C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC C Error-handling Package, SAND82-0800, Sandia C Laboratories, 1982. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 790801 DATE WRITTEN C 891214 Prologue converted to Version 4.0 format. (BAB) C 900205 Minor modifications to prologue. (WRB) C 900402 Added TYPE section. (WRB) C 910411 Added KEYWORDS section. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE J4SAVE LOGICAL ISET INTEGER IPARAM(9) SAVE IPARAM DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/ DATA IPARAM(5)/1/ DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/ C***FIRST EXECUTABLE STATEMENT J4SAVE J4SAVE = IPARAM(IWHICH) IF (ISET) IPARAM(IWHICH) = IVALUE RETURN END *DECK XERCNT SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL) C***BEGIN PROLOGUE XERCNT C***SUBSIDIARY C***PURPOSE Allow user control over handling of errors. C***LIBRARY SLATEC (XERROR) C***CATEGORY R3C C***TYPE ALL (XERCNT-A) C***KEYWORDS ERROR, XERROR C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C Abstract C Allows user control over handling of individual errors. C Just after each message is recorded, but before it is C processed any further (i.e., before it is printed or C a decision to abort is made), a call is made to XERCNT. C If the user has provided his own version of XERCNT, he C can then override the value of KONTROL used in processing C this message by redefining its value. C KONTRL may be set to any value from -2 to 2. C The meanings for KONTRL are the same as in XSETF, except C that the value of KONTRL changes only for this message. C If KONTRL is set to a value outside the range from -2 to 2, C it will be moved back into that range. C C Description of Parameters C C --Input-- C LIBRAR - the library that the routine is in. C SUBROU - the subroutine that XERMSG is being called from C MESSG - the first 20 characters of the error message. C NERR - same as in the call to XERMSG. C LEVEL - same as in the call to XERMSG. C KONTRL - the current value of the control flag as set C by a call to XSETF. C C --Output-- C KONTRL - the new value of KONTRL. If KONTRL is not C defined, it will remain at its original value. C This changed value of control affects only C the current occurrence of the current message. C C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC C Error-handling Package, SAND82-0800, Sandia C Laboratories, 1982. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 790801 DATE WRITTEN C 861211 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900206 Routine changed from user-callable to subsidiary. (WRB) C 900510 Changed calling sequence to include LIBRARY and SUBROUTINE C names, changed routine name from XERCTL to XERCNT. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE XERCNT CHARACTER*(*) LIBRAR, SUBROU, MESSG C***FIRST EXECUTABLE STATEMENT XERCNT RETURN END *DECK XERHLT SUBROUTINE XERHLT (MESSG) C***BEGIN PROLOGUE XERHLT C***SUBSIDIARY C***PURPOSE Abort program execution and print error message. C***LIBRARY SLATEC (XERROR) C***CATEGORY R3C C***TYPE ALL (XERHLT-A) C***KEYWORDS ABORT PROGRAM EXECUTION, ERROR, XERROR C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C Abstract C ***Note*** machine dependent routine C XERHLT aborts the execution of the program. C The error message causing the abort is given in the calling C sequence, in case one needs it for printing on a dayfile, C for example. C C Description of Parameters C MESSG is as in XERMSG. C C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC C Error-handling Package, SAND82-0800, Sandia C Laboratories, 1982. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 790801 DATE WRITTEN C 861211 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900206 Routine changed from user-callable to subsidiary. (WRB) C 900510 Changed calling sequence to delete length of character C and changed routine name from XERABT to XERHLT. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE XERHLT CHARACTER*(*) MESSG C***FIRST EXECUTABLE STATEMENT XERHLT STOP END *DECK XERMSG SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL) C***BEGIN PROLOGUE XERMSG C***PURPOSE Process error messages for SLATEC and other libraries. C***LIBRARY SLATEC (XERROR) C***CATEGORY R3C C***TYPE ALL (XERMSG-A) C***KEYWORDS ERROR MESSAGE, XERROR C***AUTHOR Fong, Kirby, (NMFECC at LLNL) C***DESCRIPTION C C XERMSG processes a diagnostic message in a manner determined by the C value of LEVEL and the current value of the library error control C flag, KONTRL. See subroutine XSETF for details. C C LIBRAR A character constant (or character variable) with the name C of the library. This will be 'SLATEC' for the SLATEC C Common Math Library. The error handling package is C general enough to be used by many libraries C simultaneously, so it is desirable for the routine that C detects and reports an error to identify the library name C as well as the routine name. C C SUBROU A character constant (or character variable) with the name C of the routine that detected the error. Usually it is the C name of the routine that is calling XERMSG. There are C some instances where a user callable library routine calls C lower level subsidiary routines where the error is C detected. In such cases it may be more informative to C supply the name of the routine the user called rather than C the name of the subsidiary routine that detected the C error. C C MESSG A character constant (or character variable) with the text C of the error or warning message. In the example below, C the message is a character constant that contains a C generic message. C C CALL XERMSG ('SLATEC', 'MMPY', C *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION', C *3, 1) C C It is possible (and is sometimes desirable) to generate a C specific message--e.g., one that contains actual numeric C values. Specific numeric values can be converted into C character strings using formatted WRITE statements into C character variables. This is called standard Fortran C internal file I/O and is exemplified in the first three C lines of the following example. You can also catenate C substrings of characters to construct the error message. C Here is an example showing the use of both writing to C an internal file and catenating character strings. C C CHARACTER*5 CHARN, CHARL C WRITE (CHARN,10) N C WRITE (CHARL,10) LDA C 10 FORMAT(I5) C CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN// C * ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'// C * CHARL, 3, 1) C C There are two subtleties worth mentioning. One is that C the // for character catenation is used to construct the C error message so that no single character constant is C continued to the next line. This avoids confusion as to C whether there are trailing blanks at the end of the line. C The second is that by catenating the parts of the message C as an actual argument rather than encoding the entire C message into one large character variable, we avoid C having to know how long the message will be in order to C declare an adequate length for that large character C variable. XERMSG calls XERPRN to print the message using C multiple lines if necessary. If the message is very long, C XERPRN will break it into pieces of 72 characters (as C requested by XERMSG) for printing on multiple lines. C Also, XERMSG asks XERPRN to prefix each line with ' * ' C so that the total line length could be 76 characters. C Note also that XERPRN scans the error message backwards C to ignore trailing blanks. Another feature is that C the substring '$$' is treated as a new line sentinel C by XERPRN. If you want to construct a multiline C message without having to count out multiples of 72 C characters, just use '$$' as a separator. '$$' C obviously must occur within 72 characters of the C start of each line to have its intended effect since C XERPRN is asked to wrap around at 72 characters in C addition to looking for '$$'. C C NERR An integer value that is chosen by the library routine's C author. It must be in the range -99 to 999 (three C printable digits). Each distinct error should have its C own error number. These error numbers should be described C in the machine readable documentation for the routine. C The error numbers need be unique only within each routine, C so it is reasonable for each routine to start enumerating C errors from 1 and proceeding to the next integer. C C LEVEL An integer value in the range 0 to 2 that indicates the C level (severity) of the error. Their meanings are C C -1 A warning message. This is used if it is not clear C that there really is an error, but the user's attention C may be needed. An attempt is made to only print this C message once. C C 0 A warning message. This is used if it is not clear C that there really is an error, but the user's attention C may be needed. C C 1 A recoverable error. This is used even if the error is C so serious that the routine cannot return any useful C answer. If the user has told the error package to C return after recoverable errors, then XERMSG will C return to the Library routine which can then return to C the user's routine. The user may also permit the error C package to terminate the program upon encountering a C recoverable error. C C 2 A fatal error. XERMSG will not return to its caller C after it receives a fatal error. This level should C hardly ever be used; it is much better to allow the C user a chance to recover. An example of one of the few C cases in which it is permissible to declare a level 2 C error is a reverse communication Library routine that C is likely to be called repeatedly until it integrates C across some interval. If there is a serious error in C the input such that another step cannot be taken and C the Library routine is called again without the input C error having been corrected by the caller, the Library C routine will probably be called forever with improper C input. In this case, it is reasonable to declare the C error to be fatal. C C Each of the arguments to XERMSG is input; none will be modified by C XERMSG. A routine may make multiple calls to XERMSG with warning C level messages; however, after a call to XERMSG with a recoverable C error, the routine should return to the user. Do not try to call C XERMSG with a second recoverable error after the first recoverable C error because the error package saves the error number. The user C can retrieve this error number by calling another entry point in C the error handling package and then clear the error number when C recovering from the error. Calling XERMSG in succession causes the C old error number to be overwritten by the latest error number. C This is considered harmless for error numbers associated with C warning messages but must not be done for error numbers of serious C errors. After a call to XERMSG with a recoverable error, the user C must be given a chance to call NUMXER or XERCLR to retrieve or C clear the error number. C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC C Error-handling Package, SAND82-0800, Sandia C Laboratories, 1982. C***ROUTINES CALLED FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE C***REVISION HISTORY (YYMMDD) C 880101 DATE WRITTEN C 880621 REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988. C THERE ARE TWO BASIC CHANGES. C 1. A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO C PRINT MESSAGES. THIS ROUTINE WILL BREAK LONG MESSAGES C INTO PIECES FOR PRINTING ON MULTIPLE LINES. '$$' IS C ACCEPTED AS A NEW LINE SENTINEL. A PREFIX CAN BE C ADDED TO EACH LINE TO BE PRINTED. XERMSG USES EITHER C ' ***' OR ' * ' AND LONG MESSAGES ARE BROKEN EVERY C 72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE C LENGTH OUTPUT CAN NOW BE AS GREAT AS 76. C 2. THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE C FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE C OF LOWER CASE. C 880708 REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30. C THE PRINCIPAL CHANGES ARE C 1. CLARIFY COMMENTS IN THE PROLOGUES C 2. RENAME XRPRNT TO XERPRN C 3. REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES C SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE / C CHARACTER FOR NEW RECORDS. C 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO C CLEAN UP THE CODING. C 890721 REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN C PREFIX. C 891013 REVISED TO CORRECT COMMENTS. C 891214 Prologue converted to Version 4.0 format. (WRB) C 900510 Changed test on NERR to be -9999999 < NERR < 99999999, but C NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3. Added C LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and C XERCTL to XERCNT. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE XERMSG CHARACTER*(*) LIBRAR, SUBROU, MESSG CHARACTER*8 XLIBR, XSUBR CHARACTER*72 TEMP CHARACTER*20 LFIRST C***FIRST EXECUTABLE STATEMENT XERMSG LKNTRL = J4SAVE (2, 0, .FALSE.) MAXMES = J4SAVE (4, 0, .FALSE.) C C LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL. C MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE C SHOULD BE PRINTED. C C WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN C CALLING XERMSG. THE ERROR NUMBER SHOULD BE POSITIVE, C AND THE LEVEL SHOULD BE BETWEEN 0 AND 2. C IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR. * LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' // * 'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '// * 'JOB ABORT DUE TO FATAL ERROR.', 72) CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY) CALL XERHLT (' ***XERMSG -- INVALID INPUT') RETURN ENDIF C C RECORD THE MESSAGE. C I = J4SAVE (1, NERR, .TRUE.) CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT) C C HANDLE PRINT-ONCE WARNING MESSAGES. C IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN C C ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG. C XLIBR = LIBRAR XSUBR = SUBROU LFIRST = MESSG LERR = NERR LLEVEL = LEVEL CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL) C LKNTRL = MAX(-2, MIN(2,LKNTRL)) MKNTRL = ABS(LKNTRL) C C SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS C ZERO AND THE ERROR IS NOT FATAL. C IF (LEVEL.LT.2 .AND. LKNTRL.EQ.0) GO TO 30 IF (LEVEL.EQ.0 .AND. KOUNT.GT.MAXMES) GO TO 30 IF (LEVEL.EQ.1 .AND. KOUNT.GT.MAXMES .AND. MKNTRL.EQ.1) GO TO 30 IF (LEVEL.EQ.2 .AND. KOUNT.GT.MAX(1,MAXMES)) GO TO 30 C C ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A C MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS) C AND SENDING IT OUT VIA XERPRN. PRINT ONLY IF CONTROL FLAG C IS NOT ZERO. C IF (LKNTRL .NE. 0) THEN TEMP(1:21) = 'MESSAGE FROM ROUTINE ' I = MIN(LEN(SUBROU), 16) TEMP(22:21+I) = SUBROU(1:I) TEMP(22+I:33+I) = ' IN LIBRARY ' LTEMP = 33 + I I = MIN(LEN(LIBRAR), 16) TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I) TEMP(LTEMP+I+1:LTEMP+I+1) = '.' LTEMP = LTEMP + I + 1 CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72) ENDIF C C IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE C PRINTING THE MESSAGE. THE INTRODUCTORY LINE TELLS THE CHOICE C FROM EACH OF THE FOLLOWING THREE OPTIONS. C 1. LEVEL OF THE MESSAGE C 'INFORMATIVE MESSAGE' C 'POTENTIALLY RECOVERABLE ERROR' C 'FATAL ERROR' C 2. WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE C 'PROG CONTINUES' C 'PROG ABORTED' C 3. WHETHER OR NOT A TRACEBACK WAS REQUESTED. (THE TRACEBACK C MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS C WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.) C 'TRACEBACK REQUESTED' C 'TRACEBACK NOT REQUESTED' C NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT C EXCEED 74 CHARACTERS. C WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED. C IF (LKNTRL .GT. 0) THEN C C THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL. C IF (LEVEL .LE. 0) THEN TEMP(1:20) = 'INFORMATIVE MESSAGE,' LTEMP = 20 ELSEIF (LEVEL .EQ. 1) THEN TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,' LTEMP = 30 ELSE TEMP(1:12) = 'FATAL ERROR,' LTEMP = 12 ENDIF C C THEN WHETHER THE PROGRAM WILL CONTINUE. C IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR. * (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN TEMP(LTEMP+1:LTEMP+14) = ' PROG ABORTED,' LTEMP = LTEMP + 14 ELSE TEMP(LTEMP+1:LTEMP+16) = ' PROG CONTINUES,' LTEMP = LTEMP + 16 ENDIF C C FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK. C IF (LKNTRL .GT. 0) THEN TEMP(LTEMP+1:LTEMP+20) = ' TRACEBACK REQUESTED' LTEMP = LTEMP + 20 ELSE TEMP(LTEMP+1:LTEMP+24) = ' TRACEBACK NOT REQUESTED' LTEMP = LTEMP + 24 ENDIF CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72) ENDIF C C NOW SEND OUT THE MESSAGE. C CALL XERPRN (' * ', -1, MESSG, 72) C C IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A C TRACEBACK. C IF (LKNTRL .GT. 0) THEN WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR DO 10 I=16,22 IF (TEMP(I:I) .NE. ' ') GO TO 20 10 CONTINUE C 20 CALL XERPRN (' * ', -1, TEMP(1:15) // TEMP(I:23), 72) CALL FDUMP ENDIF C C IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE. C IF (LKNTRL .NE. 0) THEN CALL XERPRN (' * ', -1, ' ', 72) CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72) CALL XERPRN (' ', 0, ' ', 72) ENDIF C C IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE C CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN. C 30 IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN C C THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A C FATAL ERROR. PRINT THE REASON FOR THE ABORT AND THE ERROR C SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT. C IF (LKNTRL.GT.0 .AND. KOUNT.LT.MAX(1,MAXMES)) THEN IF (LEVEL .EQ. 1) THEN CALL XERPRN * (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72) ELSE CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72) ENDIF CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY) CALL XERHLT (' ') ELSE CALL XERHLT (MESSG) ENDIF RETURN END *DECK XERPRN SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP) C***BEGIN PROLOGUE XERPRN C***SUBSIDIARY C***PURPOSE Print error messages processed by XERMSG. C***LIBRARY SLATEC (XERROR) C***CATEGORY R3C C***TYPE ALL (XERPRN-A) C***KEYWORDS ERROR MESSAGES, PRINTING, XERROR C***AUTHOR Fong, Kirby, (NMFECC at LLNL) C***DESCRIPTION C C This routine sends one or more lines to each of the (up to five) C logical units to which error messages are to be sent. This routine C is called several times by XERMSG, sometimes with a single line to C print and sometimes with a (potentially very long) message that may C wrap around into multiple lines. C C PREFIX Input argument of type CHARACTER. This argument contains C characters to be put at the beginning of each line before C the body of the message. No more than 16 characters of C PREFIX will be used. C C NPREF Input argument of type INTEGER. This argument is the number C of characters to use from PREFIX. If it is negative, the C intrinsic function LEN is used to determine its length. If C it is zero, PREFIX is not used. If it exceeds 16 or if C LEN(PREFIX) exceeds 16, only the first 16 characters will be C used. If NPREF is positive and the length of PREFIX is less C than NPREF, a copy of PREFIX extended with blanks to length C NPREF will be used. C C MESSG Input argument of type CHARACTER. This is the text of a C message to be printed. If it is a long message, it will be C broken into pieces for printing on multiple lines. Each line C will start with the appropriate prefix and be followed by a C piece of the message. NWRAP is the number of characters per C piece; that is, after each NWRAP characters, we break and C start a new line. In addition the characters '$$' embedded C in MESSG are a sentinel for a new line. The counting of C characters up to NWRAP starts over for each new line. The C value of NWRAP typically used by XERMSG is 72 since many C older error messages in the SLATEC Library are laid out to C rely on wrap-around every 72 characters. C C NWRAP Input argument of type INTEGER. This gives the maximum size C piece into which to break MESSG for printing on multiple C lines. An embedded '$$' ends a line, and the count restarts C at the following character. If a line break does not occur C on a blank (it would split a word) that word is moved to the C next line. Values of NWRAP less than 16 will be treated as C 16. Values of NWRAP greater than 132 will be treated as 132. C The actual line length will be NPREF + NWRAP after NPREF has C been adjusted to fall between 0 and 16 and NWRAP has been C adjusted to fall between 16 and 132. C C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC C Error-handling Package, SAND82-0800, Sandia C Laboratories, 1982. C***ROUTINES CALLED I1MACH, XGETUA C***REVISION HISTORY (YYMMDD) C 880621 DATE WRITTEN C 880708 REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF C JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK C THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE C SLASH CHARACTER IN FORMAT STATEMENTS. C 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO C STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK C LINES TO BE PRINTED. C 890721 REVISED TO ADD A NEW FEATURE. A NEGATIVE VALUE OF NPREF C CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH. C 891013 REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH. C 891214 Prologue converted to Version 4.0 format. (WRB) C 900510 Added code to break messages between words. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE XERPRN CHARACTER*(*) PREFIX, MESSG INTEGER NPREF, NWRAP CHARACTER*148 CBUFF INTEGER IU(5), NUNIT CHARACTER*2 NEWLIN PARAMETER (NEWLIN = '$$') C***FIRST EXECUTABLE STATEMENT XERPRN CALL XGETUA(IU,NUNIT) C C A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD C ERROR MESSAGE UNIT INSTEAD. I1MACH(4) RETRIEVES THE STANDARD C ERROR MESSAGE UNIT. C N = I1MACH(4) DO 10 I=1,NUNIT IF (IU(I) .EQ. 0) IU(I) = N 10 CONTINUE C C LPREF IS THE LENGTH OF THE PREFIX. THE PREFIX IS PLACED AT THE C BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING C THE REST OF THIS ROUTINE. C IF ( NPREF .LT. 0 ) THEN LPREF = LEN(PREFIX) ELSE LPREF = NPREF ENDIF LPREF = MIN(16, LPREF) IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX C C LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE C TIME FROM MESSG TO PRINT ON ONE LINE. C LWRAP = MAX(16, MIN(132, NWRAP)) C C SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS. C LENMSG = LEN(MESSG) N = LENMSG DO 20 I=1,N IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30 LENMSG = LENMSG - 1 20 CONTINUE 30 CONTINUE C C IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE. C IF (LENMSG .EQ. 0) THEN CBUFF(LPREF+1:LPREF+1) = ' ' DO 40 I=1,NUNIT WRITE(IU(I), '(A)') CBUFF(1:LPREF+1) 40 CONTINUE RETURN ENDIF C C SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING C STARTS. FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL. C WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT. C WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED. C C WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL. THE C INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE C OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH C OF THE SECOND ARGUMENT. C C THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE C FOLLOWING ORDER. WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER C OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT C POSITION NEXTC. C C LPIECE .EQ. 0 THE NEW LINE SENTINEL DOES NOT OCCUR IN THE C REMAINDER OF THE CHARACTER STRING. LPIECE C SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC, C WHICHEVER IS LESS. C C LPIECE .EQ. 1 THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC: C NEXTC). LPIECE IS EFFECTIVELY ZERO, AND WE C PRINT NOTHING TO AVOID PRODUCING UNNECESSARY C BLANK LINES. THIS TAKES CARE OF THE SITUATION C WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF C EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE C SENTINEL FOLLOWED BY MORE CHARACTERS. NEXTC C SHOULD BE INCREMENTED BY 2. C C LPIECE .GT. LWRAP+1 REDUCE LPIECE TO LWRAP. C C ELSE THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1 C RESET LPIECE = LPIECE-1. NOTE THAT THIS C PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ. C LWRAP+1. THAT IS, THE SENTINEL FALLS EXACTLY C AT THE END OF A LINE. C NEXTC = 1 50 LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN) IF (LPIECE .EQ. 0) THEN C C THERE WAS NO NEW LINE SENTINEL FOUND. C IDELTA = 0 LPIECE = MIN(LWRAP, LENMSG+1-NEXTC) IF (LPIECE .LT. LENMSG+1-NEXTC) THEN DO 52 I=LPIECE+1,2,-1 IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN LPIECE = I-1 IDELTA = 1 GOTO 54 ENDIF 52 CONTINUE ENDIF 54 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1) NEXTC = NEXTC + LPIECE + IDELTA ELSEIF (LPIECE .EQ. 1) THEN C C WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1). C DON'T PRINT A BLANK LINE. C NEXTC = NEXTC + 2 GO TO 50 ELSEIF (LPIECE .GT. LWRAP+1) THEN C C LPIECE SHOULD BE SET DOWN TO LWRAP. C IDELTA = 0 LPIECE = LWRAP DO 56 I=LPIECE+1,2,-1 IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN LPIECE = I-1 IDELTA = 1 GOTO 58 ENDIF 56 CONTINUE 58 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1) NEXTC = NEXTC + LPIECE + IDELTA ELSE C C IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1. C WE SHOULD DECREMENT LPIECE BY ONE. C LPIECE = LPIECE - 1 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1) NEXTC = NEXTC + LPIECE + 2 ENDIF C C PRINT C DO 60 I=1,NUNIT WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE) 60 CONTINUE C IF (NEXTC .LE. LENMSG) GO TO 50 RETURN END *DECK XERSVE SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, + ICOUNT) C***BEGIN PROLOGUE XERSVE C***SUBSIDIARY C***PURPOSE Record that an error has occurred. C***LIBRARY SLATEC (XERROR) C***CATEGORY R3 C***TYPE ALL (XERSVE-A) C***KEYWORDS ERROR, XERROR C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C *Usage: C C INTEGER KFLAG, NERR, LEVEL, ICOUNT C CHARACTER * (len) LIBRAR, SUBROU, MESSG C C CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT) C C *Arguments: C C LIBRAR :IN is the library that the message is from. C SUBROU :IN is the subroutine that the message is from. C MESSG :IN is the message to be saved. C KFLAG :IN indicates the action to be performed. C when KFLAG > 0, the message in MESSG is saved. C when KFLAG=0 the tables will be dumped and C cleared. C when KFLAG < 0, the tables will be dumped and C not cleared. C NERR :IN is the error number. C LEVEL :IN is the error severity. C ICOUNT :OUT the number of times this message has been seen, C or zero if the table has overflowed and does not C contain this message specifically. When KFLAG=0, C ICOUNT will not be altered. C C *Description: C C Record that this error occurred and possibly dump and clear the C tables. C C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC C Error-handling Package, SAND82-0800, Sandia C Laboratories, 1982. C***ROUTINES CALLED I1MACH, XGETUA C***REVISION HISTORY (YYMMDD) C 800319 DATE WRITTEN C 861211 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900413 Routine modified to remove reference to KFLAG. (WRB) C 900510 Changed to add LIBRARY NAME and SUBROUTINE to calling C sequence, use IF-THEN-ELSE, make number of saved entries C easily changeable, changed routine name from XERSAV to C XERSVE. (RWC) C 910626 Added LIBTAB and SUBTAB to SAVE statement. (BKS) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE XERSVE PARAMETER (LENTAB=10) INTEGER LUN(5) CHARACTER*(*) LIBRAR, SUBROU, MESSG CHARACTER*8 LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB CHARACTER*20 MESTAB(LENTAB), MES DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB) SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG DATA KOUNTX/0/, NMSG/0/ C***FIRST EXECUTABLE STATEMENT XERSVE C IF (KFLAG.LE.0) THEN C C Dump the table. C IF (NMSG.EQ.0) RETURN C C Print to each unit. C CALL XGETUA (LUN, NUNIT) DO 20 KUNIT = 1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) C C Print the table header. C WRITE (IUNIT,9000) C C Print body of table. C DO 10 I = 1,NMSG WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I), * NERTAB(I),LEVTAB(I),KOUNT(I) 10 CONTINUE C C Print number of other errors. C IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX WRITE (IUNIT,9030) 20 CONTINUE C C Clear the error tables. C IF (KFLAG.EQ.0) THEN NMSG = 0 KOUNTX = 0 ENDIF ELSE C C PROCESS A MESSAGE... C SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, C OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. C LIB = LIBRAR SUB = SUBROU MES = MESSG DO 30 I = 1,NMSG IF (LIB.EQ.LIBTAB(I) .AND. SUB.EQ.SUBTAB(I) .AND. * MES.EQ.MESTAB(I) .AND. NERR.EQ.NERTAB(I) .AND. * LEVEL.EQ.LEVTAB(I)) THEN KOUNT(I) = KOUNT(I) + 1 ICOUNT = KOUNT(I) RETURN ENDIF 30 CONTINUE C IF (NMSG.LT.LENTAB) THEN C C Empty slot found for new message. C NMSG = NMSG + 1 LIBTAB(I) = LIB SUBTAB(I) = SUB MESTAB(I) = MES NERTAB(I) = NERR LEVTAB(I) = LEVEL KOUNT (I) = 1 ICOUNT = 1 ELSE C C Table is full. C KOUNTX = KOUNTX+1 ICOUNT = 0 ENDIF ENDIF RETURN C C Formats. C 9000 FORMAT ('0 ERROR MESSAGE SUMMARY' / + ' LIBRARY SUBROUTINE MESSAGE START NERR', + ' LEVEL COUNT') 9010 FORMAT (1X,A,3X,A,3X,A,3I10) 9020 FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10) 9030 FORMAT (1X) END *DECK XGETUA SUBROUTINE XGETUA (IUNITA, N) C***BEGIN PROLOGUE XGETUA C***PURPOSE Return unit number(s) to which error messages are being C sent. C***LIBRARY SLATEC (XERROR) C***CATEGORY R3C C***TYPE ALL (XGETUA-A) C***KEYWORDS ERROR, XERROR C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C Abstract C XGETUA may be called to determine the unit number or numbers C to which error messages are being sent. C These unit numbers may have been set by a call to XSETUN, C or a call to XSETUA, or may be a default value. C C Description of Parameters C --Output-- C IUNIT - an array of one to five unit numbers, depending C on the value of N. A value of zero refers to the C default unit, as defined by the I1MACH machine C constant routine. Only IUNIT(1),...,IUNIT(N) are C defined by XGETUA. The values of IUNIT(N+1),..., C IUNIT(5) are not defined (for N .LT. 5) or altered C in any way by XGETUA. C N - the number of units to which copies of the C error messages are being sent. N will be in the C range from 1 to 5. C C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC C Error-handling Package, SAND82-0800, Sandia C Laboratories, 1982. C***ROUTINES CALLED J4SAVE C***REVISION HISTORY (YYMMDD) C 790801 DATE WRITTEN C 861211 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE XGETUA DIMENSION IUNITA(5) C***FIRST EXECUTABLE STATEMENT XGETUA N = J4SAVE(5,0,.FALSE.) DO 30 I=1,N INDEX = I+4 IF (I.EQ.1) INDEX = 3 IUNITA(I) = J4SAVE(INDEX,0,.FALSE.) 30 CONTINUE RETURN END C------------------------------------------------------------- SUBROUTINE DQAG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,NEVAL,IER, * LIMIT,LENW,LAST,IWORK,WORK) C***BEGIN PROLOGUE DQAG C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A1 C***KEYWORDS AUTOMATIC INTEGRATOR, GENERAL-PURPOSE, C INTEGRAND EXAMINATOR, GLOBALLY ADAPTIVE, C GAUSS-KRONROD C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C***PURPOSE THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT)LE.MAX(EPSABS,EPSREL*ABS(I)). C***DESCRIPTION C C COMPUTATION OF A DEFINITE INTEGRAL C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C F - DOUBLE PRECISION C FUNCTION SUBPROGAM DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE DRIVER PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C EPSABS - DOUBLE PRECISION C ABSOLUTE ACCORACY REQUESTED C EPSREL - DOUBLE PRECISION C RELATIVE ACCURACY REQUESTED C IF EPSABS.LE.0 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), C THE ROUTINE WILL END WITH IER = 6. C C KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C A GAUSS-KRONROD PAIR IS USED WITH C 7 - 15 POINTS IF KEY.LT.2, C 10 - 21 POINTS IF KEY = 2, C 15 - 31 POINTS IF KEY = 3, C 20 - 41 POINTS IF KEY = 4, C 25 - 51 POINTS IF KEY = 5, C 30 - 61 POINTS IF KEY.GT.5. C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT) C C NEVAL - INTEGER C NUMBER OF INTEGRAND EVALUATIONS C C IER - INTEGER C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE C ROUTINE. IT IS ASSUMED THAT THE REQUESTED C ACCURACY HAS BEEN ACHIEVED. C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE C THE ESTIMATES FOR RESULT AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE C SUBDIVISIONS BY INCREASING THE VALUE OF C LIMIT (AND TAKING THE ACCORDING DIMENSION C ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF C THIS YIELD NO IMPROVEMENT IT IS ADVISED C TO ANALYZE THE INTEGRAND IN ORDER TO C DETERMINE THE INTEGRATION DIFFICULATIES. C IF THE POSITION OF A LOCAL DIFFICULTY CAN C BE DETERMINED (I.E.SINGULARITY, C DISCONTINUITY WITHIN THE INTERVAL) ONE C WILL PROBABLY GAIN FROM SPLITTING UP THE C INTERVAL AT THIS POINT AND CALLING THE C INTEGRATOR ON THE SUBRANGES. IF POSSIBLE, C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR C SHOULD BE USED WHICH IS DESIGNED FOR C HANDLING THE TYPE OF DIFFICULTY INVOLVED. C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS C DETECTED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS C AT SOME POINTS OF THE INTEGRATION C INTERVAL. C = 6 THE INPUT IS INVALID, BECAUSE C (EPSABS.LE.0 AND C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)) C OR LIMIT.LT.1 OR LENW.LT.LIMIT*4. C RESULT, ABSERR, NEVAL, LAST ARE SET C TO ZERO. C EXCEPT WHEN LENW IS INVALID, IWORK(1), C WORK(LIMIT*2+1) AND WORK(LIMIT*3+1) ARE C SET TO ZERO, WORK(1) IS SET TO A AND C WORK(LIMIT+1) TO B. C C DIMENSIONING PARAMETERS C LIMIT - INTEGER C DIMENSIONING PARAMETER FOR IWORK C LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS C IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL C (A,B), LIMIT.GE.1. C IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6. C C LENW - INTEGER C DIMENSIONING PARAMETER FOR WORK C LENW MUST BE AT LEAST LIMIT*4. C IF LENW.LT.LIMIT*4, THE ROUTINE WILL END WITH C IER = 6. C C LAST - INTEGER C ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS C PRODUCED IN THE SUBDIVIOSION PROCESS, WHICH C DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS C ACTUALLY IN THE WORK ARRAYS. C C WORK ARRAYS C IWORK - INTEGER C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K C ELEMENTS OF WHICH CONTAIN POINTERS TO THE ERROR C ESTIMATES OVER THE SUBINTERVALS, SUCH THAT C WORK(LIMIT*3+IWORK(1)),... , WORK(LIMIT*3+IWORK(K)) C FORM A DECREASING SEQUENCE WITH K = LAST IF C LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST OTHERWISE C C WORK - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LENW C ON RETURN C WORK(1), ..., WORK(LAST) CONTAIN THE LEFT END C POINTS OF THE SUBINTERVALS IN THE PARTITION OF C (A,B), C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN THE C RIGHT END POINTS, C WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS, C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) CONTAIN C THE ERROR ESTIMATES. C C***REFERENCES (NONE) C***ROUTINES CALLED DQAGE,XERROR C***END PROLOGUE DQAG DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK INTEGER IER,IWORK,KEY,LAST,LENW,LIMIT,LVL,L1,L2,L3,NEVAL C DIMENSION IWORK(LIMIT),WORK(LENW) C EXTERNAL F C C CHECK VALIDITY OF LENW. C C***FIRST EXECUTABLE STATEMENT DQAG IER = 6 NEVAL = 0 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10 C C PREPARE CALL FOR DQAGE. C L1 = LIMIT+1 L2 = LIMIT+L1 L3 = LIMIT+L2 C CALL DQAGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,NEVAL, * IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST) C C CALL ERROR HANDLER IF NECESSARY. C LVL = 0 10 IF(IER.EQ.6) LVL = 1 IF(IER.NE.0) CALL XERROR('ABNORMAL RETURN FROM DQAG' ,26,IER,LVL) RETURN END SUBROUTINE DQAGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR, * NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST) C***BEGIN PROLOGUE DQAGE C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A1 C***KEYWORDS AUTOMATIC INTEGRATOR, GENERAL-PURPOSE, C INTEGRAND EXAMINATOR, GLOBALLY ADAPTIVE, C GAUSS-KRONROD C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C***PURPOSE THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)). C***DESCRIPTION C C COMPUTATION OF A DEFINITE INTEGRAL C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS C ON ENTRY C F - DOUBLE PRECISION C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE DRIVER PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C EPSABS - DOUBLE PRECISION C ABSOLUTE ACCURACY REQUESTED C EPSREL - DOUBLE PRECISION C RELATIVE ACCURACY REQUESTED C IF EPSABS.LE.0 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), C THE ROUTINE WILL END WITH IER = 6. C C KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C A GAUSS-KRONROD PAIR IS USED WITH C 7 - 15 POINTS IF KEY.LT.2, C 10 - 21 POINTS IF KEY = 2, C 15 - 31 POINTS IF KEY = 3, C 20 - 41 POINTS IF KEY = 4, C 25 - 51 POINTS IF KEY = 5, C 30 - 61 POINTS IF KEY.GT.5. C C LIMIT - INTEGER C GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS C IN THE PARTITION OF (A,B), LIMIT.GE.1. C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT) C C NEVAL - INTEGER C NUMBER OF INTEGRAND EVALUATIONS C C IER - INTEGER C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE C ROUTINE. IT IS ASSUMED THAT THE REQUESTED C ACCURACY HAS BEEN ACHIEVED. C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE C THE ESTIMATES FOR RESULT AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE C SUBDIVISIONS BY INCREASING THE VALUE C OF LIMIT. C HOWEVER, IF THIS YIELDS NO IMPROVEMENT IT C IS RATHER ADVISED TO ANALYZE THE INTEGRAND C IN ORDER TO DETERMINE THE INTEGRATION C DIFFICULTIES. IF THE POSITION OF A LOCAL C DIFFICULTY CAN BE DETERMINED(E.G. C SINGULARITY, DISCONTINUITY WITHIN THE C INTERVAL) ONE WILL PROBABLY GAIN FROM C SPLITTING UP THE INTERVAL AT THIS POINT C AND CALLING THE INTEGRATOR ON THE C SUBRANGES. IF POSSIBLE, AN APPROPRIATE C SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED C WHICH IS DESIGNED FOR HANDLING THE TYPE OF C DIFFICULTY INVOLVED. C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS C DETECTED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS C AT SOME POINTS OF THE INTEGRATION C INTERVAL. C = 6 THE INPUT IS INVALID, BECAUSE C (EPSABS.LE.0 AND C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), C RESULT, ABSERR, NEVAL, LAST, RLIST(1) , C ELIST(1) AND IORD(1) ARE SET TO ZERO. C ALIST(1) AND BLIST(1) ARE SET TO A AND B C RESPECTIVELY. C C ALIST - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE LEFT C END POINTS OF THE SUBINTERVALS IN THE PARTITION C OF THE GIVEN INTEGRATION RANGE (A,B) C C BLIST - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE RIGHT C END POINTS OF THE SUBINTERVALS IN THE PARTITION C OF THE GIVEN INTEGRATION RANGE (A,B) C C RLIST - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE C INTEGRAL APPROXIMATIONS ON THE SUBINTERVALS C C ELIST - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE MODULI OF THE C ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS C C IORD - INTEGER C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K C ELEMENTS OF WHICH ARE POINTERS TO THE C ERROR ESTIMATES OVER THE SUBINTERVALS, C SUCH THAT ELIST(IORD(1)), ..., C ELIST(IORD(K)) FORM A DECREASING SEQUENCE, C WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND C K = LIMIT+1-LAST OTHERWISE C C LAST - INTEGER C NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE C SUBDIVISION PROCESS C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH,DQK15,DQK21,DQK31, C DQK41,DQK51,DQK61,DQPSRT C***END PROLOGUE DQAGE C DOUBLE PRECISION A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B, * BLIST,B1,B2,DABS,DEFABS,DEFAB1,DEFAB2,DMAX1,D1MACH,ELIST,EPMACH, * EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F, * RESABS,RESULT,RLIST,UFLOW INTEGER IER,IORD,IROFF1,IROFF2,K,KEY,KEYF,LAST,LIMIT,MAXERR,NEVAL, * NRMAX C DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), * RLIST(LIMIT) C EXTERNAL F C C LIST OF MAJOR VARIABLES C ----------------------- C C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS C CONSIDERED UP TO NOW C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS C CONSIDERED UP TO NOW C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER C (ALIST(I),BLIST(I)) C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I) C MAXERR - POINTER TO THE INTERVAL WITH LARGEST C ERROR ESTIMATE C ERRMAX - ELIST(MAXERR) C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL* C ABS(RESULT)) C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL C LAST - INDEX FOR SUBDIVISION C C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQAGE EPMACH = D1MACH(4) UFLOW = D1MACH(1) C C TEST ON VALIDITY OF PARAMETERS C ------------------------------ C IER = 0 NEVAL = 0 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 ALIST(1) = A BLIST(1) = B RLIST(1) = 0.0D+00 ELIST(1) = 0.0D+00 IORD(1) = 0 IF(EPSABS.LE.0.0D+00.AND. * EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) IER = 6 IF(IER.EQ.6) GO TO 999 C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- C KEYF = KEY IF(KEY.LE.0) KEYF = 1 IF(KEY.GE.7) KEYF = 6 NEVAL = 0 IF(KEYF.EQ.1) CALL DQK15(F,A,B,RESULT,ABSERR,DEFABS,RESABS) IF(KEYF.EQ.2) CALL DQK21(F,A,B,RESULT,ABSERR,DEFABS,RESABS) IF(KEYF.EQ.3) CALL DQK31(F,A,B,RESULT,ABSERR,DEFABS,RESABS) IF(KEYF.EQ.4) CALL DQK41(F,A,B,RESULT,ABSERR,DEFABS,RESABS) IF(KEYF.EQ.5) CALL DQK51(F,A,B,RESULT,ABSERR,DEFABS,RESABS) IF(KEYF.EQ.6) CALL DQK61(F,A,B,RESULT,ABSERR,DEFABS,RESABS) LAST = 1 RLIST(1) = RESULT ELIST(1) = ABSERR IORD(1) = 1 C C TEST ON ACCURACY. C ERRBND = DMAX1(EPSABS,EPSREL*DABS(RESULT)) IF(ABSERR.LE.0.5D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2 IF(LIMIT.EQ.1) IER = 1 IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS) * .OR.ABSERR.EQ.0.0D+00) GO TO 60 C C INITIALIZATION C -------------- C C ERRMAX = ABSERR MAXERR = 1 AREA = RESULT ERRSUM = ABSERR NRMAX = 1 IROFF1 = 0 IROFF2 = 0 C C MAIN DO-LOOP C ------------ C DO 30 LAST = 2,LIMIT C C BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE. C A1 = ALIST(MAXERR) B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR)) A2 = B1 B2 = BLIST(MAXERR) IF(KEYF.EQ.1) CALL DQK15(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) IF(KEYF.EQ.2) CALL DQK21(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) IF(KEYF.EQ.3) CALL DQK31(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) IF(KEYF.EQ.4) CALL DQK41(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) IF(KEYF.EQ.5) CALL DQK51(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) IF(KEYF.EQ.6) CALL DQK61(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) IF(KEYF.EQ.1) CALL DQK15(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2) IF(KEYF.EQ.2) CALL DQK21(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2) IF(KEYF.EQ.3) CALL DQK31(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2) IF(KEYF.EQ.4) CALL DQK41(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2) IF(KEYF.EQ.5) CALL DQK51(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2) IF(KEYF.EQ.6) CALL DQK61(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2) C C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL C AND ERROR AND TEST FOR ACCURACY. C NEVAL = NEVAL+1 AREA12 = AREA1+AREA2 ERRO12 = ERROR1+ERROR2 ERRSUM = ERRSUM+ERRO12-ERRMAX AREA = AREA+AREA12-RLIST(MAXERR) IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 5 IF(DABS(RLIST(MAXERR)-AREA12).LE.0.1D-04*DABS(AREA12) * .AND.ERRO12.GE.0.99D+00*ERRMAX) IROFF1 = IROFF1+1 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1 5 RLIST(MAXERR) = AREA1 RLIST(LAST) = AREA2 ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA)) IF(ERRSUM.LE.ERRBND) GO TO 8 C C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG. C IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2 C C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS C EQUALS LIMIT. C IF(LAST.EQ.LIMIT) IER = 1 C C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR C AT A POINT OF THE INTEGRATION RANGE. C IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03* * EPMACH)*(DABS(A2)+0.1D+04*UFLOW)) IER = 3 C C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST. C 8 IF(ERROR2.GT.ERROR1) GO TO 10 ALIST(LAST) = A2 BLIST(MAXERR) = B1 BLIST(LAST) = B2 ELIST(MAXERR) = ERROR1 ELIST(LAST) = ERROR2 GO TO 20 10 ALIST(MAXERR) = A2 ALIST(LAST) = A1 BLIST(LAST) = B1 RLIST(MAXERR) = AREA2 RLIST(LAST) = AREA1 ELIST(MAXERR) = ERROR2 ELIST(LAST) = ERROR1 C C CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL C WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT). C 20 CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) C ***JUMP OUT OF DO-LOOP IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40 30 CONTINUE C C COMPUTE FINAL RESULT. C --------------------- C 40 RESULT = 0.0D+00 DO 50 K=1,LAST RESULT = RESULT+RLIST(K) 50 CONTINUE ABSERR = ERRSUM 60 IF(KEYF.NE.1) NEVAL = (10*KEYF+1)*(2*NEVAL+1) IF(KEYF.EQ.1) NEVAL = 30*NEVAL+15 999 RETURN END SUBROUTINE DQK15(F,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE DQK15 C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A2 C***KEYWORDS 15-POINT GAUSS-KRONROD RULES C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV - K.U.LEUVEN C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR C ESTIMATE C J = INTEGRAL OF ABS(F) OVER (A,B) C***DESCRIPTION C C INTEGRATION RULES C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS C ON ENTRY C F - DOUBLE PRECISION C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE CALLING PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE 15-POINT C KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION C OF ABSCISSAE TO THE7-POINT GAUSS RULE(RESG). C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD NOT EXCEED ABS(I-RESULT) C C RESABS - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL J C C RESASC - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) C OVER (A,B) C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH C***END PROLOGUE DQK15 C DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1, * D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK INTEGER J,JTW,JTWM1 EXTERNAL F C DIMENSION FV1(7),FV2(7),WG(4),WGK(8),XGK(8) C C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR C CORRESPONDING WEIGHTS ARE GIVEN. C C XGK - ABSCISSAE OF THE 15-POINT KRONROD RULE C XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT C GAUSS RULE C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY C ADDED TO THE 7-POINT GAUSS RULE C C WGK - WEIGHTS OF THE 15-POINT KRONROD RULE C C WG - WEIGHTS OF THE 7-POINT GAUSS RULE C C C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, C BELL LABS, NOV. 1981. C DATA WG ( 1) / 0.1294849661 6886969327 0611432679 082 D0 / DATA WG ( 2) / 0.2797053914 8927666790 1467771423 780 D0 / DATA WG ( 3) / 0.3818300505 0511894495 0369775488 975 D0 / DATA WG ( 4) / 0.4179591836 7346938775 5102040816 327 D0 / C DATA XGK ( 1) / 0.9914553711 2081263920 6854697526 329 D0 / DATA XGK ( 2) / 0.9491079123 4275852452 6189684047 851 D0 / DATA XGK ( 3) / 0.8648644233 5976907278 9712788640 926 D0 / DATA XGK ( 4) / 0.7415311855 9939443986 3864773280 788 D0 / DATA XGK ( 5) / 0.5860872354 6769113029 4144838258 730 D0 / DATA XGK ( 6) / 0.4058451513 7739716690 6606412076 961 D0 / DATA XGK ( 7) / 0.2077849550 0789846760 0689403773 245 D0 / DATA XGK ( 8) / 0.0000000000 0000000000 0000000000 000 D0 / C DATA WGK ( 1) / 0.0229353220 1052922496 3732008058 970 D0 / DATA WGK ( 2) / 0.0630920926 2997855329 0700663189 204 D0 / DATA WGK ( 3) / 0.1047900103 2225018383 9876322541 518 D0 / DATA WGK ( 4) / 0.1406532597 1552591874 5189590510 238 D0 / DATA WGK ( 5) / 0.1690047266 3926790282 6583426598 550 D0 / DATA WGK ( 6) / 0.1903505780 6478540991 3256402421 014 D0 / DATA WGK ( 7) / 0.2044329400 7529889241 4161999234 649 D0 / DATA WGK ( 8) / 0.2094821410 8472782801 2999174891 714 D0 / C C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTERVAL C HLGTH - HALF-LENGTH OF THE INTERVAL C ABSC - ABSCISSA C FVAL* - FUNCTION VALUE C RESG - RESULT OF THE 7-POINT GAUSS FORMULA C RESK - RESULT OF THE 15-POINT KRONROD FORMULA C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B), C I.E. TO I/(B-A) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQK15 EPMACH = D1MACH(4) UFLOW = D1MACH(1) C CENTR = 0.5D+00*(A+B) HLGTH = 0.5D+00*(B-A) DHLGTH = DABS(HLGTH) C C COMPUTE THE 15-POINT KRONROD APPROXIMATION TO C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. C FC = F(CENTR) RESG = FC*WG(4) RESK = FC*WGK(8) RESABS = DABS(RESK) DO 10 J=1,3 JTW = J*2 ABSC = HLGTH*XGK(JTW) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTW) = FVAL1 FV2(JTW) = FVAL2 FSUM = FVAL1+FVAL2 RESG = RESG+WG(J)*FSUM RESK = RESK+WGK(JTW)*FSUM RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2)) 10 CONTINUE DO 15 J = 1,4 JTWM1 = J*2-1 ABSC = HLGTH*XGK(JTWM1) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTWM1) = FVAL1 FV2(JTWM1) = FVAL2 FSUM = FVAL1+FVAL2 RESK = RESK+WGK(JTWM1)*FSUM RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2)) 15 CONTINUE RESKH = RESK*0.5D+00 RESASC = WGK(8)*DABS(FC-RESKH) DO 20 J=1,7 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESABS = RESABS*DHLGTH RESASC = RESASC*DHLGTH ABSERR = DABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) * ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 * ((EPMACH*0.5D+02)*RESABS,ABSERR) RETURN END SUBROUTINE DQK21(F,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE DQK21 C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A2 C***KEYWORDS 21-POINT GAUSS-KRONROD RULES C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR C ESTIMATE C J = INTEGRAL OF ABS(F) OVER (A,B) C***DESCRIPTION C C INTEGRATION RULES C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS C ON ENTRY C F - DOUBLE PRECISION C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE DRIVER PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE 21-POINT C KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION C OF ABSCISSAE TO THE 10-POINT GAUSS RULE (RESG). C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD NOT EXCEED ABS(I-RESULT) C C RESABS - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL J C C RESASC - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) C OVER (A,B) C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH C***END PROLOGUE DQK21 C DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1, * D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK INTEGER J,JTW,JTWM1 EXTERNAL F C DIMENSION FV1(10),FV2(10),WG(5),WGK(11),XGK(11) C C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR C CORRESPONDING WEIGHTS ARE GIVEN. C C XGK - ABSCISSAE OF THE 21-POINT KRONROD RULE C XGK(2), XGK(4), ... ABSCISSAE OF THE 10-POINT C GAUSS RULE C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY C ADDED TO THE 10-POINT GAUSS RULE C C WGK - WEIGHTS OF THE 21-POINT KRONROD RULE C C WG - WEIGHTS OF THE 10-POINT GAUSS RULE C C C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, C BELL LABS, NOV. 1981. C DATA WG ( 1) / 0.0666713443 0868813759 3568809893 332 D0 / DATA WG ( 2) / 0.1494513491 5058059314 5776339657 697 D0 / DATA WG ( 3) / 0.2190863625 1598204399 5534934228 163 D0 / DATA WG ( 4) / 0.2692667193 0999635509 1226921569 469 D0 / DATA WG ( 5) / 0.2955242247 1475287017 3892994651 338 D0 / C DATA XGK ( 1) / 0.9956571630 2580808073 5527280689 003 D0 / DATA XGK ( 2) / 0.9739065285 1717172007 7964012084 452 D0 / DATA XGK ( 3) / 0.9301574913 5570822600 1207180059 508 D0 / DATA XGK ( 4) / 0.8650633666 8898451073 2096688423 493 D0 / DATA XGK ( 5) / 0.7808177265 8641689706 3717578345 042 D0 / DATA XGK ( 6) / 0.6794095682 9902440623 4327365114 874 D0 / DATA XGK ( 7) / 0.5627571346 6860468333 9000099272 694 D0 / DATA XGK ( 8) / 0.4333953941 2924719079 9265943165 784 D0 / DATA XGK ( 9) / 0.2943928627 0146019813 1126603103 866 D0 / DATA XGK ( 10) / 0.1488743389 8163121088 4826001129 720 D0 / DATA XGK ( 11) / 0.0000000000 0000000000 0000000000 000 D0 / C DATA WGK ( 1) / 0.0116946388 6737187427 8064396062 192 D0 / DATA WGK ( 2) / 0.0325581623 0796472747 8818972459 390 D0 / DATA WGK ( 3) / 0.0547558965 7435199603 1381300244 580 D0 / DATA WGK ( 4) / 0.0750396748 1091995276 7043140916 190 D0 / DATA WGK ( 5) / 0.0931254545 8369760553 5065465083 366 D0 / DATA WGK ( 6) / 0.1093871588 0229764189 9210590325 805 D0 / DATA WGK ( 7) / 0.1234919762 6206585107 7958109831 074 D0 / DATA WGK ( 8) / 0.1347092173 1147332592 8054001771 707 D0 / DATA WGK ( 9) / 0.1427759385 7706008079 7094273138 717 D0 / DATA WGK ( 10) / 0.1477391049 0133849137 4841515972 068 D0 / DATA WGK ( 11) / 0.1494455540 0291690566 4936468389 821 D0 / C C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTERVAL C HLGTH - HALF-LENGTH OF THE INTERVAL C ABSC - ABSCISSA C FVAL* - FUNCTION VALUE C RESG - RESULT OF THE 10-POINT GAUSS FORMULA C RESK - RESULT OF THE 21-POINT KRONROD FORMULA C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B), C I.E. TO I/(B-A) C C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQK21 EPMACH = D1MACH(4) UFLOW = D1MACH(1) C CENTR = 0.5D+00*(A+B) HLGTH = 0.5D+00*(B-A) DHLGTH = DABS(HLGTH) C C COMPUTE THE 21-POINT KRONROD APPROXIMATION TO C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. C RESG = 0.0D+00 FC = F(CENTR) RESK = WGK(11)*FC RESABS = DABS(RESK) DO 10 J=1,5 JTW = 2*J ABSC = HLGTH*XGK(JTW) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTW) = FVAL1 FV2(JTW) = FVAL2 FSUM = FVAL1+FVAL2 RESG = RESG+WG(J)*FSUM RESK = RESK+WGK(JTW)*FSUM RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2)) 10 CONTINUE DO 15 J = 1,5 JTWM1 = 2*J-1 ABSC = HLGTH*XGK(JTWM1) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTWM1) = FVAL1 FV2(JTWM1) = FVAL2 FSUM = FVAL1+FVAL2 RESK = RESK+WGK(JTWM1)*FSUM RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2)) 15 CONTINUE RESKH = RESK*0.5D+00 RESASC = WGK(11)*DABS(FC-RESKH) DO 20 J=1,10 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESABS = RESABS*DHLGTH RESASC = RESASC*DHLGTH ABSERR = DABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) * ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 * ((EPMACH*0.5D+02)*RESABS,ABSERR) RETURN END SUBROUTINE DQK31(F,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE DQK31 C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A2 C***KEYWORDS 31-POINT GAUSS-KRONROD RULES C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B) WITH ERROR C ESTIMATE C J = INTEGRAL OF ABS(F) OVER (A,B) C***DESCRIPTION C C INTEGRATION RULES C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS C ON ENTRY C F - DOUBLE PRECISION C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE CALLING PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE 31-POINT C GAUSS-KRONROD RULE (RESK), OBTAINED BY OPTIMAL C ADDITION OF ABSCISSAE TO THE 15-POINT GAUSS C RULE (RESG). C C ABSERR - DOUBLE PRECISON C ESTIMATE OF THE MODULUS OF THE MODULUS, C WHICH SHOULD NOT EXCEED ABS(I-RESULT) C C RESABS - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL J C C RESASC - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) C OVER (A,B) C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH C***END PROLOGUE DQK31 DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1, * D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK INTEGER J,JTW,JTWM1 EXTERNAL F C DIMENSION FV1(15),FV2(15),XGK(16),WGK(16),WG(8) C C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR C CORRESPONDING WEIGHTS ARE GIVEN. C C XGK - ABSCISSAE OF THE 31-POINT KRONROD RULE C XGK(2), XGK(4), ... ABSCISSAE OF THE 15-POINT C GAUSS RULE C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY C ADDED TO THE 15-POINT GAUSS RULE C C WGK - WEIGHTS OF THE 31-POINT KRONROD RULE C C WG - WEIGHTS OF THE 15-POINT GAUSS RULE C C C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, C BELL LABS, NOV. 1981. C DATA WG ( 1) / 0.0307532419 9611726835 4628393577 204 D0 / DATA WG ( 2) / 0.0703660474 8810812470 9267416450 667 D0 / DATA WG ( 3) / 0.1071592204 6717193501 1869546685 869 D0 / DATA WG ( 4) / 0.1395706779 2615431444 7804794511 028 D0 / DATA WG ( 5) / 0.1662692058 1699393355 3200860481 209 D0 / DATA WG ( 6) / 0.1861610000 1556221102 6800561866 423 D0 / DATA WG ( 7) / 0.1984314853 2711157645 6118326443 839 D0 / DATA WG ( 8) / 0.2025782419 2556127288 0620199967 519 D0 / C DATA XGK ( 1) / 0.9980022986 9339706028 5172840152 271 D0 / DATA XGK ( 2) / 0.9879925180 2048542848 9565718586 613 D0 / DATA XGK ( 3) / 0.9677390756 7913913425 7347978784 337 D0 / DATA XGK ( 4) / 0.9372733924 0070590430 7758947710 209 D0 / DATA XGK ( 5) / 0.8972645323 4408190088 2509656454 496 D0 / DATA XGK ( 6) / 0.8482065834 1042721620 0648320774 217 D0 / DATA XGK ( 7) / 0.7904185014 4246593296 7649294817 947 D0 / DATA XGK ( 8) / 0.7244177313 6017004741 6186054613 938 D0 / DATA XGK ( 9) / 0.6509967412 9741697053 3735895313 275 D0 / DATA XGK ( 10) / 0.5709721726 0853884753 7226737253 911 D0 / DATA XGK ( 11) / 0.4850818636 4023968069 3655740232 351 D0 / DATA XGK ( 12) / 0.3941513470 7756336989 7207370981 045 D0 / DATA XGK ( 13) / 0.2991800071 5316881216 6780024266 389 D0 / DATA XGK ( 14) / 0.2011940939 9743452230 0628303394 596 D0 / DATA XGK ( 15) / 0.1011420669 1871749902 7074231447 392 D0 / DATA XGK ( 16) / 0.0000000000 0000000000 0000000000 000 D0 / C DATA WGK ( 1) / 0.0053774798 7292334898 7792051430 128 D0 / DATA WGK ( 2) / 0.0150079473 2931612253 8374763075 807 D0 / DATA WGK ( 3) / 0.0254608473 2671532018 6874001019 653 D0 / DATA WGK ( 4) / 0.0353463607 9137584622 2037948478 360 D0 / DATA WGK ( 5) / 0.0445897513 2476487660 8227299373 280 D0 / DATA WGK ( 6) / 0.0534815246 9092808726 5343147239 430 D0 / DATA WGK ( 7) / 0.0620095678 0067064028 5139230960 803 D0 / DATA WGK ( 8) / 0.0698541213 1872825870 9520077099 147 D0 / DATA WGK ( 9) / 0.0768496807 5772037889 4432777482 659 D0 / DATA WGK ( 10) / 0.0830805028 2313302103 8289247286 104 D0 / DATA WGK ( 11) / 0.0885644430 5621177064 7275443693 774 D0 / DATA WGK ( 12) / 0.0931265981 7082532122 5486872747 346 D0 / DATA WGK ( 13) / 0.0966427269 8362367850 5179907627 589 D0 / DATA WGK ( 14) / 0.0991735987 2179195933 2393173484 603 D0 / DATA WGK ( 15) / 0.1007698455 2387559504 4946662617 570 D0 / DATA WGK ( 16) / 0.1013300070 1479154901 7374792767 493 D0 / C C C LIST OF MAJOR VARIABLES C ----------------------- C CENTR - MID POINT OF THE INTERVAL C HLGTH - HALF-LENGTH OF THE INTERVAL C ABSC - ABSCISSA C FVAL* - FUNCTION VALUE C RESG - RESULT OF THE 15-POINT GAUSS FORMULA C RESK - RESULT OF THE 31-POINT KRONROD FORMULA C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B), C I.E. TO I/(B-A) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C***FIRST EXECUTABLE STATEMENT DQK31 EPMACH = D1MACH(4) UFLOW = D1MACH(1) C CENTR = 0.5D+00*(A+B) HLGTH = 0.5D+00*(B-A) DHLGTH = DABS(HLGTH) C C COMPUTE THE 31-POINT KRONROD APPROXIMATION TO C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. C FC = F(CENTR) RESG = WG(8)*FC RESK = WGK(16)*FC RESABS = DABS(RESK) DO 10 J=1,7 JTW = J*2 ABSC = HLGTH*XGK(JTW) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTW) = FVAL1 FV2(JTW) = FVAL2 FSUM = FVAL1+FVAL2 RESG = RESG+WG(J)*FSUM RESK = RESK+WGK(JTW)*FSUM RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2)) 10 CONTINUE DO 15 J = 1,8 JTWM1 = J*2-1 ABSC = HLGTH*XGK(JTWM1) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTWM1) = FVAL1 FV2(JTWM1) = FVAL2 FSUM = FVAL1+FVAL2 RESK = RESK+WGK(JTWM1)*FSUM RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2)) 15 CONTINUE RESKH = RESK*0.5D+00 RESASC = WGK(16)*DABS(FC-RESKH) DO 20 J=1,15 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESABS = RESABS*DHLGTH RESASC = RESASC*DHLGTH ABSERR = DABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) * ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 * ((EPMACH*0.5D+02)*RESABS,ABSERR) RETURN END SUBROUTINE DQK41(F,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE DQK41 C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A2 C***KEYWORDS 41-POINT GAUSS-KRONROD RULES C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR C ESTIMATE C J = INTEGRAL OF ABS(F) OVER (A,B) C***DESCRIPTION C C INTEGRATION RULES C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS C ON ENTRY C F - DOUBLE PRECISION C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE CALLING PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE 41-POINT C GAUSS-KRONROD RULE (RESK) OBTAINED BY OPTIMAL C ADDITION OF ABSCISSAE TO THE 20-POINT GAUSS C RULE (RESG). C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD NOT EXCEED ABS(I-RESULT) C C RESABS - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL J C C RESASC - DOUBLE PRECISION C APPROXIMATION TO THE INTEGAL OF ABS(F-I/(B-A)) C OVER (A,B) C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH C***END PROLOGUE DQK41 C DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1, * D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK INTEGER J,JTW,JTWM1 EXTERNAL F C DIMENSION FV1(20),FV2(20),XGK(21),WGK(21),WG(10) C C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR C CORRESPONDING WEIGHTS ARE GIVEN. C C XGK - ABSCISSAE OF THE 41-POINT GAUSS-KRONROD RULE C XGK(2), XGK(4), ... ABSCISSAE OF THE 20-POINT C GAUSS RULE C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY C ADDED TO THE 20-POINT GAUSS RULE C C WGK - WEIGHTS OF THE 41-POINT GAUSS-KRONROD RULE C C WG - WEIGHTS OF THE 20-POINT GAUSS RULE C C C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, C BELL LABS, NOV. 1981. C DATA WG ( 1) / 0.0176140071 3915211831 1861962351 853 D0 / DATA WG ( 2) / 0.0406014298 0038694133 1039952274 932 D0 / DATA WG ( 3) / 0.0626720483 3410906356 9506535187 042 D0 / DATA WG ( 4) / 0.0832767415 7670474872 4758143222 046 D0 / DATA WG ( 5) / 0.1019301198 1724043503 6750135480 350 D0 / DATA WG ( 6) / 0.1181945319 6151841731 2377377711 382 D0 / DATA WG ( 7) / 0.1316886384 4917662689 8494499748 163 D0 / DATA WG ( 8) / 0.1420961093 1838205132 9298325067 165 D0 / DATA WG ( 9) / 0.1491729864 7260374678 7828737001 969 D0 / DATA WG ( 10) / 0.1527533871 3072585069 8084331955 098 D0 / C DATA XGK ( 1) / 0.9988590315 8827766383 8315576545 863 D0 / DATA XGK ( 2) / 0.9931285991 8509492478 6122388471 320 D0 / DATA XGK ( 3) / 0.9815078774 5025025919 3342994720 217 D0 / DATA XGK ( 4) / 0.9639719272 7791379126 7666131197 277 D0 / DATA XGK ( 5) / 0.9408226338 3175475351 9982722212 443 D0 / DATA XGK ( 6) / 0.9122344282 5132590586 7752441203 298 D0 / DATA XGK ( 7) / 0.8782768112 5228197607 7442995113 078 D0 / DATA XGK ( 8) / 0.8391169718 2221882339 4529061701 521 D0 / DATA XGK ( 9) / 0.7950414288 3755119835 0638833272 788 D0 / DATA XGK ( 10) / 0.7463319064 6015079261 4305070355 642 D0 / DATA XGK ( 11) / 0.6932376563 3475138480 5490711845 932 D0 / DATA XGK ( 12) / 0.6360536807 2651502545 2836696226 286 D0 / DATA XGK ( 13) / 0.5751404468 1971031534 2946036586 425 D0 / DATA XGK ( 14) / 0.5108670019 5082709800 4364050955 251 D0 / DATA XGK ( 15) / 0.4435931752 3872510319 9992213492 640 D0 / DATA XGK ( 16) / 0.3737060887 1541956067 2548177024 927 D0 / DATA XGK ( 17) / 0.3016278681 1491300432 0555356858 592 D0 / DATA XGK ( 18) / 0.2277858511 4164507808 0496195368 575 D0 / DATA XGK ( 19) / 0.1526054652 4092267550 5220241022 678 D0 / DATA XGK ( 20) / 0.0765265211 3349733375 4640409398 838 D0 / DATA XGK ( 21) / 0.0000000000 0000000000 0000000000 000 D0 / C DATA WGK ( 1) / 0.0030735837 1852053150 1218293246 031 D0 / DATA WGK ( 2) / 0.0086002698 5564294219 8661787950 102 D0 / DATA WGK ( 3) / 0.0146261692 5697125298 3787960308 868 D0 / DATA WGK ( 4) / 0.0203883734 6126652359 8010231432 755 D0 / DATA WGK ( 5) / 0.0258821336 0495115883 4505067096 153 D0 / DATA WGK ( 6) / 0.0312873067 7703279895 8543119323 801 D0 / DATA WGK ( 7) / 0.0366001697 5820079803 0557240707 211 D0 / DATA WGK ( 8) / 0.0416688733 2797368626 3788305936 895 D0 / DATA WGK ( 9) / 0.0464348218 6749767472 0231880926 108 D0 / DATA WGK ( 10) / 0.0509445739 2372869193 2707670050 345 D0 / DATA WGK ( 11) / 0.0551951053 4828599474 4832372419 777 D0 / DATA WGK ( 12) / 0.0591114008 8063957237 4967220648 594 D0 / DATA WGK ( 13) / 0.0626532375 5478116802 5870122174 255 D0 / DATA WGK ( 14) / 0.0658345971 3361842211 1563556969 398 D0 / DATA WGK ( 15) / 0.0686486729 2852161934 5623411885 368 D0 / DATA WGK ( 16) / 0.0710544235 5344406830 5790361723 210 D0 / DATA WGK ( 17) / 0.0730306903 3278666749 5189417658 913 D0 / DATA WGK ( 18) / 0.0745828754 0049918898 6581418362 488 D0 / DATA WGK ( 19) / 0.0757044976 8455667465 9542775376 617 D0 / DATA WGK ( 20) / 0.0763778676 7208073670 5502835038 061 D0 / DATA WGK ( 21) / 0.0766007119 1799965644 5049901530 102 D0 / C C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTERVAL C HLGTH - HALF-LENGTH OF THE INTERVAL C ABSC - ABSCISSA C FVAL* - FUNCTION VALUE C RESG - RESULT OF THE 20-POINT GAUSS FORMULA C RESK - RESULT OF THE 41-POINT KRONROD FORMULA C RESKH - APPROXIMATION TO MEAN VALUE OF F OVER (A,B), I.E. C TO I/(B-A) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQK41 EPMACH = D1MACH(4) UFLOW = D1MACH(1) C CENTR = 0.5D+00*(A+B) HLGTH = 0.5D+00*(B-A) DHLGTH = DABS(HLGTH) C C COMPUTE THE 41-POINT GAUSS-KRONROD APPROXIMATION TO C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. C RESG = 0.0D+00 FC = F(CENTR) RESK = WGK(21)*FC RESABS = DABS(RESK) DO 10 J=1,10 JTW = J*2 ABSC = HLGTH*XGK(JTW) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTW) = FVAL1 FV2(JTW) = FVAL2 FSUM = FVAL1+FVAL2 RESG = RESG+WG(J)*FSUM RESK = RESK+WGK(JTW)*FSUM RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2)) 10 CONTINUE DO 15 J = 1,10 JTWM1 = J*2-1 ABSC = HLGTH*XGK(JTWM1) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTWM1) = FVAL1 FV2(JTWM1) = FVAL2 FSUM = FVAL1+FVAL2 RESK = RESK+WGK(JTWM1)*FSUM RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2)) 15 CONTINUE RESKH = RESK*0.5D+00 RESASC = WGK(21)*DABS(FC-RESKH) DO 20 J=1,20 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESABS = RESABS*DHLGTH RESASC = RESASC*DHLGTH ABSERR = DABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D+00) * ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 * ((EPMACH*0.5D+02)*RESABS,ABSERR) RETURN END SUBROUTINE DQK51(F,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE DQK51 C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A2 C***KEYWORDS 51-POINT GAUSS-KRONROD RULES C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH & PROGR. DIV. - K.U.LEUVEN C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B) WITH ERROR C ESTIMATE C J = INTEGRAL OF ABS(F) OVER (A,B) C***DESCRIPTION C C INTEGRATION RULES C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS C ON ENTRY C F - DOUBLE PRECISION C FUNCTION SUBROUTINE DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE CALLING PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE 51-POINT C KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION C OF ABSCISSAE TO THE 25-POINT GAUSS RULE (RESG). C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD NOT EXCEED ABS(I-RESULT) C C RESABS - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL J C C RESASC - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) C OVER (A,B) C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH C***END PROLOGUE DQK51 C DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1, * D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK INTEGER J,JTW,JTWM1 EXTERNAL F C DIMENSION FV1(25),FV2(25),XGK(26),WGK(26),WG(13) C C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR C CORRESPONDING WEIGHTS ARE GIVEN. C C XGK - ABSCISSAE OF THE 51-POINT KRONROD RULE C XGK(2), XGK(4), ... ABSCISSAE OF THE 25-POINT C GAUSS RULE C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY C ADDED TO THE 25-POINT GAUSS RULE C C WGK - WEIGHTS OF THE 51-POINT KRONROD RULE C C WG - WEIGHTS OF THE 25-POINT GAUSS RULE C C C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, C BELL LABS, NOV. 1981. C DATA WG ( 1) / 0.0113937985 0102628794 7902964113 235 D0 / DATA WG ( 2) / 0.0263549866 1503213726 1901815295 299 D0 / DATA WG ( 3) / 0.0409391567 0130631265 5623487711 646 D0 / DATA WG ( 4) / 0.0549046959 7583519192 5936891540 473 D0 / DATA WG ( 5) / 0.0680383338 1235691720 7187185656 708 D0 / DATA WG ( 6) / 0.0801407003 3500101801 3234959669 111 D0 / DATA WG ( 7) / 0.0910282619 8296364981 1497220702 892 D0 / DATA WG ( 8) / 0.1005359490 6705064420 2206890392 686 D0 / DATA WG ( 9) / 0.1085196244 7426365311 6093957050 117 D0 / DATA WG ( 10) / 0.1148582591 4571164833 9325545869 556 D0 / DATA WG ( 11) / 0.1194557635 3578477222 8178126512 901 D0 / DATA WG ( 12) / 0.1222424429 9031004168 8959518945 852 D0 / DATA WG ( 13) / 0.1231760537 2671545120 3902873079 050 D0 / C DATA XGK ( 1) / 0.9992621049 9260983419 3457486540 341 D0 / DATA XGK ( 2) / 0.9955569697 9049809790 8784946893 902 D0 / DATA XGK ( 3) / 0.9880357945 3407724763 7331014577 406 D0 / DATA XGK ( 4) / 0.9766639214 5951751149 8315386479 594 D0 / DATA XGK ( 5) / 0.9616149864 2584251241 8130033660 167 D0 / DATA XGK ( 6) / 0.9429745712 2897433941 4011169658 471 D0 / DATA XGK ( 7) / 0.9207471152 8170156174 6346084546 331 D0 / DATA XGK ( 8) / 0.8949919978 7827536885 1042006782 805 D0 / DATA XGK ( 9) / 0.8658470652 9327559544 8996969588 340 D0 / DATA XGK ( 10) / 0.8334426287 6083400142 1021108693 570 D0 / DATA XGK ( 11) / 0.7978737979 9850005941 0410904994 307 D0 / DATA XGK ( 12) / 0.7592592630 3735763057 7282865204 361 D0 / DATA XGK ( 13) / 0.7177664068 1308438818 6654079773 298 D0 / DATA XGK ( 14) / 0.6735663684 7346836448 5120633247 622 D0 / DATA XGK ( 15) / 0.6268100990 1031741278 8122681624 518 D0 / DATA XGK ( 16) / 0.5776629302 4122296772 3689841612 654 D0 / DATA XGK ( 17) / 0.5263252843 3471918259 9623778158 010 D0 / DATA XGK ( 18) / 0.4730027314 4571496052 2182115009 192 D0 / DATA XGK ( 19) / 0.4178853821 9303774885 1814394594 572 D0 / DATA XGK ( 20) / 0.3611723058 0938783773 5821730127 641 D0 / DATA XGK ( 21) / 0.3030895389 3110783016 7478909980 339 D0 / DATA XGK ( 22) / 0.2438668837 2098843204 5190362797 452 D0 / DATA XGK ( 23) / 0.1837189394 2104889201 5969888759 528 D0 / DATA XGK ( 24) / 0.1228646926 1071039638 7359818808 037 D0 / DATA XGK ( 25) / 0.0615444830 0568507888 6546392366 797 D0 / DATA XGK ( 26) / 0.0000000000 0000000000 0000000000 000 D0 / C DATA WGK ( 1) / 0.0019873838 9233031592 6507851882 843 D0 / DATA WGK ( 2) / 0.0055619321 3535671375 8040236901 066 D0 / DATA WGK ( 3) / 0.0094739733 8617415160 7207710523 655 D0 / DATA WGK ( 4) / 0.0132362291 9557167481 3656405846 976 D0 / DATA WGK ( 5) / 0.0168478177 0912829823 1516667536 336 D0 / DATA WGK ( 6) / 0.0204353711 4588283545 6568292235 939 D0 / DATA WGK ( 7) / 0.0240099456 0695321622 0092489164 881 D0 / DATA WGK ( 8) / 0.0274753175 8785173780 2948455517 811 D0 / DATA WGK ( 9) / 0.0307923001 6738748889 1109020215 229 D0 / DATA WGK ( 10) / 0.0340021302 7432933783 6748795229 551 D0 / DATA WGK ( 11) / 0.0371162714 8341554356 0330625367 620 D0 / DATA WGK ( 12) / 0.0400838255 0403238207 4839284467 076 D0 / DATA WGK ( 13) / 0.0428728450 2017004947 6895792439 495 D0 / DATA WGK ( 14) / 0.0455029130 4992178890 9870584752 660 D0 / DATA WGK ( 15) / 0.0479825371 3883671390 6392255756 915 D0 / DATA WGK ( 16) / 0.0502776790 8071567196 3325259433 440 D0 / DATA WGK ( 17) / 0.0523628858 0640747586 4366712137 873 D0 / DATA WGK ( 18) / 0.0542511298 8854549014 4543370459 876 D0 / DATA WGK ( 19) / 0.0559508112 2041231730 8240686382 747 D0 / DATA WGK ( 20) / 0.0574371163 6156783285 3582693939 506 D0 / DATA WGK ( 21) / 0.0586896800 2239420796 1974175856 788 D0 / DATA WGK ( 22) / 0.0597203403 2417405997 9099291932 562 D0 / DATA WGK ( 23) / 0.0605394553 7604586294 5360267517 565 D0 / DATA WGK ( 24) / 0.0611285097 1705304830 5859030416 293 D0 / DATA WGK ( 25) / 0.0614711898 7142531666 1544131965 264 D0 / C NOTE: WGK (26) WAS CALCULATED FROM THE VALUES OF WGK(1..25) DATA WGK ( 26) / 0.0615808180 6783293507 8759824240 066 D0 / C C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTERVAL C HLGTH - HALF-LENGTH OF THE INTERVAL C ABSC - ABSCISSA C FVAL* - FUNCTION VALUE C RESG - RESULT OF THE 25-POINT GAUSS FORMULA C RESK - RESULT OF THE 51-POINT KRONROD FORMULA C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B), C I.E. TO I/(B-A) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQK51 EPMACH = D1MACH(4) UFLOW = D1MACH(1) C CENTR = 0.5D+00*(A+B) HLGTH = 0.5D+00*(B-A) DHLGTH = DABS(HLGTH) C C COMPUTE THE 51-POINT KRONROD APPROXIMATION TO C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. C FC = F(CENTR) RESG = WG(13)*FC RESK = WGK(26)*FC RESABS = DABS(RESK) DO 10 J=1,12 JTW = J*2 ABSC = HLGTH*XGK(JTW) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTW) = FVAL1 FV2(JTW) = FVAL2 FSUM = FVAL1+FVAL2 RESG = RESG+WG(J)*FSUM RESK = RESK+WGK(JTW)*FSUM RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2)) 10 CONTINUE DO 15 J = 1,13 JTWM1 = J*2-1 ABSC = HLGTH*XGK(JTWM1) FVAL1 = F(CENTR-ABSC) FVAL2 = F(CENTR+ABSC) FV1(JTWM1) = FVAL1 FV2(JTWM1) = FVAL2 FSUM = FVAL1+FVAL2 RESK = RESK+WGK(JTWM1)*FSUM RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2)) 15 CONTINUE RESKH = RESK*0.5D+00 RESASC = WGK(26)*DABS(FC-RESKH) DO 20 J=1,25 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESABS = RESABS*DHLGTH RESASC = RESASC*DHLGTH ABSERR = DABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) * ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 * ((EPMACH*0.5D+02)*RESABS,ABSERR) RETURN END SUBROUTINE DQK61(F,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE DQK61 C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A1A2 C***KEYWORDS 61-POINT GAUSS-KRONROD RULES C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B) WITH ERROR C ESTIMATE C J = INTEGRAL OF DABS(F) OVER (A,B) C***DESCRIPTION C C INTEGRATION RULE C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C C PARAMETERS C ON ENTRY C F - DOUBLE PRECISION C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE CALLING PROGRAM. C C A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE 61-POINT C KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION OF C ABSCISSAE TO THE 30-POINT GAUSS RULE (RESG). C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD EQUAL OR EXCEED DABS(I-RESULT) C C RESABS - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL J C C RESASC - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL OF DABS(F-I/(B-A)) C C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH C***END PROLOGUE DQK61 C DOUBLE PRECISION A,DABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1, * D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK INTEGER J,JTW,JTWM1 EXTERNAL F C DIMENSION FV1(30),FV2(30),XGK(31),WGK(31),WG(15) C C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE C INTERVAL (-1,1). BECAUSE OF SYMMETRY ONLY THE POSITIVE C ABSCISSAE AND THEIR CORRESPONDING WEIGHTS ARE GIVEN. C C XGK - ABSCISSAE OF THE 61-POINT KRONROD RULE C XGK(2), XGK(4) ... ABSCISSAE OF THE 30-POINT C GAUSS RULE C XGK(1), XGK(3) ... OPTIMALLY ADDED ABSCISSAE C TO THE 30-POINT GAUSS RULE C C WGK - WEIGHTS OF THE 61-POINT KRONROD RULE C C WG - WEIGTHS OF THE 30-POINT GAUSS RULE C C C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, C BELL LABS, NOV. 1981. C DATA WG ( 1) / 0.0079681924 9616660561 5465883474 674 D0 / DATA WG ( 2) / 0.0184664683 1109095914 2302131912 047 D0 / DATA WG ( 3) / 0.0287847078 8332336934 9719179611 292 D0 / DATA WG ( 4) / 0.0387991925 6962704959 6801936446 348 D0 / DATA WG ( 5) / 0.0484026728 3059405290 2938140422 808 D0 / DATA WG ( 6) / 0.0574931562 1761906648 1721689402 056 D0 / DATA WG ( 7) / 0.0659742298 8218049512 8128515115 962 D0 / DATA WG ( 8) / 0.0737559747 3770520626 8243850022 191 D0 / DATA WG ( 9) / 0.0807558952 2942021535 4694938460 530 D0 / DATA WG ( 10) / 0.0868997872 0108297980 2387530715 126 D0 / DATA WG ( 11) / 0.0921225222 3778612871 7632707087 619 D0 / DATA WG ( 12) / 0.0963687371 7464425963 9468626351 810 D0 / DATA WG ( 13) / 0.0995934205 8679526706 2780282103 569 D0 / DATA WG ( 14) / 0.1017623897 4840550459 6428952168 554 D0 / DATA WG ( 15) / 0.1028526528 9355884034 1285636705 415 D0 / C DATA XGK ( 1) / 0.9994844100 5049063757 1325895705 811 D0 / DATA XGK ( 2) / 0.9968934840 7464954027 1630050918 695 D0 / DATA XGK ( 3) / 0.9916309968 7040459485 8628366109 486 D0 / DATA XGK ( 4) / 0.9836681232 7974720997 0032581605 663 D0 / DATA XGK ( 5) / 0.9731163225 0112626837 4693868423 707 D0 / DATA XGK ( 6) / 0.9600218649 6830751221 6871025581 798 D0 / DATA XGK ( 7) / 0.9443744447 4855997941 5831324037 439 D0 / DATA XGK ( 8) / 0.9262000474 2927432587 9324277080 474 D0 / DATA XGK ( 9) / 0.9055733076 9990779854 6522558925 958 D0 / DATA XGK ( 10) / 0.8825605357 9205268154 3116462530 226 D0 / DATA XGK ( 11) / 0.8572052335 4606109895 8658510658 944 D0 / DATA XGK ( 12) / 0.8295657623 8276839744 2898119732 502 D0 / DATA XGK ( 13) / 0.7997278358 2183908301 3668942322 683 D0 / DATA XGK ( 14) / 0.7677774321 0482619491 7977340974 503 D0 / DATA XGK ( 15) / 0.7337900624 5322680472 6171131369 528 D0 / DATA XGK ( 16) / 0.6978504947 9331579693 2292388026 640 D0 / DATA XGK ( 17) / 0.6600610641 2662696137 0053668149 271 D0 / DATA XGK ( 18) / 0.6205261829 8924286114 0477556431 189 D0 / DATA XGK ( 19) / 0.5793452358 2636169175 6024932172 540 D0 / DATA XGK ( 20) / 0.5366241481 4201989926 4169793311 073 D0 / DATA XGK ( 21) / 0.4924804678 6177857499 3693061207 709 D0 / DATA XGK ( 22) / 0.4470337695 3808917678 0609900322 854 D0 / DATA XGK ( 23) / 0.4004012548 3039439253 5476211542 661 D0 / DATA XGK ( 24) / 0.3527047255 3087811347 1037207089 374 D0 / DATA XGK ( 25) / 0.3040732022 7362507737 2677107199 257 D0 / DATA XGK ( 26) / 0.2546369261 6788984643 9805129817 805 D0 / DATA XGK ( 27) / 0.2045251166 8230989143 8957671002 025 D0 / DATA XGK ( 28) / 0.1538699136 0858354696 3794672743 256 D0 / DATA XGK ( 29) / 0.1028069379 6673703014 7096751318 001 D0 / DATA XGK ( 30) / 0.0514718425 5531769583 3025213166 723 D0 / DATA XGK ( 31) / 0.0000000000 0000000000 0000000000 000 D0 / C DATA WGK ( 1) / 0.0013890136 9867700762 4551591226 760 D0 / DATA WGK ( 2) / 0.0038904611 2709988405 1267201844 516 D0 / DATA WGK ( 3) / 0.0066307039 1593129217 3319826369 750 D0 / DATA WGK ( 4) / 0.0092732796 5951776342 8441146892 024 D0 / DATA WGK ( 5) / 0.0118230152 5349634174 2232898853 251 D0 / DATA WGK ( 6) / 0.0143697295 0704580481 2451432443 580 D0 / DATA WGK ( 7) / 0.0169208891 8905327262 7572289420 322 D0 / DATA WGK ( 8) / 0.0194141411 9394238117 3408951050 128 D0 / DATA WGK ( 9) / 0.0218280358 2160919229 7167485738 339 D0 / DATA WGK ( 10) / 0.0241911620 7808060136 5686370725 232 D0 / DATA WGK ( 11) / 0.0265099548 8233310161 0601709335 075 D0 / DATA WGK ( 12) / 0.0287540487 6504129284 3978785354 334 D0 / DATA WGK ( 13) / 0.0309072575 6238776247 2884252943 092 D0 / DATA WGK ( 14) / 0.0329814470 5748372603 1814191016 854 D0 / DATA WGK ( 15) / 0.0349793380 2806002413 7499670731 468 D0 / DATA WGK ( 16) / 0.0368823646 5182122922 3911065617 136 D0 / DATA WGK ( 17) / 0.0386789456 2472759295 0348651532 281 D0 / DATA WGK ( 18) / 0.0403745389 5153595911 1995279752 468 D0 / DATA WGK ( 19) / 0.0419698102 1516424614 7147541285 970 D0 / DATA WGK ( 20) / 0.0434525397 0135606931 6831728117 073 D0 / DATA WGK ( 21) / 0.0448148001 3316266319 2355551616 723 D0 / DATA WGK ( 22) / 0.0460592382 7100698811 6271735559 374 D0 / DATA WGK ( 23) / 0.0471855465 6929915394 5261478181 099 D0 / DATA WGK ( 24) / 0.0481858617 5708712914 0779492298 305 D0 / DATA WGK ( 25) / 0.0490554345 5502977888 7528165367 238 D0 / DATA WGK ( 26) / 0.0497956834 2707420635 7811569379 942 D0 / DATA WGK ( 27) / 0.0504059214 0278234684 0893085653 585 D0 / DATA WGK ( 28) / 0.0508817958 9874960649 2297473049 805 D0 / DATA WGK ( 29) / 0.0512215478 4925877217 0656282604 944 D0 / DATA WGK ( 30) / 0.0514261285 3745902593 3862879215 781 D0 / DATA WGK ( 31) / 0.0514947294 2945156755 8340433647 099 D0 / C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTERVAL C HLGTH - HALF-LENGTH OF THE INTERVAL C DABSC - ABSCISSA C FVAL* - FUNCTION VALUE C RESG - RESULT OF THE 30-POINT GAUSS RULE C RESK - RESULT OF THE 61-POINT KRONROD RULE C RESKH - APPROXIMATION TO THE MEAN VALUE OF F C OVER (A,B), I.E. TO I/(B-A) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C EPMACH = D1MACH(4) UFLOW = D1MACH(1) C CENTR = 0.5D+00*(B+A) HLGTH = 0.5D+00*(B-A) DHLGTH = DABS(HLGTH) C C COMPUTE THE 61-POINT KRONROD APPROXIMATION TO THE C INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. C C***FIRST EXECUTABLE STATEMENT DQK61 RESG = 0.0D+00 FC = F(CENTR) RESK = WGK(31)*FC RESABS = DABS(RESK) DO 10 J=1,15 JTW = J*2 DABSC = HLGTH*XGK(JTW) FVAL1 = F(CENTR-DABSC) FVAL2 = F(CENTR+DABSC) FV1(JTW) = FVAL1 FV2(JTW) = FVAL2 FSUM = FVAL1+FVAL2 RESG = RESG+WG(J)*FSUM RESK = RESK+WGK(JTW)*FSUM RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2)) 10 CONTINUE DO 15 J=1,15 JTWM1 = J*2-1 DABSC = HLGTH*XGK(JTWM1) FVAL1 = F(CENTR-DABSC) FVAL2 = F(CENTR+DABSC) FV1(JTWM1) = FVAL1 FV2(JTWM1) = FVAL2 FSUM = FVAL1+FVAL2 RESK = RESK+WGK(JTWM1)*FSUM RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2)) 15 CONTINUE RESKH = RESK*0.5D+00 RESASC = WGK(31)*DABS(FC-RESKH) DO 20 J=1,30 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESABS = RESABS*DHLGTH RESASC = RESASC*DHLGTH ABSERR = DABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) * ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 * ((EPMACH*0.5D+02)*RESABS,ABSERR) RETURN END SUBROUTINE DQPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX) C***BEGIN PROLOGUE DQPSRT C***REFER TO DQAGE,DQAGIE,DQAGPE,DQAWSE C***ROUTINES CALLED (NONE) C***REVISION DATE 810101 (YYMMDD) C***KEYWORDS SEQUENTIAL SORTING C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN C***PURPOSE THIS ROUTINE MAINTAINS THE DESCENDING ORDERING IN THE C LIST OF THE LOCAL ERROR ESTIMATED RESULTING FROM THE C INTERVAL SUBDIVISION PROCESS. AT EACH CALL TWO ERROR C ESTIMATES ARE INSERTED USING THE SEQUENTIAL SEARCH C METHOD, TOP-DOWN FOR THE LARGEST ERROR ESTIMATE AND C BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE. C***DESCRIPTION C C ORDERING ROUTINE C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS (MEANING AT OUTPUT) C LIMIT - INTEGER C MAXIMUM NUMBER OF ERROR ESTIMATES THE LIST C CAN CONTAIN C C LAST - INTEGER C NUMBER OF ERROR ESTIMATES CURRENTLY IN THE LIST C C MAXERR - INTEGER C MAXERR POINTS TO THE NRMAX-TH LARGEST ERROR C ESTIMATE CURRENTLY IN THE LIST C C ERMAX - DOUBLE PRECISION C NRMAX-TH LARGEST ERROR ESTIMATE C ERMAX = ELIST(MAXERR) C C ELIST - DOUBLE PRECISION C VECTOR OF DIMENSION LAST CONTAINING C THE ERROR ESTIMATES C C IORD - INTEGER C VECTOR OF DIMENSION LAST, THE FIRST K ELEMENTS C OF WHICH CONTAIN POINTERS TO THE ERROR C ESTIMATES, SUCH THAT C ELIST(IORD(1)),..., ELIST(IORD(K)) C FORM A DECREASING SEQUENCE, WITH C K = LAST IF LAST.LE.(LIMIT/2+2), AND C K = LIMIT+1-LAST OTHERWISE C C NRMAX - INTEGER C MAXERR = IORD(NRMAX) C C***END PROLOGUE DQPSRT C DOUBLE PRECISION ELIST,ERMAX,ERRMAX,ERRMIN INTEGER I,IBEG,IDO,IORD,ISUCC,J,JBND,JUPBN,K,LAST,LIMIT,MAXERR, * NRMAX DIMENSION ELIST(LAST),IORD(LAST) C C CHECK WHETHER THE LIST CONTAINS MORE THAN C TWO ERROR ESTIMATES. C C***FIRST EXECUTABLE STATEMENT DQPSRT IF(LAST.GT.2) GO TO 10 IORD(1) = 1 IORD(2) = 2 GO TO 90 C C THIS PART OF THE ROUTINE IS ONLY EXECUTED IF, DUE TO A C DIFFICULT INTEGRAND, SUBDIVISION INCREASED THE ERROR C ESTIMATE. IN THE NORMAL CASE THE INSERT PROCEDURE SHOULD C START AFTER THE NRMAX-TH LARGEST ERROR ESTIMATE. C 10 ERRMAX = ELIST(MAXERR) IF(NRMAX.EQ.1) GO TO 30 IDO = NRMAX-1 DO 20 I = 1,IDO ISUCC = IORD(NRMAX-1) C ***JUMP OUT OF DO-LOOP IF(ERRMAX.LE.ELIST(ISUCC)) GO TO 30 IORD(NRMAX) = ISUCC NRMAX = NRMAX-1 20 CONTINUE C C COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO BE MAINTAINED C IN DESCENDING ORDER. THIS NUMBER DEPENDS ON THE NUMBER OF C SUBDIVISIONS STILL ALLOWED. C 30 JUPBN = LAST IF(LAST.GT.(LIMIT/2+2)) JUPBN = LIMIT+3-LAST ERRMIN = ELIST(LAST) C C INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN, C STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)). C JBND = JUPBN-1 IBEG = NRMAX+1 IF(IBEG.GT.JBND) GO TO 50 DO 40 I=IBEG,JBND ISUCC = IORD(I) C ***JUMP OUT OF DO-LOOP IF(ERRMAX.GE.ELIST(ISUCC)) GO TO 60 IORD(I-1) = ISUCC 40 CONTINUE 50 IORD(JBND) = MAXERR IORD(JUPBN) = LAST GO TO 90 C C INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP. C 60 IORD(I-1) = MAXERR K = JBND DO 70 J=I,JBND ISUCC = IORD(K) C ***JUMP OUT OF DO-LOOP IF(ERRMIN.LT.ELIST(ISUCC)) GO TO 80 IORD(K+1) = ISUCC K = K-1 70 CONTINUE IORD(I) = LAST GO TO 90 80 IORD(K+1) = LAST C C SET MAXERR AND ERMAX. C 90 MAXERR = IORD(NRMAX) ERMAX = ELIST(MAXERR) RETURN END C---------------------------------------------------------------- SUBROUTINE DSET(N,DA,DX,INCX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DX(*) DO 10 I=1,N,INCX DX(I)=DA 10 CONTINUE RETURN END C------------------------------------------------------------- SUBROUTINE TSLD(A,X,R,M) INTEGER M DOUBLE PRECISION A(1),X(M),R(1) C C TSLD CALLS TSLD1 TO SOLVE THE DOUBLE PRECISION LINEAR SYSTEM C A * X = B C WITH THE T - MATRIX A . C C ON ENTRY C C A DOUBLE PRECISION(2*M - 1) C THE FIRST ROW OF THE T - MATRIX FOLLOWED BY ITS C FIRST COLUMN BEGINNING WITH THE SECOND ELEMENT . C ON RETURN A IS UNALTERED . C C X DOUBLE PRECISION(M) C THE RIGHT HAND SIDE VECTOR B . C C R DOUBLE PRECISION(2*M - 2) C A WORK VECTOR . C C M INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C X THE SOLUTION VECTOR . C C TOEPLITZ PACKAGE. THIS VERSION DATED 07/23/82 . C C SUBROUTINES AND FUNCTIONS C C TOEPLITZ PACKAGE ... TSLD1 C C CALL SUBROUTINE TSLD1 C CALL TSLD1(A,A(M+1),X,X,R,R(M),M) C RETURN END C--------------------------------------------------------------------- SUBROUTINE TSLD1(A1,A2,B,X,C1,C2,M) INTEGER M DOUBLE PRECISION A1(M),A2(1),B(M),X(M),C1(1),C2(1) C C TSLD1 SOLVES THE DOUBLE PRECISION LINEAR SYSTEM C A * X = B C WITH THE T - MATRIX A . C C ON ENTRY C C A1 DOUBLE PRECISION(M) C THE FIRST ROW OF THE T - MATRIX A . C ON RETURN A1 IS UNALTERED . C C A2 DOUBLE PRECISION(M - 1) C THE FIRST COLUMN OF THE T - MATRIX A C BEGINNING WITH THE SECOND ELEMENT . C ON RETURN A2 IS UNALTERED . C C B DOUBLE PRECISION(M) C THE RIGHT HAND SIDE VECTOR . C ON RETURN B IS UNALTERED . C C C1 DOUBLE PRECISION(M - 1) C A WORK VECTOR . C C C2 DOUBLE PRECISION(M - 1) C A WORK VECTOR . C C M INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C X DOUBLE PRECISION(M) C THE SOLUTION VECTOR. X MAY COINCIDE WITH B . C C TOEPLITZ PACKAGE. THIS VERSION DATED 07/23/82 . C C INTERNAL VARIABLES C INTEGER I1,I2,N,N1,N2 DOUBLE PRECISION R1,R2,R3,R5,R6 C C SOLVE THE SYSTEM WITH THE PRINCIPAL MINOR OF ORDER 1 . C R1 = A1(1) X(1) = B(1)/R1 IF (M .EQ. 1) GO TO 80 C C RECURRENT PROCESS FOR SOLVING THE SYSTEM C WITH THE T - MATRIX FOR N = 2, M . C DO 70 N = 2, M C C COMPUTE MULTIPLES OF THE FIRST AND LAST COLUMNS OF C THE INVERSE OF THE PRINCIPAL MINOR OF ORDER N . C N1 = N - 1 N2 = N - 2 R5 = A2(N1) R6 = A1(N) IF (N .EQ. 2) GO TO 20 C1(N1) = R2 DO 10 I1 = 1, N2 I2 = N - I1 R5 = R5 + A2(I1)*C1(I2) R6 = R6 + A1(I1+1)*C2(I1) 10 CONTINUE 20 CONTINUE R2 = -R5/R1 R3 = -R6/R1 R1 = R1 + R5*R3 IF (N .EQ. 2) GO TO 40 R6 = C2(1) C2(N1) = 0.0D0 DO 30 I1 = 2, N1 R5 = C2(I1) C2(I1) = C1(I1)*R3 + R6 C1(I1) = C1(I1) + R6*R2 R6 = R5 30 CONTINUE 40 CONTINUE C2(1) = R3 C C COMPUTE THE SOLUTION OF THE SYSTEM WITH THE C PRINCIPAL MINOR OF ORDER N . C R5 = 0.0D0 DO 50 I1 = 1, N1 I2 = N - I1 R5 = R5 + A2(I1)*X(I2) 50 CONTINUE R6 = (B(N) - R5)/R1 DO 60 I1 = 1, N1 X(I1) = X(I1) + C2(I1)*R6 60 CONTINUE X(N) = R6 70 CONTINUE 80 CONTINUE RETURN END C---------------------------------------------------------------------- SUBROUTINE XERROR(STRING,INT,IER,LVL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*50 STRING WRITE(6,*) STRING, IER RETURN END C--------------------------------------------------------------------- SUBROUTINE ZBRENT(FUNC,X1,X2,ZERO,TOL,EPS,ITMAX,ITER1) IMPLICIT DOUBLE PRECISION(A-H,O-Z) EXTERNAL FUNC A=X1 B=X2 FA=FUNC(A) FB=FUNC(B) IF((FA.GT.0D0.AND.FB.GT.0D0).OR.(FA.LT.0D0.AND.FB.LT.0D0))PAUSE *'ROOT MUST BE BRACKETED FOR ZBRENT' C=B FC=FB DO 11 ITER=1,ITMAX IF((FB.GT.0D0.AND.FC.GT.0D0).OR.(FB.LT.0D0.AND.FC.LT.0D0))THEN C=A FC=FA D=B-A E=D ENDIF IF(DABS(FC).LT.DABS(FB)) THEN A=B B=C C=A FA=FB FB=FC FC=FA ENDIF TOL1=2D0*EPS*DABS(B)+0.5D0*TOL XM=.5D0*(C-B) IF(DABS(XM).LE.TOL1 .OR. FB.EQ.0D0)THEN ZERO=B ITER1=ITER RETURN ENDIF IF(DABS(E).GE.TOL1 .AND. DABS(FA).GT.DABS(FB)) THEN S=FB/FA IF(A.EQ.C) THEN P=2.*XM*S Q=1.-S ELSE Q=FA/FC R=FB/FC P=S*(2D0*XM*Q*(Q-R)-(B-A)*(R-1D0)) Q=(Q-1D0)*(R-1D0)*(S-1D0) ENDIF IF(P.GT.0D0) Q=-Q P=DABS(P) IF(2D0*P .LT. DMIN1(3D0*XM*Q-DABS(TOL1*Q),DABS(E*Q))) THEN E=D D=P/Q ELSE D=XM E=D ENDIF ELSE D=XM E=D ENDIF A=B FA=FB IF(DABS(D) .GT. TOL1) THEN B=B+D ELSE B=B+DSIGN(TOL1,XM) ENDIF FB=FUNC(B) 11 CONTINUE PAUSE 'ZBRENT EXCEEDING MAXIMUM ITERATIONS' ZERO=B ITER1=ITER RETURN END