SUBROUTINE EXPINT(X, N, KODE, M, TOL, EN, IERR) EXP 10 C C WRITTEN BY D.E. AMOS, SANDIA LABORATORIES, ALBUQUERQUE, NM, 87185 C C REFERENCE C COMPUTATION OF EXPONENTIAL INTEGRALS BY D.E. AMOS, ACM C TRANS. MATH SOFTWARE, 1980 C C ABSTRACT C EXPINT COMPUTES M MEMBER SEQUENCES OF EXPONENTIAL INTEGRALS C E(N+K,X), K=0,1,...,M-1 FOR N.GE.1 AND X.GE.0. THE POWER C SERIES IS IMPLEMENTED FOR X.LE.XCUT AND THE CONFLUENT C HYPERGEOMETRIC REPRESENTATION C C E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X) C C IS COMPUTED FOR X.GT.XCUT. SINCE SEQUENCES ARE COMPUTED IN A C STABLE FASHION BY RECURRING AWAY FROM X, A IS SELECTED AS THE C INTEGER CLOSEST TO X WITHIN THE CONSTRAINT N.LE.A.LE.N+M-1. C FOR THE U COMPUTATION A IS FURTHER MODIFIED TO BE THE C NEAREST EVEN INTEGER. INDICES ARE CARRIED FORWARD OR C BACKWARD BY THE TWO TERM RECURSION RELATION C C K*E(K+1,X) + X*E(K,X) = EXP(-X) C C ONCE E(A,X) IS COMPUTED. THE U FUNCTION IS COMPUTED BY MEANS C OF THE BACKWARD RECURSIVE MILLER ALGORITHM APPLIED TO THE C THREE TERM CONTIGUOUS RELATION FOR U(A+K,A,X), K=0,1,... C THIS PRODUCES ACCURATE RATIOS AND DETERMINES U(A+K,A,X),AND C HENCE E(A,X), TO WITHIN A MULTIPLICATIVE CONSTANT C. C ANOTHER CONTIGUOUS RELATION APPLIED TO C*U(A,A,X) AND C C*U(A+1,A,X) GETS C*U(A+1,A+1,X), A QUANTITY PROPORTIONAL TO C E(A+1,X). THE NORMALIZING CONSTANT C IS OBTAINED FROM THE C TWO TERM RECURSION RELATION ABOVE WITH K=A. C C MACHINE DEPENDENT PARAMETERS - XCUT, XLIM, ETOL, EULER, DIGAM C C EXPINT WRITES ERROR DIAGNOSTICS TO LOGICAL UNIT 3 C C DESCRIPTION OF ARGUMENTS C C INPUT C X X.GT.0.0 FOR N=1 AND X.GE.0.0 FOR N.GE.2 C N ORDER OF THE FIRST MEMBER OF THE SEQUENCE, N.GE.1 C KODE A SELECTION PARAMETER FOR SCALED VALUES C KODE=1 RETURNS E(N+K,X), K=0,1,...,M-1. C =2 RETURNS EXP(X)*E(N+K,X), K=0,1,...,M-1. C M NUMBER OF EXPONENTIAL INTEGRALS IN THE SEQUENCE, C M.GE.1 C TOL RELATIVE ACCURACY WANTED, ETOL.LE.TOL.LE.0.1 C ETOL=1.E-12 C C OUTPUT C EN A VECTOR OF DIMENSION AT LEAST M CONTAINING VALUES C EN(K) = E(N+K-1,X) OR EXP(X)*E(N+K-1,X), K=1,M C DEPENDING ON KODE C IERR UNDERFLOW INDICATOR C IERR=0 A NORMAL RETURN C =1 X EXCEEDS XLIM AND AN UNDERFLOW OCCURS. C EN(K)=0.0 , K=1,M RETURNED ON KODE=1 C XLIM=667. C C ERROR CONDITIONS C AN IMPROPER INPUT PARAMETER IS A FATAL ERROR C UNDERFLOW IS A NON FATAL ERROR. ZERO ANSWERS ARE RETURNED. C DIMENSION EN(1), A(99), B(99), Y(2) C DATA XCUT, XLIM, ETOL /2.0E0,667.0E0,1.0E-12/ DATA EULER /-5.77215664901533E-01/ DATA LUN /3/ C IF (N.LT.1) GO TO 260 IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 270 IF (M.LT.1) GO TO 280 IF (TOL.LT.ETOL .OR. TOL.GT.0.1E0) GO TO 290 C IERR = 0 IF (X.GT.XCUT) GO TO 100 IF (X.LT.0.0E0) GO TO 300 IF (X.EQ.0.0E0 .AND. N.EQ.1) GO TO 310 IF (X.EQ.0.0E0 .AND. N.GT.1) GO TO 80 C C SERIES FOR E(N,X) FOR X.LE.XCUT C IX = INT(X+0.5E0) C ICASE=1 MEANS INTEGER CLOSEST TO X IS 2 AND N=1 C ICASE=2 MEANS INTEGER CLOSEST TO X IS 0,1, OR 2 AND N.GE.2 ICASE = 2 IF (IX.GT.N) ICASE = 1 NM = N - ICASE + 1 ND = NM + 1 IND = 3 - ICASE MU = M - IND ML = 1 KS = ND FNM = FLOAT(NM) S = 0.0E0 XTOL = 3.0E0*TOL IF (ND.EQ.1) GO TO 10 XTOL = 0.3333E0*TOL S = 1.0E0/FNM 10 CONTINUE AA = 1.0E0 AK = 1.0E0 DO 50 I=1,35 AA = -AA*X/AK IF (I.EQ.NM) GO TO 30 S = S - AA/(AK-FNM) IF (ABS(AA).LE.XTOL*ABS(S)) GO TO 20 AK = AK + 1.0E0 GO TO 50 20 CONTINUE IF (I.LT.2) GO TO 40 IF (ND-2.GT.I .OR. I.GT.ND-1) GO TO 60 AK = AK + 1.0E0 GO TO 50 30 S = S + AA*(-ALOG(X)+DIGAM(ND)) XTOL = 3.0E0*TOL 40 AK = AK + 1.0E0 50 CONTINUE GO TO 320 60 IF (ND.EQ.1) S = S + (-ALOG(X)+EULER) IF (KODE.EQ.2) S = S*EXP(X) EN(1) = S EMX = 1.0E0 IF (M.EQ.1) GO TO 70 EN(IND) = S AA = FLOAT(KS) IF (KODE.EQ.1) EMX = EXP(-X) GO TO (220, 240), ICASE 70 IF (ICASE.EQ.2) RETURN IF (KODE.EQ.1) EMX = EXP(-X) EN(1) = (EMX-S)/X RETURN 80 CONTINUE DO 90 I=1,M EN(I) = 1.0E0/FLOAT(N+I-2) 90 CONTINUE RETURN C C BACKWARD RECURSIVE MILLER ALGORITHM FOR C E(N,X)=EXP(-X)*(X**(N-1))*U(N,N,X) C WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO X. C U(A,B,X) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION C 100 CONTINUE EMX = 1.0E0 IF (KODE.EQ.2) GO TO 130 IF (X.LE.XLIM) GO TO 120 IERR = 1 DO 110 I=1,M EN(I) = 0.0E0 110 CONTINUE RETURN 120 EMX = EXP(-X) 130 CONTINUE IX = INT(X+0.5E0) KN = N + M - 1 IF (KN.LE.IX) GO TO 140 IF (N.LT.IX .AND. IX.LT.KN) GO TO 170 IF (N.GE.IX) GO TO 160 GO TO 340 140 ICASE = 1 KS = KN ML = M - 1 MU = -1 IND = M IF (KN.GT.1) GO TO 180 150 KS = 2 ICASE = 3 GO TO 180 160 ICASE = 2 IND = 1 KS = N MU = M - 1 IF (N.GT.1) GO TO 180 IF (KN.EQ.1) GO TO 150 IX = 2 170 ICASE = 1 KS = IX ML = IX - N IND = ML + 1 MU = KN - IX 180 CONTINUE IK = KS/2 AH = FLOAT(IK) JSET = 1 + KS - (IK+IK) C START COMPUTATION FOR C EN(IND) = C*U( A , A ,X) JSET=1 C EN(IND) = C*U(A+1,A+1,X) JSET=2 C FOR AN EVEN INTEGER A. IC = 0 AA = AH + AH AAMS = AA - 1.0E0 AAMS = AAMS*AAMS TX = X + X FX = TX + TX AK = AH XTOL = TOL IF (TOL.LE.1.0E-3) XTOL = 20.0E0*TOL CT = AAMS + FX*AH EM = (AH+1.0E0)/((X+AA)*XTOL*SQRT(CT)) BK = AA CC = AH*AH C FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD C RECURSION P1 = 0.0E0 P2 = 1.0E0 190 CONTINUE IF (IC.EQ.99) GO TO 330 IC = IC + 1 AK = AK + 1.0E0 AT = BK/(BK+AK+CC+FLOAT(IC)) BK = BK + AK + AK A(IC) = AT BT = (AK+AK+X)/(AK+1.0E0) B(IC) = BT PT = P2 P2 = BT*P2 - AT*P1 P1 = PT CT = CT + FX EM = EM*AT*(1.0E0-TX/CT) IF (EM*(AK+1.0E0).GT.P1*P1) GO TO 190 ICT = IC KK = IC + 1 BT = TX/(CT+FX) Y2 = (BK/(BK+CC+FLOAT(KK)))*(P1/P2)*(1.0E0-BT+0.375E0*BT*BT) Y1 = 1.0E0 C BACKWARD RECURRENCE FOR C Y1= C*U( A ,A,X) C Y2= C*(A/(1+A/2))*U(A+1,A,X) DO 200 K=1,ICT KK = KK - 1 YT = Y1 Y1 = (B(KK)*Y1-Y2)/A(KK) Y2 = YT 200 CONTINUE C THE CONTIGUOUS RELATION C X*U(B,C+1,X)=(C-B)*U(B,C,X)+U(B-1,C,X) C WITH B=A+1 , C=A IS USED FOR C Y(2) = C * U(A+1,A+1,X) C X IS INCORPORATED INTO THE NORMALIZING RELATION FOR CNORM. PT=Y2/Y1 CNORM=1.0E0-PT*(AH+1.0E0)/AA Y(1)=1.0E0/(CNORM*AA+X) Y(2)=CNORM*Y(1) IF (ICASE.EQ.3) GO TO 210 EN(IND) = EMX*Y(JSET) IF (M.EQ.1) RETURN AA = FLOAT(KS) GO TO (220, 240), ICASE C C RECURSION SECTION N*E(N+1,X) + X*E(N,X)=EMX C 210 EN(1) = EMX*(1.0E0-Y(1))/X RETURN 220 K = IND - 1 DO 230 I=1,ML AA = AA - 1.0E0 EN(K) = (EMX-AA*EN(K+1))/X K = K - 1 230 CONTINUE IF (MU.LE.0) RETURN AA = FLOAT(KS) 240 K = IND DO 250 I=1,MU EN(K+1) = (EMX-X*EN(K))/AA AA = AA + 1.0E0 K = K + 1 250 CONTINUE RETURN C C 260 WRITE (LUN,99999) RETURN 270 WRITE (LUN,99998) RETURN 280 WRITE (LUN,99997) RETURN 290 WRITE (LUN,99996) RETURN 300 WRITE (LUN,99995) RETURN 310 WRITE (LUN,99994) RETURN 320 WRITE (LUN,99993) RETURN 330 WRITE (LUN,99992) RETURN 340 WRITE (LUN,99991) RETURN 99999 FORMAT (32H IN EXPINT, N NOT GREATER THAN 0) 99998 FORMAT (27H IN EXPINT, KODE NOT 1 OR 2) 99997 FORMAT (32H IN EXPINT, M NOT GREATER THAN 0) 99996 FORMAT (33H IN EXPINT, TOL NOT WITHIN LIMITS) 99995 FORMAT (37H IN EXPINT, X IS NOT ZERO OR POSITIVE) 99994 FORMAT (46H IN EXPINT, THE EXPONENTIAL INTEGRAL IS NOT DE, * 21HFINED FOR X=0 AND N=1) 99993 FORMAT (46H IN EXPINT, RELATIVE ERROR TEST FOR SERIES TER, * 28HMINATION NOT MET IN 36 TERMS) 99992 FORMAT (46H IN EXPINT, TERMINATION TEST FOR MILLER ALGORI, * 23HTHM NOT MET IN 99 STEPS) 99991 FORMAT (46H IN EXPINT, AN ERROR IN PLACING INT(X+0.5) WIT, * 47HH RESPECT TO N AND N+M-1 OCCURRED FOR X.GT.XCUT) END FUNCTION DIGAM(N) DIG 10 C C THIS SUBROUTINE RETURNS VALUES OF PSI(X)=DERIVATIVE OF LOG C GAMMA(X), X.GT.0.0 AT INTEGER ARGUMENTS. A TABLE LOOK-UP IS C PERFORMED FOR N.LE.100, AND THE ASYMPTOTIC EXPANSION IS C EVALUATED FOR N.GT.100. C DIMENSION B(4), C(100), C1(32), C2(27), C3(22), C4(19) EQUIVALENCE (C(1),C1(1)) EQUIVALENCE (C(33),C2(1)) EQUIVALENCE (C(60),C3(1)) EQUIVALENCE (C(82),C4(1)) C DATA C1 /-5.7721566490153E-01,4.22784335098467E-01, * 9.22784335098467E-01,1.25611766843180E+00,1.50611766843180E+00, * 1.70611766843180E+00,1.87278433509847E+00,2.01564147795561E+00, * 2.14064147795561E+00,2.25175258906672E+00,2.35175258906672E+00, * 2.44266167997581E+00,2.52599501330915E+00,2.60291809023222E+00, * 2.67434666166079E+00,2.74101332832746E+00,2.80351332832746E+00, * 2.86233685773923E+00,2.91789241329478E+00,2.97052399224215E+00, * 3.02052399224215E+00,3.06814303986120E+00,3.11359758531574E+00, * 3.15707584618531E+00,3.19874251285197E+00,3.23874251285197E+00, * 3.27720405131351E+00,3.31424108835055E+00,3.34995537406484E+00, * 3.38443813268552E+00,3.41777146601886E+00,3.45002953053499E+00/ DATA C2 /3.48127953053499E+00,3.51158256083802E+00, * 3.54099432554390E+00,3.56956575411533E+00,3.59734353189311E+00, * 3.62437055892013E+00,3.65068634839382E+00,3.67632737403484E+00, * 3.70132737403484E+00,3.72571761793728E+00,3.74952714174681E+00, * 3.77278295570029E+00,3.79551022842757E+00,3.81773245064979E+00, * 3.83947158108457E+00,3.86074817682925E+00,3.88158151016259E+00, * 3.90198967342789E+00,3.92198967342789E+00,3.94159751656515E+00, * 3.96082828579592E+00,3.97969621032422E+00,3.99821472884274E+00, * 4.01639654702455E+00,4.03425368988170E+00,4.05179754953082E+00, * 4.06903892884117E+00/ DATA C3 /4.08598808138354E+00,4.10265474805020E+00, * 4.11904819067316E+00,4.13517722293122E+00,4.15105023880424E+00, * 4.16667523880424E+00,4.18205985418885E+00,4.19721136934037E+00, * 4.21213674247470E+00,4.22684262482764E+00,4.24133537845082E+00, * 4.25562109273654E+00,4.26970559977879E+00,4.28359448866768E+00, * 4.29729311880467E+00,4.31080663231818E+00,4.32413996565151E+00, * 4.33729786038836E+00,4.35028487337537E+00,4.36310538619588E+00, * 4.37576361404398E+00,4.38826361404398E+00/ DATA C4 /4.40060929305633E+00,4.41280441500755E+00, * 4.42485260777863E+00,4.43675736968340E+00,4.44852207556575E+00, * 4.46014998254249E+00,4.47164423541606E+00,4.48300787177969E+00, * 4.49424382683587E+00,4.50535493794698E+00,4.51634394893599E+00, * 4.52721351415338E+00,4.53796620232543E+00,4.54860450019777E+00, * 4.55913081598724E+00,4.56954748265391E+00,4.57985676100442E+00, * 4.59006084263708E+00,4.60016185273809E+00/ C DATA B /1.66666666666667E-01,-3.33333333333333E-02, * 2.38095238095238E-02,-3.33333333333333E-02/ C IF (N.GT.100) GO TO 10 DIGAM = C(N) RETURN 10 FN = N AX = 1.0E0 AK = 2.0E0 S = -0.5E0/FN IF (FN.GT.1.E+8) GO TO 30 FN2 = FN*FN DO 20 K=1,3 AX = AX*FN2 S = S - B(K)/(AX*AK) AK = AK + 2.0E0 20 CONTINUE 30 CONTINUE DIGAM = S + ALOG(FN) RETURN END C PROGRAM TSTEXP(INPUT,OUTPUT,TAPE3=OUTPUT) 00000010 C 00000020 C PROGRAM TO TEST SUBROUTINE EXPINT AGAINST AN ADAPTIVE QUADRATURE. 00000030 C PARAMETER VALUES ARE PRINTED AND, IN THE EVENT THAT THE RELATIVE 00000040 C ERROR TEST IS NOT SATISFIED, X, ERROR, N, AND KODE ARE ALSO 00000050 C PRINTED. AN OUTPUT WITH ONLY PARAMETER VALUES INDICATES THAT ALL 00000060 C TESTS WERE PASSED. GAUS8 COMPUTES THE QUADRATURES. 00000070 C 00000080 DIMENSION XTOL(10), EN(50), EV(50) 00000090 IOUT = 3 00000100 NL = 1 00000110 NU = 16 00000120 NINC = 5 00000130 ML = 1 00000140 MU = 25 00000150 MINC = 8 00000160 KM = 5 00000170 JL = 1 00000180 JU = 40 00000190 JINC = 3 00000200 XTOL(1) = 1.0E-2 00000210 DO 10 I=2,3 00000220 XTOL(I) = XTOL(I-1)*1.0E-3 00000230 10 CONTINUE 00000240 DO 80 IT=1,3 00000250 TOL = XTOL(IT) 00000260 TOLA = AMAX1(1.0E-12,TOL/10.0E0) 00000270 BTOL = TOL 00000280 WRITE (IOUT,99999) TOL 00000290 DO 70 M=ML,MU,MINC 00000300 WRITE (IOUT,99998) M 00000310 DO 60 N=NL,NU,NINC 00000320 WRITE (IOUT,99997) N 00000330 DO 50 J=JL,JU,JINC 00000340 X = FLOAT(J-1)/5.0E0 00000350 EX = EXP(-X) 00000360 IF (X.EQ.0. .AND. N.EQ.1) GO TO 50 00000370 CALL EXPINT(X, N, 1, M, TOL, EN, IERR) 00000380 CALL EXPINT(X, N, 2, M, TOL, EV, IERR) 00000390 DO 40 K=1,M,KM 00000400 IF (X.GT.0.) GO TO 20 00000410 IF (N+K.EQ.2) GO TO 40 00000420 Y = 1.0E0/FLOAT(N+K-2) 00000430 YY = Y 00000440 GO TO 30 00000450 20 CONTINUE 00000460 NN = N + K - 1 00000470 YY = EINT(NN,X,TOLA,2) 00000480 Y = YY*EX 00000490 30 CONTINUE 00000500 ER = ABS((Y-EN(K))/Y) 00000510 KODE = 1 00000520 IF (ER.GT.BTOL) WRITE (IOUT,99996) X, ER, NN, KODE 00000530 KODE = 2 00000540 ERR = ABS((YY-EV(K))/YY) 00000550 IF (ERR.GT.BTOL) WRITE (IOUT,99996) X, ERR, NN, KODE 00000560 40 CONTINUE 00000570 50 CONTINUE 00000580 60 CONTINUE 00000590 70 CONTINUE 00000600 80 CONTINUE 00000610 STOP 00000620 99999 FORMAT (1H0, 5H TOL=, E15.4/) 00000630 99998 FORMAT (1H0, 2HM=, I5/) 00000640 99997 FORMAT (3X, 2HN=, I5) 00000650 99996 FORMAT (2E15.6, 2I5) 00000660 END 00000670 FUNCTION EINT(N, X, TOL, KODE) 00000680 COMMON /GEINT/ XX, FN EXTERNAL FEINT XX = X FN = N SIG = 1.0E0 S = 0.0E0 TOLA = TOL B = X 10 CONTINUE A = B REL = TOL B = B + SIG CALL GAUS8(FEINT, A, B, REL, ANS, IERR) S = S + ANS IF (ABS(ANS).LT.S*TOLA) GO TO 20 GO TO 10 20 EINT = S*EXP((FN-1.0E0)*ALOG(X)-FLOAT(2-KODE)*X) RETURN END FUNCTION FEINT(T) 00000880 COMMON /GEINT/ XX, FN FEINT = EXP(-T+XX-FN*ALOG(T)) RETURN END SUBROUTINE GAUS8 (FUN,A,B,ERR,ANS,IERR) 00000930 C C BY RONDALL E JONES, SANDIA LABORATORIES C SALIENT FEATURES -- INTERVAL BISECTION, COMBINED RELATIVE/ABSOLUTE C ERROR CONTROL, COMPUTED MAXIMUM REFINEMENT LEVEL WHEN A IS C CLOSE TO B. C C ABSTRACT C GAUS8 INTEGRATES REAL FUNCTIONS OF ONE VARIABLE OVER FINITE C INTERVALS, USING AN ADAPTIVE 8-POINT LEGENDRE-GAUSS ALGORITHM. C GAUS8 IS INTENDED PRIMARILY FOR HIGH ACCURACY INTEGRATION C OR INTEGRATION OF SMOOTH FUNCTIONS. FOR LOWER ACCURACY C INTEGRATION OF FUNCTIONS WHICH ARE NOT VERY SMOOTH, C EITHER QNC3 OR QNC7 MAY BE MORE EFFICIENT. C C DESCRIPTION OF ARGUMENTS C C INPUT-- C FUN - NAME OF EXTERNAL FUNCTION TO BE INTEGRATED. THIS NAME C MUST BE IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM. C FUN MUST BE A FUNCTION OF ONE REAL ARGUMENT. THE VALUE C OF THE ARGUMENT TO FUN IS THE VARIABLE OF INTEGRATION C WHICH RANGES FROM A TO B. C A - LOWER LIMIT OF INTEGRAL C B - UPPER LIMIT OF INTEGRAL (MAY BE LESS THAN A) C ERR - IS A REQUESTED ERROR TOLERANCE. NORMALLY PICK A VALUE OF C ABS(ERR).LT.1.E-3. ANS WILL NORMALLY HAVE NO MORE ERROR C THAN ABS(ERR) TIMES THE INTEGRAL OF THE ABSOLUTE VALUE C OF FUN(X). USUALLY, SMALLER VALUES FOR ERR YIELD C MORE ACCURACY AND REQUIRE MORE FUNCTION EVALUATIONS. C A NEGATIVE VALUE FOR ERR CAUSES AN ESTIMATE OF THE C ABSOLUTE ERROR IN ANS TO BE RETURNED IN ERR. C C OUTPUT-- C ERR - WILL BE AN ESTIMATE OF THE ERROR IN ANS IF THE INPUT C VALUE OF ERR WAS NEGATIVE. THE ESTIMATED ERROR IS SOLELY C FOR INFORMATION TO THE USER AND SHOULD NOT BE USED AS C A CORRECTION TO THE COMPUTED INTEGRAL. C ANS - COMPUTED VALUE OF INTEGRAL C IERR- A STATUS CODE C --NORMAL CODES C 1 ANS MOST LIKELY MEETS REQUESTED ERROR TOLERANCE, C OR A=B. C -1 A AND B ARE TOO NEARLY EQUAL TO ALLOW NORMAL C INTEGRATION. ANS IS SET TO ZERO. C --ABNORMAL CODE C 2 ANS PROBABLY DOES NOT MEET REQUESTED ERROR TOLERANCE. C C C C GAUS8 USES SUBROUTINES ERRCHK, ERRGET, ERRPRT, ERXSET, ERSTGT C COMPILE DECKS GAUS8, ERRCHK C DIMENSION AA(30),HH(30),LR(30),VL(30),GR(30) DATA X1,X2,X3,X4/0.18343 46424 95650 , 0.52553 24099 16329 , 1 0.79666 64774 13627 , 0.96028 98564 97536 / DATA W1,W2,W3,W4/0.36268 37833 78362 , 0.31370 66458 77887 , 1 0.22238 10344 53374 , 0.10122 85362 90376 / DATA SQ2/1.41421356/,ICALL/0/ DATA NLMN/1/,NLMX/30/,KMX/5000/,KML/6/,NBITS/48/ G8(X,H) = H*( (W1*(FUN(X-X1*H)+FUN(X+X1*H)) 1 +W2*(FUN(X-X2*H)+FUN(X+X2*H))) 2 +(W3*(FUN(X-X3*H)+FUN(X+X3*H)) 3 +W4*(FUN(X-X4*H)+FUN(X+X4*H))) ) C C INITIALIZE C IF(ICALL.NE.0)CALL ERRCHK(-71,71H*****GAUS8 CALLED RECURSIVELY. R 1ECURSIVE CALLS ARE ILLEGAL IN FORTRAN. ) ICALL = 1 ANS = 0.0 IERR = 1 CE = 0.0 IF (A.EQ.B) GO TO 35 LMX = NLMX LMN = NLMN IF (B.EQ.0.0) GO TO 4 IF (SIGN(1.0,B)*A.LE.0.0) GO TO 4 C = ABS(1.0-A/B) IF (C.GT.0.1) GO TO 4 IF (C.LE.0.0) GO TO 35 NIB = 0.5-ALOG(C)/ALOG(2.0) LMX = MIN0(NLMX , NBITS-NIB-7) IF (LMX.LT.1) GO TO 32 LMN = MIN0(LMN,LMX) 4 TOL = AMAX1(ABS(ERR),2.0**(5-NBITS))/2.0 IF (ERR.EQ.0.0) TOL = 0.5E-6 EPS = TOL HH(1) = (B-A)/4.0 AA(1) = A LR(1) = 1 L = 1 EST = G8(AA(L)+2.0*HH(L),2.0*HH(L)) K = 8 AREA = ABS(EST) EF = 0.5 MXL = 0 C C COMPUTE REFINED ESTIMATES, ESTIMATE THE ERROR, ETC. C 5 GL = G8(AA(L)+HH(L),HH(L)) GR(L) = G8(AA(L)+3.0*HH(L),HH(L)) K = K+16 AREA = AREA+(ABS(GL)+ABS(GR(L))-ABS(EST)) C IF (L.LT.LMN) GO TO 11 GLR = GL+GR(L) EE = ABS(EST-GLR)*EF AE = AMAX1(EPS*AREA,TOL*ABS(GLR)) IF (EE-AE) 8,8,10 7 MXL = 1 8 CE = CE + (EST-GLR) IF (LR(L)) 15,15,20 C C CONSIDER THE LEFT HALF OF THIS LEVEL C 10 IF (K.GT.KMX) LMX = KML IF (L.GE.LMX) GO TO 7 11 L = L+1 EPS = EPS*0.5 EF = EF/SQ2 HH(L) = HH(L-1)*0.5 LR(L) = -1 AA(L) = AA(L-1) EST = GL GO TO 5 C C PROCEED TO RIGHT HALF AT THIS LEVEL C 15 VL(L) = GLR 16 EST = GR(L-1) LR(L) = 1 AA(L) = AA(L)+4.0*HH(L) GO TO 5 C C RETURN ONE LEVEL C 20 VR = GLR 22 IF (L.LE.1) GO TO 30 L = L-1 EPS = EPS*2.0 EF = EF*SQ2 IF (LR(L)) 24,24,26 24 VL(L) = VL(L+1)+VR GO TO 16 26 VR = VL(L+1)+VR GO TO 22 C C EXIT C 30 ANS = VR IF ((MXL.EQ.0).OR.(ABS(CE).LE.2.0*TOL*AREA)) GO TO 35 IERR = 2 CALL ERRCHK(51,51HIN GAUS8 , ANS IS PROBABLY INSUFFICIENTLY ACCURA 1TE.) GO TO 35 32 IERR =-1 CALL ONECHK(-70,70HTHE FOLLOWING TEMPORARY INFORMATIVE DIAGNOSTIC + WILL APPEAR ONLY ONCE. ) CALLONECHK(-102,102HIN GAUS8 , A AND B ARE TOO NEARLY EQUAL TO ALL 1OW NORMAL INTEGRATION. ANS IS SET TO ZERO, AND IERR=-1.) 35 ICALL = 0 IF (ERR.LT.0.0) ERR = CE RETURN END SUBROUTINE ERRCHK(NCHARS,NARRAY) 00002570 C C SANDIA MATHEMATICAL PROGRAM LIBRARY C APPLIED MATHEMATICS DIVISION 2642 C SANDIA LABORATORIES C ALBUQUERQUE, NEW MEXICO 87115 C C SIMPLIFIED VERSION FOR STAND-ALONE USE. APRIL 1977 C C ABSTRACT C THE ROUTINES ERRCHK, ERXSET, AND ERRGET TOGETHER PROVIDE C A UNIFORM METHOD WITH SEVERAL OPTIONS FOR THE PROCESSING C OF DIAGNOSTICS AND WARNING MESSAGES WHICH ORIGINATE C IN THE MATHEMATICAL PROGRAM LIBRARY ROUTINES. C ERRCHK IS THE CENTRAL ROUTINE, WHICH ACTUALLY PROCESSES C MESSAGES. C C DESCRIPTION OF ARGUMENTS C NCHARS - NUMBER OF CHARACTERS IN HOLLERITH MESSAGE. C IF NCHARS IS NEGATED, ERRCHK WILL UNCONDITIONALLY C PRINT THE MESSAGE AND STOP EXECUTION. OTHERWISE, C THE BEHAVIOR OF ERRCHK MAY BE CONTROLLED BY C AN APPROPRIATE CALL TO ERXSET. C NARRAY - NAME OF ARRAY OR VARIABLE CONTAINING THE MESSAGE, C OR ELSE A LITERAL HOLLERITH CONSTANT CONTAINING C THE MESSAGE. BY CONVENTION, ALL MESSAGES SHOULD C BEGIN WITH *IN SUBNAM, ...*, WHERE SUBNAM IS THE C NAME OF THE ROUTINE CALLING ERRCHK. C C EXAMPLES C 1. TO ALLOW CONTROL BY CALLING ERXSET, USE C CALL ERRCHK(30,30HIN QUAD, INVALID VALUE OF ERR.) C 2. TO UNCONDITIONALLY PRINT A MESSAGE AND STOP EXECUTION, USE C CALL ERRCHK(-30,30HIN QUAD, INVALID VALUE OF ERR.) C C C C ERRCHK USES SUBROUTINES ERRGET, ERRPRT, ERXSET, ERSTGT C COMPILE DECKS ERRCHK C DIMENSION NARRAY(14) C IOUT=6 CALL ERRGET(NF,NT) C IF ERRCHK WAS CALLED WITH NEGATIVE CHARACTER COUNT, SET FATAL FLAG IF (NCHARS.LT.0) NF = -1 C IF MESSAGES ARE TO BE SUPPRESSED, RETURN IF (NF.EQ.0) RETURN C IF CHARACTER COUNT IS INVALID, STOP C IF (NCHARS.EQ.0) PRINT 5 IF (NCHARS .EQ. 0) WRITE (IOUT,5) 5 FORMAT(/31H ERRCHK WAS CALLED INCORRECTLY.) IF (NCHARS.EQ.0) STOP C PRINT MESSAGE CALL ERRPRT(IABS(NCHARS),NARRAY) C IF LAST MESSAGE, SAY SO C IF (NF.EQ.1) PRINT 10 IF (NF .EQ. 1) WRITE (IOUT,10) 10 FORMAT (30H ERRCHK MESSAGE LIMIT REACHED.) C PRINT TRACE-BACK IF ASKED TO C IF ((NT.GT.0).OR.(NF.LT.0)) CALL SYSTEM ROUTINE FOR TRACEBACK C DECREMENT MESSAGE COUNT IF (NF.GT.0) NF = NF-1 CALL ERXSET(NF,NT) C IF ALL IS WELL, RETURN IF (NF.GE.0) RETURN C IF THIS MESSAGE IS SUPPRESSABLE BY AN ERXSET CALL, C THEN EXPLAIN ERXSET USAGE. C IF (NCHARS.GT.0) PRINT 15 IF (NCHARS .GT. 0) WRITE (IOUT,15) 15 FORMAT (/13H *** NOTE *** 1/53H TO MAKE THE ERROR MESSAGE PRINTED ABOVE BE NONFATAL, 2/39H OR TO SUPPRESS THE MESSAGE COMPLETELY, 3/37H INSERT AN APPROPRIATE CALL TO ERXSET 4,30H AT THE START OF YOUR PROGRAM. 5/62H FOR EXAMPLE, TO PRINT UP TO 10 NONFATAL WARNING MESSAGES, USE 6/27H CALL ERXSET(10,0) ) C PRINT 20 WRITE (IOUT,20) 20 FORMAT (/28H PROGRAM ABORT DUE TO ERROR.) STOP END SUBROUTINE ONECHK(NCHARS,NARRAY) 00003390 C C ABSTRACT C ONECHK IS A COMPANION ROUTINE OF ERRCHK. IT IS CALLED C JUST LIKE ERRCHK, AND MESSAGES FROM IT MAY BE SUPPRESSED C BY AN APPROPRIATE CALL TO ERXSET. IT DIFFERS FROM ERRCHK C IN THAT EACH CALL TO ONECHK WILL PRODUCE NO MORE THAN ONE C PRINTED MESSAGE, REGARDLESS OF HOW MANY TIMES THAT CALL IS C EXECUTED, AND ONECHK NEVER TERMINATES EXECUTION. C ITS PURPOSE IS TO PROVIDE ONE-TIME-ONLY INFORMATIVE C DIAGNOSTICS. C C DESCRIPTION OF ARGUMENTS C NCHARS - NUMBER OF CHARACTERS IN THE MESSAGE. C IF NEGATED, THE MESSAGE WILL BE PRINTED (ONCE) EVEN C IF NFATAL HAS BEEN SET TO 0 (SEE ERXSET). C NARRAY - SAME AS IN ERRCHK C C C C ONECHK USES SUBROUTINES ERRGET, ERRPRT, ERXSET, ERSTGT C COMPILE DECKS ERRCHK C DIMENSION NARRAY(14) DATA NFLAG/4H.$,*/ IF (NARRAY(1).EQ.NFLAG) RETURN CALL ERRGET(NF,NT) IF ((NF.EQ.0).AND.(NCHARS.GT.0)) RETURN CALL ERRPRT (59,59HTHE FOLLOWING INFORMATIVE DIAGNOSTIC WILL APPEA 1R ONLY ONCE.) CALL ERRPRT(IABS(NCHARS),NARRAY) IF (NF.GT.0) NF = NF-1 CALL ERXSET(NF,NT) NARRAY(1) = NFLAG RETURN END SUBROUTINE ERRPRT(NCHARS,NARRAY) 00003750 C C UTILITY ROUTINE TO SIMPLY PRINT THE HOLLERITH MESSAGE IN NARRAY, C WHOSE LENGTH IS NCHARS CHARACTERS. C DIMENSION NARRAY(14) C C NOTE - NCH MUST BE THE NUMBER OF HOLLERITH CHARACTERS STORED C PER WORD. IF NCH IS CHANGED, FORMAT 1 MUST ALSO BE C CHANGED CORRESPONDINGLY. C IOUT=6 NCH = 10 C FOR LINE PRINTERS, USE 1 FORMAT (1X,13A10) C FOR DATA TERMINALS, USE C 1 FORMAT (1X,7A10) NWORDS = (NCHARS+NCH-1)/NCH C PRINT 1,(NARRAY(I),I=1,NWORDS) WRITE (IOUT,1) (NARRAY(I),I=1,NWORDS) RETURN END SUBROUTINE ERXSET(NFATAL,NTRACE) 00003970 C C ABSTRACT C ERXSET IS A COMPANION ROUTINE TO SUBROUTINE ERRCHK. C ERXSET ASSIGNS THE VALUES OF NFATAL AND NTRACE RESPECTIVELY C TO NF AND NT IN COMMON BLOCK MLBLK0 THEREBY SPECIFYING THE C STATE OF THE OPTIONS WHICH CONTROL THE EXECUTION OF ERRCHK. C C DESCRIPTION OF ARGUMENTS C BOTH ARGUMENTS ARE INPUT ARGUMENTS OF DATA TYPE INTEGER. C NFATAL - IS A FATAL-ERROR / MESSAGE-LIMIT FLAG. A NEGATIVE C VALUE DENOTES THAT DETECTED DIFFICULTIES ARE TO BE C TREATED AS FATAL ERRORS. NONNEGATIVE MEANS NONFATAL. C A NONNEGATIVE VALUE IS THE MAXIMUM NUMBER OF NONFATAL C WARNING MESSAGES WHICH WILL BE PRINTED BY ERRCHK, C AFTER WHICH NONFATAL MESSAGES WILL NOT BE PRINTED. C (DEFAULT VALUE IS -1.) C NTRACE - .GE.1 WILL CAUSE A TRACE-BACK TO BE GIVEN, C IF THIS FEATURE IS IMPLEMENTED ON THIS SYSTEM. C .LE.0 WILL SUPPRESS ANY TRACE-BACK, EXCEPT FOR C CASES WHEN EXECUTION IS TERMINATED. C (DEFAULT VALUE IS 0.) C C *NOTE* -- SOME CALLS TO ERRCHK WILL CAUSE UNCONDITIONAL C TERMINATION OF EXECUTION. ERXSET HAS NO EFFECT ON SUCH CALLS. C C EXAMPLES C 1. TO PRINT UP TO 100 MESSAGES AS NONFATAL WARNINGS USE C CALL ERXSET(100,0) C 2. TO SUPPRESS ALL MATHLIB WARNING MESSAGES USE C CALL ERXSET(0,0) C C C C ERXSET USES SUBROUTINES ERSTGT C COMPILE DECKS ERRCHK C CALL ERSTGT(0,NFATAL,NTRACE) RETURN END SUBROUTINE ERRGET(NFATAL,NTRACE) 00004370 C C ABSTRACT C ERRGET IS A COMPANION ROUTINE TO SUBROUTINE ERRCHK. C ERRGET ASSIGNS TO NFATAL AND NTRACE RESPECTIVELY THE VALUES C OF NF AND NT IN COMMON BLOCK MLBLK0 THEREBY ASCERTAINING THE C STATE OF THE OPTIONS WHICH CONTROL THE EXECUTION OF ERRCHK. C C DESCRIPTION OF ARGUMENTS C DESCRIPTION OF ARGUMENTS C BOTH ARGUMENTS ARE OUTPUT ARGUMENTS OF DATA TYPE INTEGER. C NFATAL - CURRENT VALUE OF NF (SEE DESCRIPTION OF ERXSET.) C NTRACE - CURRENT VALUE OF NT (SEE DESCRIPTION OF ERXSET.) C CALL ERSTGT(1,NFATAL,NTRACE) RETURN END SUBROUTINE ERSTGT(K,NFATAL,NTRACE) 00004530 C C THIS ROUTINE IS A SLAVE TO ERRGET AND ERRSET WHICH KEEPS C THE FLAGS AS LOCAL VARIABLES. C C *** IF LOCAL VARIABLES ARE NOT NORMALLY RETAINED BETWEEN C CALLS ON THIS SYSTEM, THE VARIABLES LNF AND LNT CAN BE C PLACED IN A COMMON BLOCK AND PRESET TO THE FOLLOWING C VALUES IN THE MAIN PROGRAM. C DATA LNF/-1/,LNT/0/ IF (K.LE.0) LNF = NFATAL IF (K.LE.0) LNT = NTRACE IF (K.GT.0) NFATAL = LNF IF (K.GT.0) NTRACE = LNT RETURN END .