C ALGORITHM 609, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 4, C DEC., 1983, P. 480-493. C PROGRAM QCKIN(INPUT,OUTPUT,TAPE6=OUTPUT) MAN 10 C MAN 20 C ABSTRACT MAN 30 C PROGRAM QCKIN IS A QUICK CHECK DRIVER WHICH EXERCISES THE MAJOR MAN 40 C LOOPS IN SUBROUTINE BSKIN (X,N,KODE,M,Y,IERR) FOR BICKLEY MAN 50 C FUNCTIONS KI(J,X). MORE PRECISELY, QCKIN DOES CONSISTENCY CHECKS MAN 60 C ON THE OUTPUT FROM BSKIN BY COMPARING SINGLE EVALUATIONS (M=1) MAN 70 C AGAINST SELECTED MEMBERS OF SEQUENCES WHICH ARE GENERATED BY MAN 80 C RECURSION. IF THE RELATIVE ERROR IS LESS THAN 1000 TIMES UNIT MAN 90 C ROUND OFF, THEN THE TEST IS PASSED - IF NOT, THEN X, THE VALUES MAN 100 C TO BE COMPARED, THE RELATIVE ERROR AND PARAMETERS KODE, N, M AND KMAN 110 C ARE WRITTEN ON LOGICAL UNIT 6 WHERE K IS THE MEMBER OF THE MAN 120 C SEQUENCE OF LENGTH M WHICH FAILED THE TEST. THAT IS, THE INDEX MAN 130 C OF THE FUNCTION WHICH FAILED THE TEST IS J=N+K-1. UNDERFLOW MAN 140 C TESTS ARE MADE AND ERROR CONDITIONS ARE TRIGGERED. MAN 150 C MAN 160 C FUNCTIONS I1MACH AND R1MACH MUST BE INITIALIZED ACCORDING TO THE MAN 170 C PROLOGUE IN EACH FUNCTION FOR THE MACHINE ENVIRONMENT BEFORE MAN 180 C QCKIN OR BSKIN CAN BE EXECUTED. FIFTEEN MACHINE ENVIRONMENTS MAN 190 C CAN BE DEFINED IN I1MACH AND R1MACH. MAN 200 C MAN 210 INTEGER I, IERR, IFLG, IX, I1M12, J, K, KODE, LUN, M, MDEL, MM, MAN 220 * N, NDEL, NN MAN 230 INTEGER I1MACH MAN 240 REAL AIX, ER, TOL, V, X, XINC, Y MAN 250 REAL R1MACH MAN 260 DIMENSION V(1), Y(10) MAN 270 DATA LUN /6/ MAN 280 C-----------------------------------------------------------------------MAN 290 C CALLS TO SLATEC PACKAGE TO OVERRIDE FATAL ERRORS AND WRITE MAN 300 C MESSSAGES TO UNIT LUN MAN 310 C CALL XSETF(-1) MAN 320 C CALL XSETUN(LUN) MAN 330 C-----------------------------------------------------------------------MAN 340 CALL XSETF(-1) MAN 350 CALL XSETUN(LUN) MAN 360 TOL = 1000.0E0*AMAX1(R1MACH(4),1.0E-18) MAN 370 IFLG = 0 MAN 380 WRITE (LUN,99999) MAN 390 99999 FORMAT (1H1//34H QUICK CHECK DIAGNOSTICS FOR BSKIN//) MAN 400 DO 70 KODE=1,2 MAN 410 N = 0 MAN 420 DO 60 NN=1,7 MAN 430 M = 1 MAN 440 DO 50 MM=1,4 MAN 450 X = 0.0E0 MAN 460 DO 40 IX=1,6 MAN 470 IF (N.EQ.0 .AND. IX.EQ.1) GO TO 30 MAN 480 CALL BSKIN(X, N, KODE, M, Y, IERR) MAN 490 DO 20 K=1,M,2 MAN 500 J = N + K - 1 MAN 510 CALL BSKIN(X, J, KODE, 1, V, IERR) MAN 520 ER = ABS((V(1)-Y(K))/V(1)) MAN 530 IF (ER.LE.TOL) GO TO 20 MAN 540 IF (IFLG.NE.0) GO TO 10 MAN 550 WRITE (LUN,99998) MAN 560 99998 FORMAT (8X, 1HX, 13X, 4HV(1), 11X, 4HY(K), 9X, 6HREL ER,MAN 570 * 1HR, 5X, 4HKODE, 3X, 1HN, 4X, 1HM, 4X, 1HK) MAN 580 10 CONTINUE MAN 590 IFLG = IFLG + 1 MAN 600 WRITE (LUN,99997) X, V(1), Y(K), ER, KODE, N, M, K MAN 610 99997 FORMAT (4E15.6, 4I5) MAN 620 IF (IFLG.GT.200) GO TO 130 MAN 630 20 CONTINUE MAN 640 30 CONTINUE MAN 650 AIX = FLOAT(2*IX-3) MAN 660 XINC = AMAX1(1.0E0,AIX) MAN 670 X = X + XINC MAN 680 40 CONTINUE MAN 690 MDEL = MAX0(1,MM-1) MAN 700 M = M + MDEL MAN 710 50 CONTINUE MAN 720 NDEL = MAX0(1,2*N-2) MAN 730 N = N + NDEL MAN 740 60 CONTINUE MAN 750 70 CONTINUE MAN 760 C-----------------------------------------------------------------------MAN 770 C TEST UNDERFLOW MAN 780 C-----------------------------------------------------------------------MAN 790 KODE = 1 MAN 800 M = 10 MAN 810 N = 10 MAN 820 I1M12 = I1MACH(12) MAN 830 X = -2.302E0*R1MACH(5)*FLOAT(I1M12) MAN 840 CALL BSKIN(X, N, KODE, M, Y, IERR) MAN 850 IF (IERR.EQ.1) GO TO 80 MAN 860 WRITE (LUN,99996) MAN 870 99996 FORMAT (//32H IERR IN UNDERFLOW TEST IS NOT 1//) MAN 880 IFLG = IFLG + 1 MAN 890 GO TO 110 MAN 900 80 CONTINUE MAN 910 DO 90 I=1,M MAN 920 IF (Y(I).NE.0.0E0) GO TO 100 MAN 930 90 CONTINUE MAN 940 GO TO 110 MAN 950 100 CONTINUE MAN 960 IFLG = IFLG + 1 MAN 970 WRITE (LUN,99995) MAN 980 99995 FORMAT (//43H SOME Y VALUE IN UNDERFLOW TEST IS NOT ZERO//) MAN 990 110 CONTINUE MAN 1000 IF (IFLG.NE.0) GO TO 120 MAN 1010 WRITE (LUN,99994) MAN 1020 99994 FORMAT (//16H QUICK CHECKS OK//) MAN 1030 120 CONTINUE MAN 1040 C-----------------------------------------------------------------------MAN 1050 C TRIGGER ERROR DIAGNOSTICS MAN 1060 C-----------------------------------------------------------------------MAN 1070 WRITE (LUN,99993) MAN 1080 99993 FORMAT (//27H TRIGGER 4 ERROR CONDITIONS//) MAN 1090 X = -1.0E0 MAN 1100 N = -1 MAN 1110 KODE = -1 MAN 1120 M = -1 MAN 1130 CALL BSKIN(X, N, KODE, M, Y, IERR) MAN 1140 STOP MAN 1150 130 CONTINUE MAN 1160 WRITE (LUN,99992) MAN 1170 99992 FORMAT (//52H PROCESSING OF MAIN LOOPS TERMINATED BECAUSE THE NUM,MAN 1180 * 36HBER OF DIAGNOSTIC PRINTS EXCEEDS 200//) MAN 1190 GO TO 120 MAN 1200 END MAN 1210 SUBROUTINE BSKIN(X, N, KODE, M, Y, IERR) BSK 10 C C WRITTEN BY D.E. AMOS, JUNE, 1982 C C REFERENCES C UNIFORM ASYMPTOTIC EXPANSIONS FOR EXPONENTIAL INTEGRALS E(N,X) C AND BICKLEY FUNCTIONS KI(N,X) BY D.E. AMOS, ACM TRANS. MATH C SOFTWARE, 1983 C C COMPUTATION OF BICKLEY FUNCTIONS BY D.E. AMOS, ACM TRANS. C MATH SOFTWARE, 1983 C C ABSTRACT C BSKIN COMPUTES M MEMBER SEQUENCES OF BICKLEY FUNCTIONS, C REPEATED INTEGRALS OF THE K0 BESSEL FUNCTION, Y(K)=KI(N+K-1,X) C ON KODE=1 OR SCALED FUNCTIONS Y(K)=EXP(X)*KI(N+K-1,X) ON C KODE=2, K=1,...,M, FOR N.GE.0 AND X.GE.0 (N AND X CANNOT BE C ZERO SIMULTANEOUSLY). C C NUMERICAL RECURRENCE ON C C (L-1)KI(L,X)=X(KI(L-3,X)-KI(L-1,X))+(L-2)KI(L-2,X) C C IS STABLE WHERE RECURRENCE IS CARRIED FORWARD OR BACKWARD C AWAY FROM INT(X+0.5). THE POWER SERIES FOR INDICES 0,1 AND 2 C ON 0.LE.X.LE.2 STARTS A STABLE RECURRENCE FOR INDICES C GREATER THAN 2. IF N IS SUFFICIENTLY LARGE (N.GT.NLIM), THE C UNIFORM ASYMPTOTIC EXPANSION FOR N TO INFINITY IS MORE C ECONOMICAL. ON X.GT.2 THE RECURSION IS STARTED BY EVALUATING C THE UNIFORM EXPANSION FOR THE THREE MEMBERS WHOSE INDICES ARE C CLOSEST TO INT(X+0.5) WITHIN THE SET N,...,N+M-1. FORWARD C RECURRENCE, BACKWARD RECURRENCE OR BOTH COMPLETE THE C SEQUENCE DEPENDING ON THE RELATION OF INT(X+0.5) TO THE C INDICES N,...,N+M-1. C C FUNCTIONS I1MACH AND R1MACH MUST BE INITIALIZED ACCORDING C TO PROLOGUE INSTRUCTIONS. FIFTEEN MACHINE ENVIRONMENTS C CAN BE DEFINED BY MAKING APPROPRIATE COMMENT CARDS ACTIVE C IN I1MACH AND R1MACH. C C THE NOMINAL COMPUTATIONAL ACCURACY IS THE MAXIMUM OF UNIT C ROUNDOFF (=R1MACH(4)) AND 1.0E-18 SINCE CRITICAL CONSTANTS C ARE GIVEN TO ONLY 18 DIGITS. C C BSKIN USES BKIAS,BKISR,EXINT,GAMRN,PSIXN,HKSEQ,BDIFF, C I1MACH,R1MACH,XERROR PACKAGE(18 ROUTINES) C C DBSKIN IS THE DOUBLE PRECISION VERSION OF BSKIN. C C DESCRIPTION OF ARGUMENTS C INPUT C X - ARGUMENT, X.GE.0.0E0 C N - ORDER OF FIRST MEMBER OF THE SEQUENCE N.GE.0 C KODE - SELECTION PARAMETER C KODE=1 RETURNS Y(K)= KI(N+K-1,X), K=1,M C =2 RETURNS Y(K)=EXP(X)*KI(N+K-1,X), K=1,M C M - NUMBER OF MEMBERS IN THE SEQUENCE, M.GE.1 C C OUTPUT C Y - A VECTOR OF DIMENSION AT LEAST M CONTAINING THE C SEQUENCE SELECTED BY KODE C IERR - UNDERFLOW FLAG C IERR=0 MEANS COMPUTATION COMPLETED C =1 MEANS AN EXPONENTIAL UNDERFLOW OCCURRED ON C KODE=1 AND Y(K)=0.0E0, K=1,...,M IS RETURNED C C ERROR CONDITIONS C AN IMPROPER INPUT IS A FATAL ERROR C***END PROLOGUE INTEGER I, ICASE, IERR, II, IL, ISET, I1M, K, KK, KODE, KTRMS, M, * M3, N, NE, NFLG, NL, NLIM, NN, NP, NS, NT INTEGER I1MACH REAL A, ENLIM, EXI, FN, GR, H, HN, HRTPI, SS, TOL, T1, T2, W, X, * XLIM, XNLIM, XP, Y, YS, YSS REAL GAMRN, R1MACH DIMENSION EXI(102), A(50), YS(3), YSS(3), H(31), Y(1) C----------------------------------------------------------------------- C COEFFICIENTS IN SERIES OF EXPONENTIAL INTEGRALS C----------------------------------------------------------------------- DATA A(1), A(2), A(3), A(4), A(5), A(6), A(7), A(8), A(9), A(10), * A(11), A(12), A(13), A(14), A(15), A(16), A(17), A(18), A(19), * A(20), A(21), A(22), A(23), A(24) /1.00000000000000000E+00, * 5.00000000000000000E-01,3.75000000000000000E-01, * 3.12500000000000000E-01,2.73437500000000000E-01, * 2.46093750000000000E-01,2.25585937500000000E-01, * 2.09472656250000000E-01,1.96380615234375000E-01, * 1.85470581054687500E-01,1.76197052001953125E-01, * 1.68188095092773438E-01,1.61180257797241211E-01, * 1.54981017112731934E-01,1.49445980787277222E-01, * 1.44464448094367981E-01,1.39949934091418982E-01, * 1.35833759559318423E-01,1.32060599571559578E-01, * 1.28585320635465905E-01,1.25370687619579257E-01, * 1.22385671247684513E-01,1.19604178719328047E-01, * 1.17004087877603524E-01/ DATA A(25), A(26), A(27), A(28), A(29), A(30), A(31), A(32), * A(33), A(34), A(35), A(36), A(37), A(38), A(39), A(40), A(41), * A(42), A(43), A(44), A(45), A(46), A(47), A(48) * /1.14566502713486784E-01,1.12275172659217048E-01, * 1.10116034723462874E-01,1.08076848895250599E-01, * 1.06146905164978267E-01,1.04316786110409676E-01, * 1.02578173008569515E-01,1.00923686347140974E-01, * 9.93467537479668965E-02,9.78414999033007314E-02, * 9.64026543164874854E-02,9.50254735405376642E-02, * 9.37056752969190855E-02,9.24393823875012600E-02, * 9.12230747245078224E-02,9.00535481254756708E-02, * 8.89278787739072249E-02,8.78433924473961612E-02, * 8.67976377754033498E-02,8.57883629175498224E-02, * 8.48134951571231199E-02,8.38711229887106408E-02, * 8.29594803475290034E-02,8.20769326842574183E-02/ DATA A(49), A(50) /8.12219646354630702E-02,8.03931690779583449E-02 * / C----------------------------------------------------------------------- C SQRT(PI)/2 C----------------------------------------------------------------------- DATA HRTPI /8.86226925452758014E-01/ C NFLG = 0 ISET = 1 IF (X.LT.0.0E0) GO TO 340 10 CONTINUE ISET = 2 IF (N.LT.0) GO TO 340 20 CONTINUE ISET = 3 IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 340 30 CONTINUE ISET = 4 IF (M.LT.1) GO TO 340 40 CONTINUE ISET = 5 IF (X.EQ.0.0E0 .AND. N.EQ.0) GO TO 340 50 CONTINUE IF (NFLG.NE.0) GO TO 410 IERR = 0 IF (X.EQ.0.0E0) GO TO 300 I1M = -I1MACH(12) T1 = 2.3026E0*R1MACH(5)*FLOAT(I1M) XLIM = T1 - 3.228086E0 T2 = T1 + FLOAT(N+M-1) IF (T2.GT.1000.0E0) XLIM = T1 - 0.5E0*(ALOG(T2)-0.451583E0) IF (X.GT.XLIM .AND. KODE.EQ.1) GO TO 320 TOL = AMAX1(R1MACH(4),1.0E-18) I1M = I1MACH(11) C----------------------------------------------------------------------- C LN(NLIM) = 0.125*LN(EPS), NLIM = 2*KTRMS+N C----------------------------------------------------------------------- XNLIM = 0.287823E0*FLOAT(I1M-1)*R1MACH(5) ENLIM = EXP(XNLIM) NLIM = INT(ENLIM) + 2 NLIM = MIN0(100,NLIM) NLIM = MAX0(20,NLIM) M3 = MIN0(M,3) NL = N + M - 1 IF (X.GT.2.0E0) GO TO 130 IF (N.GT.NLIM) GO TO 280 C----------------------------------------------------------------------- C COMPUTATION BY SERIES FOR 0.LE.X.LE.2 C----------------------------------------------------------------------- NFLG = 0 NN = N IF (NL.LE.2) GO TO 60 M3 = 3 NN = 0 NFLG = 1 60 CONTINUE XP = 1.0E0 IF (KODE.EQ.2) XP = EXP(X) DO 80 I=1,M3 CALL BKISR(X, NN, W) W = W*XP IF (NN.LT.N) GO TO 70 KK = NN - N + 1 Y(KK) = W 70 CONTINUE YS(I) = W NN = NN + 1 80 CONTINUE IF (NFLG.EQ.0) RETURN NS = NN XP = 1.0E0 90 CONTINUE C----------------------------------------------------------------------- C FORWARD RECURSION SCALED BY EXP(X) ON ICASE=0,1,2 C----------------------------------------------------------------------- FN = FLOAT(NS-1) IL = NL - NS + 1 IF (IL.LE.0) GO TO 120 DO 110 I=1,IL T1 = YS(2) T2 = YS(3) YS(3) = (X*(YS(1)-YS(3))+(FN-1.0E0)*YS(2))/FN YS(2) = T2 YS(1) = T1 FN = FN + 1.0E0 IF (NS.LT.N) GO TO 100 KK = NS - N + 1 Y(KK) = YS(3)*XP 100 CONTINUE NS = NS + 1 110 CONTINUE 120 CONTINUE RETURN C----------------------------------------------------------------------- C COMPUTATION BY ASYMPTOTIC EXPANSION FOR X.GT.2 C----------------------------------------------------------------------- 130 CONTINUE W = X + 0.5E0 NT = INT(W) IF (NL.GT.NT) GO TO 270 C----------------------------------------------------------------------- C CASE NL.LE.NT, ICASE=0 C----------------------------------------------------------------------- ICASE = 0 NN = NL NFLG = MIN0(M-M3,1) 140 CONTINUE KK = (NLIM-NN)/2 KTRMS = MAX0(0,KK) NS = NN + 1 NP = NN - M3 + 1 XP = 1.0E0 IF (KODE.EQ.1) XP = EXP(-X) DO 150 I=1,M3 KK = I CALL BKIAS(X, NP, KTRMS, A, W, KK, NE, GR, H) YS(I) = W NP = NP + 1 150 CONTINUE C----------------------------------------------------------------------- C SUM SERIES OF EXPONENTIAL INTEGRALS BACKWARD C----------------------------------------------------------------------- IF (KTRMS.EQ.0) GO TO 160 NE = KTRMS + KTRMS + 1 NP = NN - M3 + 2 CALL EXINT(X, NP, 2, NE, TOL, EXI, IERR) IF (IERR.EQ.1) GO TO 320 160 CONTINUE DO 190 I=1,M3 SS = 0.0E0 IF (KTRMS.EQ.0) GO TO 180 KK = I + KTRMS + KTRMS - 2 IL = KTRMS DO 170 K=1,KTRMS SS = SS + A(IL)*EXI(KK) KK = KK - 2 IL = IL - 1 170 CONTINUE 180 CONTINUE YS(I) = YS(I) + SS 190 CONTINUE IF (ICASE.EQ.1) GO TO 200 IF (NFLG.NE.0) GO TO 220 200 CONTINUE DO 210 I=1,M3 Y(I) = YS(I)*XP 210 CONTINUE IF (ICASE.EQ.1 .AND. NFLG.EQ.1) GO TO 90 RETURN 220 CONTINUE C----------------------------------------------------------------------- C BACKWARD RECURSION SCALED BY EXP(X) ICASE=0,2 C----------------------------------------------------------------------- KK = NN - N + 1 K = M3 DO 230 I=1,M3 Y(KK) = YS(K)*XP YSS(I) = YS(I) KK = KK - 1 K = K - 1 230 CONTINUE IL = KK IF (IL.LE.0) GO TO 250 FN = FLOAT(NN-3) DO 240 I=1,IL T1 = YS(2) T2 = YS(1) YS(1) = YS(2) + ((FN+2.0E0)*YS(3)-(FN+1.0E0)*YS(1))/X YS(2) = T2 YS(3) = T1 Y(KK) = YS(1)*XP KK = KK - 1 FN = FN - 1.0E0 240 CONTINUE 250 CONTINUE IF (ICASE.NE.2) RETURN DO 260 I=1,M3 YS(I) = YSS(I) 260 CONTINUE GO TO 90 270 CONTINUE IF (N.LT.NT) GO TO 290 C----------------------------------------------------------------------- C ICASE=1, NT.LE.N.LE.NL WITH FORWARD RECURSION C----------------------------------------------------------------------- 280 CONTINUE NN = N + M3 - 1 NFLG = MIN0(M-M3,1) ICASE = 1 GO TO 140 C----------------------------------------------------------------------- C ICASE=2, N.LT.NT.LT.NL WITH BOTH FORWARD AND BACKWARD RECURSION C----------------------------------------------------------------------- 290 CONTINUE NN = NT + 1 NFLG = MIN0(M-M3,1) ICASE = 2 GO TO 140 C----------------------------------------------------------------------- C X=0 CASE C----------------------------------------------------------------------- 300 CONTINUE FN = FLOAT(N) HN = 0.5E0*FN GR = GAMRN(HN) Y(1) = HRTPI*GR IF (M.EQ.1) RETURN Y(2) = HRTPI/(HN*GR) IF (M.EQ.2) RETURN DO 310 K=3,M Y(K) = FN*Y(K-2)/(FN+1.0E0) FN = FN + 1.0E0 310 CONTINUE RETURN C----------------------------------------------------------------------- C UNDERFLOW ON KODE=1, X.GT.XLIM C----------------------------------------------------------------------- 320 CONTINUE IERR = 1 DO 330 I=1,M Y(I) = 0.0E0 330 CONTINUE RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 340 IF (NFLG.NE.0) GO TO 350 CALL XGETF(II) CALL XSETF(-1) 350 GO TO (360, 370, 380, 390, 400), ISET 360 CALL XERROR(33H BSKIN, X IS NOT ZERO OR POSITIVE, 33, 2, 1) NFLG = 1 GO TO 10 370 CALL XERROR(33H BSKIN, N IS NOT ZERO OR POSITIVE, 33, 2, 1) NFLG = 1 GO TO 20 380 CALL XERROR(26H BSKIN, KODE IS NOT 1 OR 2, 26, 2, 1) NFLG = 1 GO TO 30 390 CALL XERROR(34H BSKIN, M IS NOT GREATER THAN ZERO, 34, 2, 1) NFLG = 1 GO TO 40 400 CALL XERROR(47H BSKIN, X AND N CANNOT BE ZERO AT THE SAME TIME, * 47, 2, 1) NFLG = 1 GO TO 50 410 CALL XSETF(II) CALL XERROR(34H BSKIN, END INPUT ERRORS FOR BSKIN, 34, 2, 1) RETURN END SUBROUTINE BKISR(X, N, SUM) BKI 10 C C BKISR COMPUTES REPEATED INTEGRALS OF THE K0 BESSEL FUNCTION C BY THE SERIES FOR N=0,1, AND 2. C INTEGER I, K, KK, KKN, K1, N, NP REAL AK, ATOL, BK, C, FK, FN, HX, HXS, POL, PR, SUM, TKP, TOL, * TRM, X, XLN REAL PSIXN, R1MACH DIMENSION C(2) DATA C(1), C(2) /1.57079632679489662E+00,1.0E0/ C TOL = AMAX1(R1MACH(4),1.0E-18) IF (X.LT.TOL) GO TO 50 PR = 1.0E0 POL = 0.0E0 IF (N.EQ.0) GO TO 20 DO 10 I=1,N POL = -POL*X + C(I) PR = PR*X/FLOAT(I) 10 CONTINUE 20 CONTINUE HX = X*0.5E0 HXS = HX*HX XLN = ALOG(HX) NP = N + 1 TKP = 3.0E0 FK = 2.0E0 FN = FLOAT(N) BK = 4.0E0 AK = 2.0E0/((FN+1.0E0)*(FN+2.0E0)) SUM = AK*(PSIXN(N+3)-PSIXN(3)+PSIXN(2)-XLN) ATOL = SUM*TOL*0.75E0 DO 30 K=2,20 AK = AK*(HXS/BK)*((TKP+1.0E0)/(TKP+FN+1.0E0))*(TKP/(TKP+FN)) K1 = K + 1 KK = K1 + K KKN = KK + N TRM = (PSIXN(K1)+PSIXN(KKN)-PSIXN(KK)-XLN)*AK SUM = SUM + TRM IF (ABS(TRM).LE.ATOL) GO TO 40 TKP = TKP + 2.0E0 BK = BK + TKP FK = FK + 1.0E0 30 CONTINUE GO TO 80 40 CONTINUE SUM = (SUM*HXS+PSIXN(NP)-XLN)*PR IF (N.EQ.1) SUM = -SUM SUM = POL + SUM RETURN C----------------------------------------------------------------------- C SMALL X CASE, X.LT.WORD TOLERANCE C----------------------------------------------------------------------- 50 CONTINUE IF (X.EQ.0.0E0 .AND. N.EQ.0) GO TO 70 IF (N.GT.0) GO TO 60 HX = X*0.5E0 SUM = PSIXN(1) - ALOG(HX) RETURN 60 CONTINUE SUM = C(N) RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 70 CALL XERROR(43H BKISR, X IS NOT GREATER THAN ZERO WHEN N=0, 43, * 2, 1) RETURN 80 CALL XERROR(50H BKISR, X IS TOO LARGE FOR CONVERGENCE IN 20 TERMS, * 50, 5, 1) RETURN END SUBROUTINE BKIAS(X, N, KTRMS, T, ANS, IND, MS, GMRN, H) BKI 10 C C BKIAS COMPUTES REPEATED INTEGRALS OF THE K0 BESSEL FUNCTION C BY THE ASYMPTOTIC EXPANSION C INTEGER I, II, IND, J, JMI, JN, K, KK, KM, KTRMS, MM, MP, MS, N REAL ANS, B, BND, DEN1, DEN2, DEN3, ER, ERR, FJ, FK, FLN, FM1, * GMRN, G1, GS, H, HN, HRTPI, RAT, RG1, RXP, RZ, RZX, S, SS, SUMI, * SUMJ, T, TOL, V, W, X, XP, Z REAL GAMRN, R1MACH DIMENSION B(120), XP(16), S(31), H(1), V(52), W(52), T(50), * BND(15) C----------------------------------------------------------------------- C COEFFICIENTS OF POLYNOMIAL P(J-1,X), J=1,15 C----------------------------------------------------------------------- DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19), * B(20), B(21), B(22), B(23), B(24) /1.00000000000000000E+00, * 1.00000000000000000E+00,-2.00000000000000000E+00, * 1.00000000000000000E+00,-8.00000000000000000E+00, * 6.00000000000000000E+00,1.00000000000000000E+00, * -2.20000000000000000E+01,5.80000000000000000E+01, * -2.40000000000000000E+01,1.00000000000000000E+00, * -5.20000000000000000E+01,3.28000000000000000E+02, * -4.44000000000000000E+02,1.20000000000000000E+02, * 1.00000000000000000E+00,-1.14000000000000000E+02, * 1.45200000000000000E+03,-4.40000000000000000E+03, * 3.70800000000000000E+03,-7.20000000000000000E+02, * 1.00000000000000000E+00,-2.40000000000000000E+02, * 5.61000000000000000E+03/ DATA B(25), B(26), B(27), B(28), B(29), B(30), B(31), B(32), * B(33), B(34), B(35), B(36), B(37), B(38), B(39), B(40), B(41), * B(42), B(43), B(44), B(45), B(46), B(47), B(48) * /-3.21200000000000000E+04,5.81400000000000000E+04, * -3.39840000000000000E+04,5.04000000000000000E+03, * 1.00000000000000000E+00,-4.94000000000000000E+02, * 1.99500000000000000E+04,-1.95800000000000000E+05, * 6.44020000000000000E+05,-7.85304000000000000E+05, * 3.41136000000000000E+05,-4.03200000000000000E+04, * 1.00000000000000000E+00,-1.00400000000000000E+03, * 6.72600000000000000E+04,-1.06250000000000000E+06, * 5.76550000000000000E+06,-1.24400640000000000E+07, * 1.10262960000000000E+07,-3.73392000000000000E+06, * 3.62880000000000000E+05,1.00000000000000000E+00, * -2.02600000000000000E+03,2.18848000000000000E+05/ DATA B(49), B(50), B(51), B(52), B(53), B(54), B(55), B(56), * B(57), B(58), B(59), B(60), B(61), B(62), B(63), B(64), B(65), * B(66), B(67), B(68), B(69), B(70), B(71), B(72) * /-5.32616000000000000E+06,4.47650000000000000E+07, * -1.55357384000000000E+08,2.38904904000000000E+08, * -1.62186912000000000E+08,4.43390400000000000E+07, * -3.62880000000000000E+06,1.00000000000000000E+00, * -4.07200000000000000E+03,6.95038000000000000E+05, * -2.52439040000000000E+07,3.14369720000000000E+08, * -1.64838430400000000E+09,4.00269508800000000E+09, * -4.64216395200000000E+09,2.50748121600000000E+09, * -5.68356480000000000E+08,3.99168000000000000E+07, * 1.00000000000000000E+00,-8.16600000000000000E+03, * 2.17062600000000000E+06,-1.14876376000000000E+08, * 2.05148277600000000E+09,-1.55489607840000000E+10/ DATA B(73), B(74), B(75), B(76), B(77), B(78), B(79), B(80), * B(81), B(82), B(83), B(84), B(85), B(86), B(87), B(88), B(89), * B(90), B(91), B(92), B(93), B(94), B(95), B(96) * /5.60413987840000000E+10,-1.01180433024000000E+11, * 9.21997902240000000E+10,-4.07883018240000000E+10, * 7.82771904000000000E+09,-4.79001600000000000E+08, * 1.00000000000000000E+00,-1.63560000000000000E+04, * 6.69969600000000000E+06,-5.07259276000000000E+08, * 1.26698177760000000E+10,-1.34323420224000000E+11, * 6.87720046384000000E+11,-1.81818864230400000E+12, * 2.54986547342400000E+12,-1.88307966182400000E+12, * 6.97929436800000000E+11,-1.15336085760000000E+11, * 6.22702080000000000E+09,1.00000000000000000E+00, * -3.27380000000000000E+04,2.05079880000000000E+07, * -2.18982980800000000E+09,7.50160522280000000E+10/ DATA B(97), B(98), B(99), B(100), B(101), B(102), B(103), B(104), * B(105), B(106), B(107), B(108), B(109), B(110), B(111), B(112), * B(113), B(114), B(115), B(116), B(117), B(118) * /-1.08467651241600000E+12,7.63483214939200000E+12, * -2.82999100661120000E+13,5.74943734645920000E+13, * -6.47283751398720000E+13,3.96895780558080000E+13, * -1.25509040179200000E+13,1.81099255680000000E+12, * -8.71782912000000000E+10,1.00000000000000000E+00, * -6.55040000000000000E+04,6.24078900000000000E+07, * -9.29252692000000000E+09,4.29826006340000000E+11, * -8.30844432796800000E+12,7.83913848313120000E+13, * -3.94365587815520000E+14,1.11174747256968000E+15, * -1.79717122069056000E+15,1.66642448627145600E+15, * -8.65023253219584000E+14,2.36908271543040000E+14/ DATA B(119), B(120) /-3.01963769856000000E+13, * 1.30767436800000000E+12/ C----------------------------------------------------------------------- C BOUNDS B(M,K) , K=M-3 C----------------------------------------------------------------------- DATA BND(1), BND(2), BND(3), BND(4), BND(5), BND(6), BND(7), * BND(8), BND(9), BND(10), BND(11), BND(12), BND(13), BND(14), * BND(15) /1.0E0,1.0E0,1.0E0,1.0E0,3.10E0,5.18E0,11.7E0,29.8E0, * 90.4E0,297.0E0,1070.0E0,4290.0E0,18100.0E0,84700.0E0,408000.0E0/ DATA HRTPI /8.86226925452758014E-01/ C TOL = AMAX1(R1MACH(4),1.0E-18) FLN = FLOAT(N) RZ = 1.0E0/(X+FLN) RZX = X*RZ Z = 0.5E0*(X+FLN) IF (IND.GT.1) GO TO 10 GMRN = GAMRN(Z) 10 CONTINUE GS = HRTPI*GMRN G1 = GS + GS RG1 = 1.0E0/G1 GMRN = (RZ+RZ)/GMRN IF (IND.GT.1) GO TO 70 C----------------------------------------------------------------------- C EVALUATE ERROR FOR M=MS C----------------------------------------------------------------------- HN = 0.5E0*FLN DEN2 = KTRMS + KTRMS + N DEN3 = DEN2 - 2.0E0 DEN1 = X + DEN2 ERR = RG1*(X+X)/(DEN1-1.0E0) IF (N.EQ.0) GO TO 20 RAT = 1.0E0/(FLN*FLN) 20 CONTINUE IF (KTRMS.EQ.0) GO TO 30 FJ = FLOAT(KTRMS) RAT = 0.25E0/(HRTPI*DEN3*SQRT(FJ)) 30 CONTINUE ERR = ERR*RAT FJ = -3.0E0 DO 50 J=1,15 IF (J.LE.5) ERR = ERR/DEN1 FM1 = AMAX1(1.0E0,FJ) FJ = FJ + 1.0E0 ER = BND(J)*ERR IF (KTRMS.EQ.0) GO TO 40 ER = ER/FM1 IF (ER.LT.TOL) GO TO 60 IF (J.GE.5) ERR = ERR/DEN3 GO TO 50 40 CONTINUE ER = ER*(1.0E0+HN/FM1) IF (ER.LT.TOL) GO TO 60 IF (J.GE.5) ERR = ERR/FLN 50 CONTINUE GO TO 200 60 CONTINUE MS = J 70 CONTINUE MM = MS + MS MP = MM + 1 C----------------------------------------------------------------------- C H(K)=(-Z)**(K)*(PSI(K-1,Z)-PSI(K-1,Z+0.5))/GAMMA(K) , K=1,2,...,MM C----------------------------------------------------------------------- IF (IND.GT.1) GO TO 80 CALL HKSEQ(Z, MM, H) GO TO 100 80 CONTINUE RAT = Z/(Z-0.5E0) RXP = RAT DO 90 I=1,MM H(I) = RXP*(1.0E0-H(I)) RXP = RXP*RAT 90 CONTINUE 100 CONTINUE C----------------------------------------------------------------------- C SCALED S SEQUENCE C----------------------------------------------------------------------- S(1) = 1.0E0 FK = 1.0E0 DO 120 K=2,MP SS = 0.0E0 KM = K - 1 I = KM DO 110 J=1,KM SS = SS + S(J)*H(I) I = I - 1 110 CONTINUE S(K) = SS/FK FK = FK + 1.0E0 120 CONTINUE C----------------------------------------------------------------------- C SCALED S-TILDA SEQUENCE C----------------------------------------------------------------------- IF (KTRMS.EQ.0) GO TO 160 FK = 0.0E0 SS = 0.0E0 RG1 = RG1/Z DO 130 K=1,KTRMS V(K) = Z/(Z+FK) W(K) = T(K)*V(K) SS = SS + W(K) FK = FK + 1.0E0 130 CONTINUE S(1) = S(1) - SS*RG1 DO 150 I=2,MP SS = 0.0E0 DO 140 K=1,KTRMS W(K) = W(K)*V(K) SS = SS + W(K) 140 CONTINUE S(I) = S(I) - SS*RG1 150 CONTINUE 160 CONTINUE C----------------------------------------------------------------------- C SUM ON J C----------------------------------------------------------------------- SUMJ = 0.0E0 JN = 1 RXP = 1.0E0 XP(1) = 1.0E0 DO 190 J=1,MS JN = JN + J - 1 XP(J+1) = XP(J)*RZX RXP = RXP*RZ C----------------------------------------------------------------------- C SUM ON I C----------------------------------------------------------------------- SUMI = 0.0E0 II = JN DO 180 I=1,J JMI = J - I + 1 KK = J + I + 1 DO 170 K=1,JMI V(K) = S(KK)*XP(K) KK = KK + 1 170 CONTINUE CALL BDIFF(JMI, V) SUMI = SUMI + B(II)*V(JMI)*XP(I+1) II = II + 1 180 CONTINUE SUMJ = SUMJ + SUMI*RXP 190 CONTINUE ANS = GS*(S(1)-SUMJ) RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 200 CALL XERROR(66H BKIAS, ERROR TEST NOT MET. INCREASE TOL IN BKIAS O *R NLIM IN BSKIN, 66, 5, 1) RETURN END SUBROUTINE BDIFF(L, V) BDI 10 C C BDIFF COMPUTES THE SUM OF B(L,K)*V(K)*(-1)**K WHERE B(L,K) C ARE THE BINOMIAL COEFFICIENTS. TRUNCATED SUMS ARE COMPUTED BY C SETTING LAST PART OF THE V VECTOR TO ZERO. ON RETURN, THE BINOMIAL C SUM IS IN V(L). C INTEGER I, J, K, L REAL V DIMENSION V(1) IF (L.EQ.1) RETURN DO 20 J=2,L K = L DO 10 I=J,L V(K) = V(K-1) - V(K) K = K - 1 10 CONTINUE 20 CONTINUE RETURN END FUNCTION GAMRN(X) GAM 10 C C WRITTEN BY D.E. AMOS, JANUARY, 1981. C C REFERENCES C THE SPECIAL FUNCTIONS AND THEIR APPROXIMATIONS, VOL. 1 BY C Y.L. LUKE, MATH IN SCI. AND ENG. SERIES 53, ACADEMIC PRESS, C NEW YORK, 1969, PP 34-35. C C ABSTRACT C GAMRN COMPUTES THE GAMMA FUNCTION RATIO GAMMA(X)/GAMMA(X+0.5) C FOR REAL X.GT.0. IF X.GE.XMIN, AN ASYMPTOTIC EXPANSION IS C EVALUATED. IF X.LT.XMIN, AN INTEGER IS ADDED TO X TO FORM A C NEW VALUE OF X.GE.XMIN AND THE ASYMPTOTIC EXPANSION IS EVAL- C UATED FOR THIS NEW VALUE OF X. SUCCESSIVE APPLICATION OF THE C RECURRENCE RELATION C C W(X)=W(X+1)*(1+0.5/X) C C REDUCES THE ARGUMENT TO ITS ORIGINAL VALUE. XMIN AND COMP- C UTATIONAL TOLERANCES ARE COMPUTED AS A FUNCTION OF THE NUMBER C OF DIGITS CARRIED IN A WORD BY CALLS TO I1MACH AND R1MACH. C HOWEVER, THE COMPUTATIONAL ACCURACY IS LIMITED TO THE MAX- C IMUM OF UNIT ROUNDOFF (=R1MACH(4)) AND 1.0E-18 SINCE CRITICAL C CONSTANTS ARE GIVEN TO ONLY 18 DIGITS. C C GAMRN CALLS R1MACH, I1MACH, XERROR. C C DESCRIPTION OF ARGUMENTS C C INPUT C X - ARGUMENT, X.GT.0.0 C C OUTPUT C GAMRN - RATIO GAMMA(X)/GAMMA(X+0.5) C C ERROR CONDITIONS C X.LE.0 IS A FATAL ERROR C INTEGER I, I1M11, K, MX, NX INTEGER I1MACH REAL FLN, GR, RLN, S, TOL, TRM, X, XDMY, XINC, XM, XMIN, XP, XSQ REAL R1MACH DIMENSION GR(12) C DATA GR(1), GR(2), GR(3), GR(4), GR(5), GR(6), GR(7), GR(8), * GR(9), GR(10), GR(11), GR(12) /1.00000000000000000E+00, * -1.56250000000000000E-02,2.56347656250000000E-03, * -1.27983093261718750E-03,1.34351104497909546E-03, * -2.43289663922041655E-03,6.75423753364157164E-03, * -2.66369606131178216E-02,1.41527455519564332E-01, * -9.74384543032201613E-01,8.43686251229783675E+00, * -8.97258321640552515E+01/ C IF (X.LE.0.0E0) GO TO 60 NX = INT(X) TOL = AMAX1(R1MACH(4),1.0E-18) I1M11 = I1MACH(11) RLN = R1MACH(5)*FLOAT(I1M11) FLN = AMIN1(RLN,20.0E0) FLN = AMAX1(FLN,3.0E0) FLN = FLN - 3.0E0 XM = 2.0E0 + FLN*(0.2366E0+0.01723E0*FLN) MX = INT(XM) + 1 XMIN = FLOAT(MX) XDMY = X - 0.25E0 XINC = 0.0E0 IF (X.GE.XMIN) GO TO 10 XINC = XMIN - FLOAT(NX) XDMY = XDMY + XINC 10 CONTINUE S = 1.0E0 IF (XDMY*TOL.GT.1.0E0) GO TO 30 XSQ = 1.0E0/(XDMY*XDMY) XP = XSQ DO 20 K=2,12 TRM = GR(K)*XP IF (ABS(TRM).LT.TOL) GO TO 30 S = S + TRM XP = XP*XSQ 20 CONTINUE 30 CONTINUE S = S/SQRT(XDMY) IF (XINC.NE.0.0E0) GO TO 40 GAMRN = S RETURN 40 CONTINUE NX = INT(XINC) XP = 0.0E0 DO 50 I=1,NX S = S*(1.0E0+0.5E0/(X+XP)) XP = XP + 1.0E0 50 CONTINUE GAMRN = S RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 60 CALL XERROR(34H GAMRN, X IS NOT GREATER THAN XERO, 34, 2, 1) RETURN END SUBROUTINE HKSEQ(X, M, H) HKS 10 C C HKSEQ IS AN ADAPTATION OF SUBROUTINE PSIFN, ACM TRANS. MATH. C SOFTWARE, 1983 TO GENERATE H(K,X)=(-X)**(K+1)*(PSI(K,X)- C -PSI(K,X+0.5))/GAMMA(K+1) FOR K=0,...,M. C INTEGER I, J, K, M, MX, NX INTEGER I1MACH REAL B, FK, FLN, FN, FNP, H, HRX, RLN, RXSQ, R1M5, S, SLOPE, T, * TK, TRM, TRMH, TRMR, TST, U, V, WDTOL, X, XDMY, XH, XINC, XM, * XMIN, YINT REAL R1MACH DIMENSION B(22), TRM(22), TRMR(25), TRMH(25), U(25), V(25), H(1) C----------------------------------------------------------------------- C SCALED BERNOULLI NUMBERS 2.0*B(2K)*(1-2**(-2K)) C----------------------------------------------------------------------- DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19), * B(20), B(21), B(22) /1.00000000000000000E+00, * -5.00000000000000000E-01,2.50000000000000000E-01, * -6.25000000000000000E-02,4.68750000000000000E-02, * -6.64062500000000000E-02,1.51367187500000000E-01, * -5.06103515625000000E-01,2.33319091796875000E+00, * -1.41840972900390625E+01,1.09941936492919922E+02, * -1.05824747562408447E+03,1.23842434241771698E+04, * -1.73160495905935764E+05,2.85103429084961116E+06, * -5.45964619322445132E+07,1.20316174668075304E+09, * -3.02326315271452307E+10,8.59229286072319606E+11, * -2.74233104097776039E+13,9.76664637943633248E+14, * -3.85931586838450360E+16/ C C WDTOL = AMAX1(R1MACH(4),1.0E-18) FN = FLOAT(M-1) FNP = FN + 1.0E0 C----------------------------------------------------------------------- C COMPUTE XMIN C----------------------------------------------------------------------- R1M5 = R1MACH(5) RLN = R1M5*FLOAT(I1MACH(11)) RLN = AMIN1(RLN,18.06E0) FLN = AMAX1(RLN,3.0E0) - 3.0E0 YINT = 3.50E0 + 0.40E0*FLN SLOPE = 0.21E0 + FLN*(0.0006038E0*FLN+0.008677E0) XM = YINT + SLOPE*FN MX = INT(XM) + 1 XMIN = FLOAT(MX) C----------------------------------------------------------------------- C GENERATE H(M-1,XDMY)*XDMY**(M) BY THE ASYMPTOTIC EXPANSION C----------------------------------------------------------------------- XDMY = X XINC = 0.0E0 IF (X.GE.XMIN) GO TO 10 NX = INT(X) XINC = XMIN - FLOAT(NX) XDMY = X + XINC 10 CONTINUE RXSQ = 1.0E0/(XDMY*XDMY) HRX = 0.5E0/XDMY TST = 0.5E0*WDTOL T = FNP*HRX C----------------------------------------------------------------------- C INITIALIZE COEFFICIENT ARRAY C----------------------------------------------------------------------- S = T*B(3) IF (ABS(S).LT.TST) GO TO 30 TK = 2.0E0 DO 20 K=4,22 T = T*((TK+FN+1.0E0)/(TK+1.0E0))*((TK+FN)/(TK+2.0E0))*RXSQ TRM(K) = T*B(K) IF (ABS(TRM(K)).LT.TST) GO TO 30 S = S + TRM(K) TK = TK + 2.0E0 20 CONTINUE GO TO 110 30 CONTINUE H(M) = S + 0.5E0 IF (M.EQ.1) GO TO 70 C----------------------------------------------------------------------- C GENERATE LOWER DERIVATIVES, I.LT.M-1 C----------------------------------------------------------------------- DO 60 I=2,M FNP = FN FN = FN - 1.0E0 S = FNP*HRX*B(3) IF (ABS(S).LT.TST) GO TO 50 FK = FNP + 3.0E0 DO 40 K=4,22 TRM(K) = TRM(K)*FNP/FK IF (ABS(TRM(K)).LT.TST) GO TO 50 S = S + TRM(K) FK = FK + 2.0E0 40 CONTINUE GO TO 110 50 CONTINUE MX = M - I + 1 H(MX) = S + 0.5E0 60 CONTINUE 70 CONTINUE IF (XINC.EQ.0.0E0) RETURN C----------------------------------------------------------------------- C RECUR BACKWARD FROM XDMY TO X C----------------------------------------------------------------------- XH = X + 0.5E0 S = 0.0E0 NX = INT(XINC) DO 80 I=1,NX TRMR(I) = X/(X+FLOAT(NX-I)) U(I) = TRMR(I) TRMH(I) = X/(XH+FLOAT(NX-I)) V(I) = TRMH(I) S = S + U(I) - V(I) 80 CONTINUE MX = NX + 1 TRMR(MX) = X/XDMY U(MX) = TRMR(MX) H(1) = H(1)*TRMR(MX) + S IF (M.EQ.1) RETURN DO 100 J=2,M S = 0.0E0 DO 90 I=1,NX TRMR(I) = TRMR(I)*U(I) TRMH(I) = TRMH(I)*V(I) S = S + TRMR(I) - TRMH(I) 90 CONTINUE TRMR(MX) = TRMR(MX)*U(MX) H(J) = H(J)*TRMR(MX) + S 100 CONTINUE RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 110 CALL XERROR(44H HKSEQ, ERROR TEST IN ASYMPTOTIC SUM NOT MET, 44, * 2, 1) RETURN END SUBROUTINE EXINT(X, N, KODE, M, TOL, EN, IERR) EXI 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, 6, 1980, PP. 365-377, 420-428 C C ABSTRACT C EXINT 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 EXINT CALLS I1MACH, R1MACH, PSIXN, XERROR 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 = SINGLE PRECISION UNIT ROUNDOFF = R1MACH(4) 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.0E0 , K=1,M RETURNED ON KODE=1 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***END PROLOGUE C C REAL A, AA, AAMS, AH, AK, AT, B, BK, BT, CC, CNORM, CT, EM, EMX, * EN, ETOL, FNM, FX, PT, P1, P2, S, TOL, TX, X, XCUT, XLIM, XTOL, * Y, YT, Y1, Y2 REAL R1MACH, PSIXN INTEGER I, IC, ICASE, ICT, IERR, IK, IND, IX, I1M, JSET, K, KK, * KN, KODE, KS, M, ML, MU, N, ND, NM INTEGER I1MACH DIMENSION EN(1), A(99), B(99), Y(2) C ETOL = AMAX1(R1MACH(4),0.5E-18) I1M = -I1MACH(12) PT = 2.3026E0*R1MACH(5)*FLOAT(I1M) XLIM = PT - 6.907755E0 BT = PT + FLOAT(N+M-1) IF (BT.GT.1000.0E0) XLIM = PT - ALOG(BT) 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 XCUT = 2.0E0 IF (ETOL.GT.2.0E-7) XCUT = 1.0E0 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----------------------------------------------------------------------- TX = X + 0.5E0 IX = INT(TX) C----------------------------------------------------------------------- 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 C----------------------------------------------------------------------- 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 IC = 35 IF (X.LT.ETOL) IC = 1 DO 50 I=1,IC 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)+PSIXN(ND)) XTOL = 3.0E0*TOL 40 AK = AK + 1.0E0 50 CONTINUE IF (IC.NE.1) GO TO 320 60 IF (ND.EQ.1) S = S + (-ALOG(X)+PSIXN(1)) 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----------------------------------------------------------------------- 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. C----------------------------------------------------------------------- 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----------------------------------------------------------------------- C FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD C RECURSION C----------------------------------------------------------------------- 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----------------------------------------------------------------------- C BACKWARD RECURRENCE FOR C Y1= C*U( A ,A,X) C Y2= C*(A/(1+A/2))*U(A+1,A,X) C----------------------------------------------------------------------- DO 200 K=1,ICT KK = KK - 1 YT = Y1 Y1 = (B(KK)*Y1-Y2)/A(KK) Y2 = YT 200 CONTINUE C----------------------------------------------------------------------- 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 C----------------------------------------------------------------------- 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 ERROR MESSAGES C----------------------------------------------------------------------- 260 CALL XERROR(31HIN EXINT, N NOT GREATER THAN 0, 31, 2, 1) RETURN 270 CALL XERROR(26HIN EXINT, KODE NOT 1 OR 2, 26, 2, 1) RETURN 280 CALL XERROR(31HIN EXINT, M NOT GREATER THAN 0, 31, 2, 1) RETURN 290 CALL XERROR(32HIN EXINT, TOL NOT WITHIN LIMITS, 32, 2, 1) RETURN 300 CALL XERROR(36HIN EXINT, X IS NOT ZERO OR POSITIVE, 36, 2, 1) RETURN 310 CALL XERROR(66HIN EXINT, THE EXPONENTIAL INTEGRAL IS NOT DEFINED *FOR X=0 AND N=1, 66, 2, 1) RETURN 320 CALL XERROR(73HIN EXINT, RELATIVE ERROR TEST FOR SERIES TERMINATI *ON NOT MET IN 36 TERMS, 73, 3, 1) RETURN 330 CALL XERROR(68HIN EXINT, TERMINATION TEST FOR MILLER ALGORITHM NO *T MET IN 99 STEPS, 68, 3, 1) RETURN 340 CALL XERROR(92HIN EXINT, AN ERROR IN PLACING INT(X+0.5) WITH RESP *ECT TO N AND N+M-1 OCCURRED FOR X.GT.XCUT, 92, 8, 1) RETURN END FUNCTION PSIXN(N) PSI 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 C INTEGER N, K REAL AX, B, C, FN, RFN2, TRM, S, WDTOL REAL R1MACH DIMENSION B(6), C(100) C----------------------------------------------------------------------- C PSIXN(N), N = 1,100 C----------------------------------------------------------------------- DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), 2 C(19), C(20), C(21), C(22), C(23), C(24)/ 3 -5.77215664901532861E-01, 4.22784335098467139E-01, 4 9.22784335098467139E-01, 1.25611766843180047E+00, 5 1.50611766843180047E+00, 1.70611766843180047E+00, 6 1.87278433509846714E+00, 2.01564147795561000E+00, 7 2.14064147795561000E+00, 2.25175258906672111E+00, 8 2.35175258906672111E+00, 2.44266167997581202E+00, 9 2.52599501330914535E+00, 2.60291809023222227E+00, A 2.67434666166079370E+00, 2.74101332832746037E+00, B 2.80351332832746037E+00, 2.86233685773922507E+00, C 2.91789241329478063E+00, 2.97052399224214905E+00, D 3.02052399224214905E+00, 3.06814303986119667E+00, E 3.11359758531574212E+00, 3.15707584618530734E+00/ DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ 3 3.19874251285197401E+00, 3.23874251285197401E+00, 4 3.27720405131351247E+00, 3.31424108835054951E+00, 5 3.34995537406483522E+00, 3.38443813268552488E+00, 6 3.41777146601885821E+00, 3.45002953053498724E+00, 7 3.48127953053498724E+00, 3.51158256083801755E+00, 8 3.54099432554389990E+00, 3.56956575411532847E+00, 9 3.59734353189310625E+00, 3.62437055892013327E+00, A 3.65068634839381748E+00, 3.67632737403484313E+00, B 3.70132737403484313E+00, 3.72571761793728215E+00, C 3.74952714174680596E+00, 3.77278295570029433E+00, D 3.79551022842756706E+00, 3.81773245064978928E+00, E 3.83947158108457189E+00, 3.86074817682925274E+00/ DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), 2 C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ 3 3.88158151016258607E+00, 3.90198967342789220E+00, 4 3.92198967342789220E+00, 3.94159751656514710E+00, 5 3.96082828579591633E+00, 3.97969621032421822E+00, 6 3.99821472884273674E+00, 4.01639654702455492E+00, 7 4.03425368988169777E+00, 4.05179754953082058E+00, 8 4.06903892884116541E+00, 4.08598808138353829E+00, 9 4.10265474805020496E+00, 4.11904819067315578E+00, A 4.13517722293122029E+00, 4.15105023880423617E+00, B 4.16667523880423617E+00, 4.18205985418885155E+00, C 4.19721136934036670E+00, 4.21213674247469506E+00, D 4.22684262482763624E+00, 4.24133537845082464E+00, E 4.25562109273653893E+00, 4.26970559977879245E+00/ DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ 3 4.28359448866768134E+00, 4.29729311880466764E+00, 4 4.31080663231818115E+00, 4.32413996565151449E+00, 5 4.33729786038835659E+00, 4.35028487337536958E+00, 6 4.36310538619588240E+00, 4.37576361404398366E+00, 7 4.38826361404398366E+00, 4.40060929305632934E+00, 8 4.41280441500754886E+00, 4.42485260777863319E+00, 9 4.43675736968339510E+00, 4.44852207556574804E+00, A 4.46014998254249223E+00, 4.47164423541605544E+00, B 4.48300787177969181E+00, 4.49424382683587158E+00, C 4.50535493794698269E+00, 4.51634394893599368E+00, D 4.52721351415338499E+00, 4.53796620232542800E+00, E 4.54860450019776842E+00, 4.55913081598724211E+00/ DATA C(97), C(98), C(99), C(100)/ 1 4.56954748265390877E+00, 4.57985676100442424E+00, 2 4.59006084263707730E+00, 4.60016185273808740E+00/ C----------------------------------------------------------------------- C COEFFICIENTS OF ASYMPTOTIC EXPANSION C----------------------------------------------------------------------- DATA B(1), B(2), B(3), B(4), B(5), B(6)/ 1 8.33333333333333333E-02, -8.33333333333333333E-03, 2 3.96825396825396825E-03, -4.16666666666666666E-03, 3 7.57575757575757576E-03, -2.10927960927960928E-02/ C IF (N.GT.100) GO TO 10 PSIXN = C(N) RETURN 10 CONTINUE WDTOL = AMAX1(R1MACH(4),1.0E-18) FN = FLOAT(N) AX = 1.0E0 S = -0.5E0/FN IF (ABS(S).LE.WDTOL) GO TO 30 RFN2 = 1.0E0/(FN*FN) DO 20 K=1,6 AX = AX*RFN2 TRM = -B(K)*AX IF (ABS(TRM).LT.WDTOL) GO TO 30 S = S + TRM 20 CONTINUE 30 CONTINUE PSIXN = S + ALOG(FN) RETURN END C PROGRAM DQCKIN(INPUT,OUTPUT,TAPE6=OUTPUT) MAN 10 C MAN 20 C ABSTRACT MAN 30 C PROGRAM DQCKIN IS A QUICK CHECK DRIVER WHICH EXERCISES THE MAJOR MAN 40 C LOOPS IN SUBROUTINE DBSKIN (X,N,KODE,M,Y,IERR) FOR BICKLEY MAN 50 C FUNCTIONS KI(J,X). MORE PRECISELY, DQCKIN DOES CONSISTENCY CHECKSMAN 60 C ON THE OUTPUT FROM DBSKIN BY COMPARING SINGLE EVALUATIONS (M=1) MAN 70 C AGAINST SELECTED MEMBERS OF SEQUENCES WHICH ARE GENERATED BY MAN 80 C RECURSION. IF THE RELATIVE ERROR IS LESS THAN 1000 TIMES UNIT MAN 90 C ROUND OFF, THEN THE TEST IS PASSED - IF NOT, THEN X, THE VALUES MAN 100 C TO BE COMPARED, THE RELATIVE ERROR AND PARAMETERS KODE, N, M AND KMAN 110 C ARE WRITTEN ON LOGICAL UNIT 6 WHERE K IS THE MEMBER OF THE MAN 120 C SEQUENCE OF LENGTH M WHICH FAILED THE TEST. THAT IS, THE INDEX MAN 130 C OF THE FUNCTION WHICH FAILED THE TEST IS J=N+K-1. UNDERFLOW MAN 140 C TESTS ARE MADE AND ERROR CONDITIONS ARE TRIGGERED. MAN 150 C MAN 160 C FUNCTIONS I1MACH AND D1MACH MUST BE INITIALIZED ACCORDING TO THE MAN 170 C PROLOGUE IN EACH FUNCTION FOR THE MACHINE ENVIRONMENT BEFORE MAN 180 C DQCKIN OR DBSKIN CAN BE EXECUTED. FIFTEEN MACHINE ENVIRONMENTS MAN 190 C CAN BE DEFINED IN I1MACH AND D1MACH. MAN 200 C MAN 210 INTEGER I, IERR, IFLG, IX, I1M12, J, K, KODE, LUN, M, MDEL, MM, MAN 220 * N, NDEL, NN MAN 230 INTEGER I1MACH MAN 240 DOUBLE PRECISION AIX, ER, TOL, V, X, XINC, Y MAN 250 DOUBLE PRECISION D1MACH MAN 260 DIMENSION V(1), Y(10) MAN 270 DATA LUN /6/ MAN 280 C-----------------------------------------------------------------------MAN 290 C CALLS TO SLATEC PACKAGE TO OVERRIDE FATAL ERRORS AND WRITE MAN 300 C MESSSAGES TO UNIT LUN MAN 310 C CALL XSETF(-1) MAN 320 C CALL XSETUN(LUN) MAN 330 C-----------------------------------------------------------------------MAN 340 CALL XSETF(-1) MAN 350 CALL XSETUN(LUN) MAN 360 TOL = 1000.0D0*DMAX1(D1MACH(4),1.0D-18) MAN 370 IFLG = 0 MAN 380 WRITE (LUN,99999) MAN 390 99999 FORMAT (1H1//35H QUICK CHECK DIAGNOSTICS FOR DBSKIN//) MAN 400 DO 70 KODE=1,2 MAN 410 N = 0 MAN 420 DO 60 NN=1,7 MAN 430 M = 1 MAN 440 DO 50 MM=1,4 MAN 450 X = 0.0D0 MAN 460 DO 40 IX=1,6 MAN 470 IF (N.EQ.0 .AND. IX.EQ.1) GO TO 30 MAN 480 CALL DBSKIN(X, N, KODE, M, Y, IERR) MAN 490 DO 20 K=1,M,2 MAN 500 J = N + K - 1 MAN 510 CALL DBSKIN(X, J, KODE, 1, V, IERR) MAN 520 ER = DABS((V(1)-Y(K))/V(1)) MAN 530 IF (ER.LE.TOL) GO TO 20 MAN 540 IF (IFLG.NE.0) GO TO 10 MAN 550 WRITE (LUN,99998) MAN 560 99998 FORMAT (8X, 1HX, 13X, 4HV(1), 11X, 4HY(K), 9X, 6HREL ER,MAN 570 * 1HR, 5X, 4HKODE, 3X, 1HN, 4X, 1HM, 4X, 1HK) MAN 580 10 CONTINUE MAN 590 IFLG = IFLG + 1 MAN 600 WRITE (LUN,99997) X, V(1), Y(K), ER, KODE, N, M, K MAN 610 99997 FORMAT (4E15.6, 4I5) MAN 620 IF (IFLG.GT.200) GO TO 130 MAN 630 20 CONTINUE MAN 640 30 CONTINUE MAN 650 AIX = DBLE(FLOAT(2*IX-3)) MAN 660 XINC = DMAX1(1.0D0,AIX) MAN 670 X = X + XINC MAN 680 40 CONTINUE MAN 690 MDEL = MAX0(1,MM-1) MAN 700 M = M + MDEL MAN 710 50 CONTINUE MAN 720 NDEL = MAX0(1,2*N-2) MAN 730 N = N + NDEL MAN 740 60 CONTINUE MAN 750 70 CONTINUE MAN 760 C-----------------------------------------------------------------------MAN 770 C TEST UNDERFLOW MAN 780 C-----------------------------------------------------------------------MAN 790 KODE = 1 MAN 800 M = 10 MAN 810 N = 10 MAN 820 I1M12 = I1MACH(15) MAN 830 X = -2.302D0*D1MACH(5)*DBLE(FLOAT(I1M12)) MAN 840 CALL DBSKIN(X, N, KODE, M, Y, IERR) MAN 850 IF (IERR.EQ.1) GO TO 80 MAN 860 WRITE (LUN,99996) MAN 870 99996 FORMAT (//32H IERR IN UNDERFLOW TEST IS NOT 1//) MAN 880 IFLG = IFLG + 1 MAN 890 GO TO 110 MAN 900 80 CONTINUE MAN 910 DO 90 I=1,M MAN 920 IF (Y(I).NE.0.0D0) GO TO 100 MAN 930 90 CONTINUE MAN 940 GO TO 110 MAN 950 100 CONTINUE MAN 960 IFLG = IFLG + 1 MAN 970 WRITE (LUN,99995) MAN 980 99995 FORMAT (//43H SOME Y VALUE IN UNDERFLOW TEST IS NOT ZERO//) MAN 990 110 CONTINUE MAN 1000 IF (IFLG.NE.0) GO TO 120 MAN 1010 WRITE (LUN,99994) MAN 1020 99994 FORMAT (//16H QUICK CHECKS OK//) MAN 1030 120 CONTINUE MAN 1040 C-----------------------------------------------------------------------MAN 1050 C TRIGGER ERROR DIAGNOSTICS MAN 1060 C-----------------------------------------------------------------------MAN 1070 WRITE (LUN,99993) MAN 1080 99993 FORMAT (//27H TRIGGER 4 ERROR CONDITIONS//) MAN 1090 X = -1.0D0 MAN 1100 N = -1 MAN 1110 KODE = -1 MAN 1120 M = -1 MAN 1130 CALL DBSKIN(X, N, KODE, M, Y, IERR) MAN 1140 STOP MAN 1150 130 CONTINUE MAN 1160 WRITE (LUN,99992) MAN 1170 99992 FORMAT (//52H PROCESSING OF MAIN LOOPS TERMINATED BECAUSE THE NUM,MAN 1180 * 36HBER OF DIAGNOSTIC PRINTS EXCEEDS 200//) MAN 1190 GO TO 120 MAN 1200 END MAN 1210 SUBROUTINE DBSKIN(X, N, KODE, M, Y, IERR) DBS 10 C C WRITTEN BY D.E. AMOS, JUNE, 1982 C C REFERENCES C UNIFORM ASYMPTOTIC EXPANSIONS FOR EXPONENTIAL INTEGRALS E(N,X) C AND BICKLEY FUNCTIONS KI(N,X) BY D.E. AMOS, ACM TRANS. MATH C SOFTWARE, 1983 C C COMPUTATION OF BICKLEY FUNCTIONS BY D.E. AMOS, ACM TRANS. C MATH SOFTWARE, 1983 C C ABSTRACT ***A DOUBLE PRECISION ROUTINE*** C DBSKIN COMPUTES M MEMBER SEQUENCES OF BICKLEY FUNCTIONS, C REPEATED INTEGRALS OF THE K0 BESSEL FUNCTION, Y(K)=KI(N+K-1,X) C ON KODE=1 OR SCALED FUNCTIONS Y(K)=EXP(X)*KI(N+K-1,X) ON C KODE=2, K=1,...,M, FOR N.GE.0 AND X.GE.0 (N AND X CANNOT BE C ZERO SIMULTANEOUSLY). C C NUMERICAL RECURRENCE ON C C (L-1)KI(L,X)=X(KI(L-3,X)-KI(L-1,X))+(L-2)KI(L-2,X) C C IS STABLE WHERE RECURRENCE IS CARRIED FORWARD OR BACKWARD C AWAY FROM INT(X+0.5). THE POWER SERIES FOR INDICES 0,1 AND 2 C ON 0.LE.X.LE.2 STARTS A STABLE RECURRENCE FOR INDICES C GREATER THAN 2. IF N IS SUFFICIENTLY LARGE (N.GT.NLIM), THE C UNIFORM ASYMPTOTIC EXPANSION FOR N TO INFINITY IS MORE C ECONOMICAL. ON X.GT.2 THE RECURSION IS STARTED BY EVALUATING C THE UNIFORM EXPANSION FOR THE THREE MEMBERS WHOSE INDICES ARE C CLOSEST TO INT(X+0.5) WITHIN THE SET N,...,N+M-1. FORWARD C RECURRENCE, BACKWARD RECURRENCE OR BOTH COMPLETE THE C SEQUENCE DEPENDING ON THE RELATION OF INT(X+0.5) TO THE C INDICES N,...,N+M-1. C C FUNCTIONS I1MACH AND D1MACH MUST BE INITIALIZED ACCORDING C TO PROLOGUE INSTRUCTIONS. FIFTEEN MACHINE ENVIRONMENTS C CAN BE DEFINED BY MAKING APPROPRIATE COMMENT CARDS ACTIVE C IN I1MACH AND D1MACH. C C THE NOMINAL COMPUTATIONAL ACCURACY IS THE MAXIMUM OF UNIT C ROUNDOFF (=D1MACH(4)) AND 1.0D-18 SINCE CRITICAL CONSTANTS C ARE GIVEN TO ONLY 18 DIGITS. C C DBSKIN USES DBKIAS,DBKISR,DEXINT,DGAMRN,DPSIXN,DHKSEQ,DBDIFF, C I1MACH,D1MACH,XERROR PACKAGE(18 ROUTINES) C C BSKIN IS THE SINGLE PRECISION VERSION OF DBSKIN. C C DESCRIPTION OF ARGUMENTS C INPUT X IS DOUBLE PRECISION C X - ARGUMENT, X.GE.0.0D0 C N - ORDER OF FIRST MEMBER OF THE SEQUENCE N.GE.0 C KODE - SELECTION PARAMETER C KODE=1 RETURNS Y(K)= KI(N+K-1,X), K=1,M C =2 RETURNS Y(K)=EXP(X)*KI(N+K-1,X), K=1,M C M - NUMBER OF MEMBERS IN THE SEQUENCE, M.GE.1 C C OUTPUT Y IS DOUBLE PRECISION C Y - A VECTOR OF DIMENSION AT LEAST M CONTAINING THE C SEQUENCE SELECTED BY KODE C IERR - UNDERFLOW FLAG C IERR=0 MEANS COMPUTATION COMPLETED C =1 MEANS AN EXPONENTIAL UNDERFLOW OCCURRED ON C KODE=1 AND Y(K)=0.0D0, K=1,...,M IS RETURNED C C ERROR CONDITIONS C AN IMPROPER INPUT IS A FATAL ERROR C***END PROLOGUE INTEGER I, ICASE, IERR, II, IL, ISET, I1M, K, KK, KODE, KTRMS, M, * M3, N, NE, NFLG, NL, NLIM, NN, NP, NS, NT INTEGER I1MACH DOUBLE PRECISION A, ENLIM, EXI, FN, GR, H, HN, HRTPI, SS, TOL, * T1, T2, W, X, XLIM, XNLIM, XP, Y, YS, YSS DOUBLE PRECISION DGAMRN, D1MACH DIMENSION EXI(102), A(50), YS(3), YSS(3), H(31), Y(1) C----------------------------------------------------------------------- C COEFFICIENTS IN SERIES OF EXPONENTIAL INTEGRALS C----------------------------------------------------------------------- DATA A(1), A(2), A(3), A(4), A(5), A(6), A(7), A(8), A(9), A(10), * A(11), A(12), A(13), A(14), A(15), A(16), A(17), A(18), A(19), * A(20), A(21), A(22), A(23), A(24) /1.00000000000000000D+00, * 5.00000000000000000D-01,3.75000000000000000D-01, * 3.12500000000000000D-01,2.73437500000000000D-01, * 2.46093750000000000D-01,2.25585937500000000D-01, * 2.09472656250000000D-01,1.96380615234375000D-01, * 1.85470581054687500D-01,1.76197052001953125D-01, * 1.68188095092773438D-01,1.61180257797241211D-01, * 1.54981017112731934D-01,1.49445980787277222D-01, * 1.44464448094367981D-01,1.39949934091418982D-01, * 1.35833759559318423D-01,1.32060599571559578D-01, * 1.28585320635465905D-01,1.25370687619579257D-01, * 1.22385671247684513D-01,1.19604178719328047D-01, * 1.17004087877603524D-01/ DATA A(25), A(26), A(27), A(28), A(29), A(30), A(31), A(32), * A(33), A(34), A(35), A(36), A(37), A(38), A(39), A(40), A(41), * A(42), A(43), A(44), A(45), A(46), A(47), A(48) * /1.14566502713486784D-01,1.12275172659217048D-01, * 1.10116034723462874D-01,1.08076848895250599D-01, * 1.06146905164978267D-01,1.04316786110409676D-01, * 1.02578173008569515D-01,1.00923686347140974D-01, * 9.93467537479668965D-02,9.78414999033007314D-02, * 9.64026543164874854D-02,9.50254735405376642D-02, * 9.37056752969190855D-02,9.24393823875012600D-02, * 9.12230747245078224D-02,9.00535481254756708D-02, * 8.89278787739072249D-02,8.78433924473961612D-02, * 8.67976377754033498D-02,8.57883629175498224D-02, * 8.48134951571231199D-02,8.38711229887106408D-02, * 8.29594803475290034D-02,8.20769326842574183D-02/ DATA A(49), A(50) /8.12219646354630702D-02,8.03931690779583449D-02 * / C----------------------------------------------------------------------- C DSQRT(PI)/2 C----------------------------------------------------------------------- DATA HRTPI /8.86226925452758014D-01/ C NFLG = 0 ISET = 1 IF (X.LT.0.0D0) GO TO 340 10 CONTINUE ISET = 2 IF (N.LT.0) GO TO 340 20 CONTINUE ISET = 3 IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 340 30 CONTINUE ISET = 4 IF (M.LT.1) GO TO 340 40 CONTINUE ISET = 5 IF (X.EQ.0.0D0 .AND. N.EQ.0) GO TO 340 50 CONTINUE IF (NFLG.NE.0) GO TO 410 IERR = 0 IF (X.EQ.0.0D0) GO TO 300 I1M = -I1MACH(15) T1 = 2.3026D0*D1MACH(5)*DBLE(FLOAT(I1M)) XLIM = T1 - 3.228086D0 T2 = T1 + DBLE(FLOAT(N+M-1)) IF (T2.GT.1000.0D0) XLIM = T1 - 0.5D0*(DLOG(T2)-0.451583D0) IF (X.GT.XLIM .AND. KODE.EQ.1) GO TO 320 TOL = DMAX1(D1MACH(4),1.0D-18) I1M = I1MACH(14) C----------------------------------------------------------------------- C LN(NLIM) = 0.125*LN(EPS), NLIM = 2*KTRMS+N C----------------------------------------------------------------------- XNLIM = 0.287823D0*DBLE(FLOAT(I1M-1))*D1MACH(5) ENLIM = DEXP(XNLIM) NLIM = INT(SNGL(ENLIM)) + 2 NLIM = MIN0(100,NLIM) NLIM = MAX0(20,NLIM) M3 = MIN0(M,3) NL = N + M - 1 IF (X.GT.2.0D0) GO TO 130 IF (N.GT.NLIM) GO TO 280 C----------------------------------------------------------------------- C COMPUTATION BY SERIES FOR 0.LE.X.LE.2 C----------------------------------------------------------------------- NFLG = 0 NN = N IF (NL.LE.2) GO TO 60 M3 = 3 NN = 0 NFLG = 1 60 CONTINUE XP = 1.0D0 IF (KODE.EQ.2) XP = DEXP(X) DO 80 I=1,M3 CALL DBKISR(X, NN, W) W = W*XP IF (NN.LT.N) GO TO 70 KK = NN - N + 1 Y(KK) = W 70 CONTINUE YS(I) = W NN = NN + 1 80 CONTINUE IF (NFLG.EQ.0) RETURN NS = NN XP = 1.0D0 90 CONTINUE C----------------------------------------------------------------------- C FORWARD RECURSION SCALED BY DEXP(X) ON ICASE=0,1,2 C----------------------------------------------------------------------- FN = DBLE(FLOAT(NS-1)) IL = NL - NS + 1 IF (IL.LE.0) GO TO 120 DO 110 I=1,IL T1 = YS(2) T2 = YS(3) YS(3) = (X*(YS(1)-YS(3))+(FN-1.0D0)*YS(2))/FN YS(2) = T2 YS(1) = T1 FN = FN + 1.0D0 IF (NS.LT.N) GO TO 100 KK = NS - N + 1 Y(KK) = YS(3)*XP 100 CONTINUE NS = NS + 1 110 CONTINUE 120 CONTINUE RETURN C----------------------------------------------------------------------- C COMPUTATION BY ASYMPTOTIC EXPANSION FOR X.GT.2 C----------------------------------------------------------------------- 130 CONTINUE W = X + 0.5D0 NT = INT(SNGL(W)) IF (NL.GT.NT) GO TO 270 C----------------------------------------------------------------------- C CASE NL.LE.NT, ICASE=0 C----------------------------------------------------------------------- ICASE = 0 NN = NL NFLG = MIN0(M-M3,1) 140 CONTINUE KK = (NLIM-NN)/2 KTRMS = MAX0(0,KK) NS = NN + 1 NP = NN - M3 + 1 XP = 1.0D0 IF (KODE.EQ.1) XP = DEXP(-X) DO 150 I=1,M3 KK = I CALL DBKIAS(X, NP, KTRMS, A, W, KK, NE, GR, H) YS(I) = W NP = NP + 1 150 CONTINUE C----------------------------------------------------------------------- C SUM SERIES OF EXPONENTIAL INTEGRALS BACKWARD C----------------------------------------------------------------------- IF (KTRMS.EQ.0) GO TO 160 NE = KTRMS + KTRMS + 1 NP = NN - M3 + 2 CALL DEXINT(X, NP, 2, NE, TOL, EXI, IERR) IF (IERR.EQ.1) GO TO 320 160 CONTINUE DO 190 I=1,M3 SS = 0.0D0 IF (KTRMS.EQ.0) GO TO 180 KK = I + KTRMS + KTRMS - 2 IL = KTRMS DO 170 K=1,KTRMS SS = SS + A(IL)*EXI(KK) KK = KK - 2 IL = IL - 1 170 CONTINUE 180 CONTINUE YS(I) = YS(I) + SS 190 CONTINUE IF (ICASE.EQ.1) GO TO 200 IF (NFLG.NE.0) GO TO 220 200 CONTINUE DO 210 I=1,M3 Y(I) = YS(I)*XP 210 CONTINUE IF (ICASE.EQ.1 .AND. NFLG.EQ.1) GO TO 90 RETURN 220 CONTINUE C----------------------------------------------------------------------- C BACKWARD RECURSION SCALED BY DEXP(X) ICASE=0,2 C----------------------------------------------------------------------- KK = NN - N + 1 K = M3 DO 230 I=1,M3 Y(KK) = YS(K)*XP YSS(I) = YS(I) KK = KK - 1 K = K - 1 230 CONTINUE IL = KK IF (IL.LE.0) GO TO 250 FN = DBLE(FLOAT(NN-3)) DO 240 I=1,IL T1 = YS(2) T2 = YS(1) YS(1) = YS(2) + ((FN+2.0D0)*YS(3)-(FN+1.0D0)*YS(1))/X YS(2) = T2 YS(3) = T1 Y(KK) = YS(1)*XP KK = KK - 1 FN = FN - 1.0D0 240 CONTINUE 250 CONTINUE IF (ICASE.NE.2) RETURN DO 260 I=1,M3 YS(I) = YSS(I) 260 CONTINUE GO TO 90 270 CONTINUE IF (N.LT.NT) GO TO 290 C----------------------------------------------------------------------- C ICASE=1, NT.LE.N.LE.NL WITH FORWARD RECURSION C----------------------------------------------------------------------- 280 CONTINUE NN = N + M3 - 1 NFLG = MIN0(M-M3,1) ICASE = 1 GO TO 140 C----------------------------------------------------------------------- C ICASE=2, N.LT.NT.LT.NL WITH BOTH FORWARD AND BACKWARD RECURSION C----------------------------------------------------------------------- 290 CONTINUE NN = NT + 1 NFLG = MIN0(M-M3,1) ICASE = 2 GO TO 140 C----------------------------------------------------------------------- C X=0 CASE C----------------------------------------------------------------------- 300 CONTINUE FN = DBLE(FLOAT(N)) HN = 0.5D0*FN GR = DGAMRN(HN) Y(1) = HRTPI*GR IF (M.EQ.1) RETURN Y(2) = HRTPI/(HN*GR) IF (M.EQ.2) RETURN DO 310 K=3,M Y(K) = FN*Y(K-2)/(FN+1.0D0) FN = FN + 1.0D0 310 CONTINUE RETURN C----------------------------------------------------------------------- C UNDERFLOW ON KODE=1, X.GT.XLIM C----------------------------------------------------------------------- 320 CONTINUE IERR = 1 DO 330 I=1,M Y(I) = 0.0D0 330 CONTINUE RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 340 IF (NFLG.NE.0) GO TO 350 CALL XGETF(II) CALL XSETF(-1) 350 GO TO (360, 370, 380, 390, 400), ISET 360 CALL XERROR(34H DBSKIN, X IS NOT ZERO OR POSITIVE, 34, 2, 1) NFLG = 1 GO TO 10 370 CALL XERROR(34H DBSKIN, N IS NOT ZERO OR POSITIVE, 34, 2, 1) NFLG = 1 GO TO 20 380 CALL XERROR(27H DBSKIN, KODE IS NOT 1 OR 2, 27, 2, 1) NFLG = 1 GO TO 30 390 CALL XERROR(35H DBSKIN, M IS NOT GREATER THAN ZERO, 35, 2, 1) NFLG = 1 GO TO 40 400 CALL XERROR(48H DBSKIN, X AND N CANNOT BE ZERO AT THE SAME TIME, * 48, 2, 1) NFLG = 1 GO TO 50 410 CALL XSETF(II) CALL XERROR(36H DBSKIN, END INPUT ERRORS FOR DBSKIN, 36, 2, 1) RETURN END SUBROUTINE DBKISR(X, N, SUM) DBK 10 C C DBKISR COMPUTES REPEATED INTEGRALS OF THE K0 BESSEL FUNCTION C BY THE SERIES FOR N=0,1, AND 2. C INTEGER I, K, KK, KKN, K1, N, NP DOUBLE PRECISION AK, ATOL, BK, C, FK, FN, HX, HXS, POL, PR, SUM, * TKP, TOL, TRM, X, XLN DOUBLE PRECISION DPSIXN, D1MACH DIMENSION C(2) DATA C(1), C(2) /1.57079632679489662D+00,1.0D0/ C TOL = DMAX1(D1MACH(4),1.0D-18) IF (X.LT.TOL) GO TO 50 PR = 1.0D0 POL = 0.0D0 IF (N.EQ.0) GO TO 20 DO 10 I=1,N POL = -POL*X + C(I) PR = PR*X/DBLE(FLOAT(I)) 10 CONTINUE 20 CONTINUE HX = X*0.5D0 HXS = HX*HX XLN = DLOG(HX) NP = N + 1 TKP = 3.0D0 FK = 2.0D0 FN = DBLE(FLOAT(N)) BK = 4.0D0 AK = 2.0D0/((FN+1.0D0)*(FN+2.0D0)) SUM = AK*(DPSIXN(N+3)-DPSIXN(3)+DPSIXN(2)-XLN) ATOL = SUM*TOL*0.75D0 DO 30 K=2,20 AK = AK*(HXS/BK)*((TKP+1.0D0)/(TKP+FN+1.0D0))*(TKP/(TKP+FN)) K1 = K + 1 KK = K1 + K KKN = KK + N TRM = (DPSIXN(K1)+DPSIXN(KKN)-DPSIXN(KK)-XLN)*AK SUM = SUM + TRM IF (DABS(TRM).LE.ATOL) GO TO 40 TKP = TKP + 2.0D0 BK = BK + TKP FK = FK + 1.0D0 30 CONTINUE GO TO 80 40 CONTINUE SUM = (SUM*HXS+DPSIXN(NP)-XLN)*PR IF (N.EQ.1) SUM = -SUM SUM = POL + SUM RETURN C----------------------------------------------------------------------- C SMALL X CASE, X.LT.WORD TOLERANCE C----------------------------------------------------------------------- 50 CONTINUE IF (X.EQ.0.0D0 .AND. N.EQ.0) GO TO 70 IF (N.GT.0) GO TO 60 HX = X*0.5D0 SUM = DPSIXN(1) - DLOG(HX) RETURN 60 CONTINUE SUM = C(N) RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 70 CALL XERROR(44H DBKISR, X IS NOT GREATER THAN ZERO WHEN N=0, 44, * 2, 1) RETURN 80 CALL XERROR( * 51H DBKISR, X IS TOO LARGE FOR CONVERGENCE IN 20 TERMS, 51, 5, 1) RETURN END SUBROUTINE DBKIAS(X, N, KTRMS, T, ANS, IND, MS, GMRN, H) DBK 10 C C DBKIAS COMPUTES REPEATED INTEGRALS OF THE K0 BESSEL FUNCTION C BY THE ASYMPTOTIC EXPANSION C INTEGER I, II, IND, J, JMI, JN, K, KK, KM, KTRMS, MM, MP, MS, N DOUBLE PRECISION ANS, B, BND, DEN1, DEN2, DEN3, ER, ERR, FJ, FK, * FLN, FM1, GMRN, G1, GS, H, HN, HRTPI, RAT, RG1, RXP, RZ, RZX, S, * SS, SUMI, SUMJ, T, TOL, V, W, X, XP, Z DOUBLE PRECISION DGAMRN, D1MACH DIMENSION B(120), XP(16), S(31), H(1), V(52), W(52), T(50), * BND(15) C----------------------------------------------------------------------- C COEFFICIENTS OF POLYNOMIAL P(J-1,X), J=1,15 C----------------------------------------------------------------------- DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19), * B(20), B(21), B(22), B(23), B(24) /1.00000000000000000D+00, * 1.00000000000000000D+00,-2.00000000000000000D+00, * 1.00000000000000000D+00,-8.00000000000000000D+00, * 6.00000000000000000D+00,1.00000000000000000D+00, * -2.20000000000000000D+01,5.80000000000000000D+01, * -2.40000000000000000D+01,1.00000000000000000D+00, * -5.20000000000000000D+01,3.28000000000000000D+02, * -4.44000000000000000D+02,1.20000000000000000D+02, * 1.00000000000000000D+00,-1.14000000000000000D+02, * 1.45200000000000000D+03,-4.40000000000000000D+03, * 3.70800000000000000D+03,-7.20000000000000000D+02, * 1.00000000000000000D+00,-2.40000000000000000D+02, * 5.61000000000000000D+03/ DATA B(25), B(26), B(27), B(28), B(29), B(30), B(31), B(32), * B(33), B(34), B(35), B(36), B(37), B(38), B(39), B(40), B(41), * B(42), B(43), B(44), B(45), B(46), B(47), B(48) * /-3.21200000000000000D+04,5.81400000000000000D+04, * -3.39840000000000000D+04,5.04000000000000000D+03, * 1.00000000000000000D+00,-4.94000000000000000D+02, * 1.99500000000000000D+04,-1.95800000000000000D+05, * 6.44020000000000000D+05,-7.85304000000000000D+05, * 3.41136000000000000D+05,-4.03200000000000000D+04, * 1.00000000000000000D+00,-1.00400000000000000D+03, * 6.72600000000000000D+04,-1.06250000000000000D+06, * 5.76550000000000000D+06,-1.24400640000000000D+07, * 1.10262960000000000D+07,-3.73392000000000000D+06, * 3.62880000000000000D+05,1.00000000000000000D+00, * -2.02600000000000000D+03,2.18848000000000000D+05/ DATA B(49), B(50), B(51), B(52), B(53), B(54), B(55), B(56), * B(57), B(58), B(59), B(60), B(61), B(62), B(63), B(64), B(65), * B(66), B(67), B(68), B(69), B(70), B(71), B(72) * /-5.32616000000000000D+06,4.47650000000000000D+07, * -1.55357384000000000D+08,2.38904904000000000D+08, * -1.62186912000000000D+08,4.43390400000000000D+07, * -3.62880000000000000D+06,1.00000000000000000D+00, * -4.07200000000000000D+03,6.95038000000000000D+05, * -2.52439040000000000D+07,3.14369720000000000D+08, * -1.64838430400000000D+09,4.00269508800000000D+09, * -4.64216395200000000D+09,2.50748121600000000D+09, * -5.68356480000000000D+08,3.99168000000000000D+07, * 1.00000000000000000D+00,-8.16600000000000000D+03, * 2.17062600000000000D+06,-1.14876376000000000D+08, * 2.05148277600000000D+09,-1.55489607840000000D+10/ DATA B(73), B(74), B(75), B(76), B(77), B(78), B(79), B(80), * B(81), B(82), B(83), B(84), B(85), B(86), B(87), B(88), B(89), * B(90), B(91), B(92), B(93), B(94), B(95), B(96) * /5.60413987840000000D+10,-1.01180433024000000D+11, * 9.21997902240000000D+10,-4.07883018240000000D+10, * 7.82771904000000000D+09,-4.79001600000000000D+08, * 1.00000000000000000D+00,-1.63560000000000000D+04, * 6.69969600000000000D+06,-5.07259276000000000D+08, * 1.26698177760000000D+10,-1.34323420224000000D+11, * 6.87720046384000000D+11,-1.81818864230400000D+12, * 2.54986547342400000D+12,-1.88307966182400000D+12, * 6.97929436800000000D+11,-1.15336085760000000D+11, * 6.22702080000000000D+09,1.00000000000000000D+00, * -3.27380000000000000D+04,2.05079880000000000D+07, * -2.18982980800000000D+09,7.50160522280000000D+10/ DATA B(97), B(98), B(99), B(100), B(101), B(102), B(103), B(104), * B(105), B(106), B(107), B(108), B(109), B(110), B(111), B(112), * B(113), B(114), B(115), B(116), B(117), B(118) * /-1.08467651241600000D+12,7.63483214939200000D+12, * -2.82999100661120000D+13,5.74943734645920000D+13, * -6.47283751398720000D+13,3.96895780558080000D+13, * -1.25509040179200000D+13,1.81099255680000000D+12, * -8.71782912000000000D+10,1.00000000000000000D+00, * -6.55040000000000000D+04,6.24078900000000000D+07, * -9.29252692000000000D+09,4.29826006340000000D+11, * -8.30844432796800000D+12,7.83913848313120000D+13, * -3.94365587815520000D+14,1.11174747256968000D+15, * -1.79717122069056000D+15,1.66642448627145600D+15, * -8.65023253219584000D+14,2.36908271543040000D+14/ DATA B(119), B(120) /-3.01963769856000000D+13, * 1.30767436800000000D+12/ C----------------------------------------------------------------------- C BOUNDS B(M,K) , K=M-3 C----------------------------------------------------------------------- DATA BND(1), BND(2), BND(3), BND(4), BND(5), BND(6), BND(7), * BND(8), BND(9), BND(10), BND(11), BND(12), BND(13), BND(14), * BND(15) /1.0D0,1.0D0,1.0D0,1.0D0,3.10D0,5.18D0,11.7D0,29.8D0, * 90.4D0,297.0D0,1070.0D0,4290.0D0,18100.0D0,84700.0D0,408000.0D0/ DATA HRTPI /8.86226925452758014D-01/ C TOL = DMAX1(D1MACH(4),1.0D-18) FLN = DBLE(FLOAT(N)) RZ = 1.0D0/(X+FLN) RZX = X*RZ Z = 0.5D0*(X+FLN) IF (IND.GT.1) GO TO 10 GMRN = DGAMRN(Z) 10 CONTINUE GS = HRTPI*GMRN G1 = GS + GS RG1 = 1.0D0/G1 GMRN = (RZ+RZ)/GMRN IF (IND.GT.1) GO TO 70 C----------------------------------------------------------------------- C EVALUATE ERROR FOR M=MS C----------------------------------------------------------------------- HN = 0.5D0*FLN DEN2 = KTRMS + KTRMS + N DEN3 = DEN2 - 2.0D0 DEN1 = X + DEN2 ERR = RG1*(X+X)/(DEN1-1.0D0) IF (N.EQ.0) GO TO 20 RAT = 1.0D0/(FLN*FLN) 20 CONTINUE IF (KTRMS.EQ.0) GO TO 30 FJ = DBLE(FLOAT(KTRMS)) RAT = 0.25D0/(HRTPI*DEN3*DSQRT(FJ)) 30 CONTINUE ERR = ERR*RAT FJ = -3.0D0 DO 50 J=1,15 IF (J.LE.5) ERR = ERR/DEN1 FM1 = DMAX1(1.0D0,FJ) FJ = FJ + 1.0D0 ER = BND(J)*ERR IF (KTRMS.EQ.0) GO TO 40 ER = ER/FM1 IF (ER.LT.TOL) GO TO 60 IF (J.GE.5) ERR = ERR/DEN3 GO TO 50 40 CONTINUE ER = ER*(1.0D0+HN/FM1) IF (ER.LT.TOL) GO TO 60 IF (J.GE.5) ERR = ERR/FLN 50 CONTINUE GO TO 200 60 CONTINUE MS = J 70 CONTINUE MM = MS + MS MP = MM + 1 C----------------------------------------------------------------------- C H(K)=(-Z)**(K)*(PSI(K-1,Z)-PSI(K-1,Z+0.5))/GAMMA(K) , K=1,2,...,MM C----------------------------------------------------------------------- IF (IND.GT.1) GO TO 80 CALL DHKSEQ(Z, MM, H) GO TO 100 80 CONTINUE RAT = Z/(Z-0.5D0) RXP = RAT DO 90 I=1,MM H(I) = RXP*(1.0D0-H(I)) RXP = RXP*RAT 90 CONTINUE 100 CONTINUE C----------------------------------------------------------------------- C SCALED S SEQUENCE C----------------------------------------------------------------------- S(1) = 1.0D0 FK = 1.0D0 DO 120 K=2,MP SS = 0.0D0 KM = K - 1 I = KM DO 110 J=1,KM SS = SS + S(J)*H(I) I = I - 1 110 CONTINUE S(K) = SS/FK FK = FK + 1.0D0 120 CONTINUE C----------------------------------------------------------------------- C SCALED S-TILDA SEQUENCE C----------------------------------------------------------------------- IF (KTRMS.EQ.0) GO TO 160 FK = 0.0D0 SS = 0.0D0 RG1 = RG1/Z DO 130 K=1,KTRMS V(K) = Z/(Z+FK) W(K) = T(K)*V(K) SS = SS + W(K) FK = FK + 1.0D0 130 CONTINUE S(1) = S(1) - SS*RG1 DO 150 I=2,MP SS = 0.0D0 DO 140 K=1,KTRMS W(K) = W(K)*V(K) SS = SS + W(K) 140 CONTINUE S(I) = S(I) - SS*RG1 150 CONTINUE 160 CONTINUE C----------------------------------------------------------------------- C SUM ON J C----------------------------------------------------------------------- SUMJ = 0.0D0 JN = 1 RXP = 1.0D0 XP(1) = 1.0D0 DO 190 J=1,MS JN = JN + J - 1 XP(J+1) = XP(J)*RZX RXP = RXP*RZ C----------------------------------------------------------------------- C SUM ON I C----------------------------------------------------------------------- SUMI = 0.0D0 II = JN DO 180 I=1,J JMI = J - I + 1 KK = J + I + 1 DO 170 K=1,JMI V(K) = S(KK)*XP(K) KK = KK + 1 170 CONTINUE CALL DBDIFF(JMI, V) SUMI = SUMI + B(II)*V(JMI)*XP(I+1) II = II + 1 180 CONTINUE SUMJ = SUMJ + SUMI*RXP 190 CONTINUE ANS = GS*(S(1)-SUMJ) RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 200 CALL XERROR(69H DBKIAS, ERROR TEST NOT MET. INCREASE TOL IN DBKIAS * OR NLIM IN DBSKIN, 69, 5, 1) RETURN END SUBROUTINE DBDIFF(L, V) DBD 10 C C DBDIFF COMPUTES THE SUM OF B(L,K)*V(K)*(-1)**K WHERE B(L,K) C ARE THE BINOMIAL COEFFICIENTS. TRUNCATED SUMS ARE COMPUTED BY C SETTING LAST PART OF THE V VECTOR TO ZERO. ON RETURN, THE BINOMIAL C SUM IS IN V(L). C INTEGER I, J, K, L DOUBLE PRECISION V DIMENSION V(1) IF (L.EQ.1) RETURN DO 20 J=2,L K = L DO 10 I=J,L V(K) = V(K-1) - V(K) K = K - 1 10 CONTINUE 20 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DGAMRN(X) DGA 10 C C WRITTEN BY D.E. AMOS, JANUARY, 1981. C C REFERENCES C THE SPECIAL FUNCTIONS AND THEIR APPROXIMATIONS, VOL. 1 BY C Y.L. LUKE, MATH IN SCI. AND ENG. SERIES 53, ACADEMIC PRESS, C NEW YORK, 1969, PP 34-35. C C ABSTRACT *** A DOUBLE PRECISION ROUTINE*** C DGAMRN COMPUTES THE GAMMA FUNCTION RATIO GAMMA(X)/GAMMA(X+0.5) C FOR REAL X.GT.0. IF X.GE.XMIN, AN ASYMPTOTIC EXPANSION IS C EVALUATED. IF X.LT.XMIN, AN INTEGER IS ADDED TO X TO FORM A C NEW VALUE OF X.GE.XMIN AND THE ASYMPTOTIC EXPANSION IS EVAL- C UATED FOR THIS NEW VALUE OF X. SUCCESSIVE APPLICATION OF THE C RECURRENCE RELATION C C W(X)=W(X+1)*(1+0.5/X) C C REDUCES THE ARGUMENT TO ITS ORIGINAL VALUE. XMIN AND COMP- C UTATIONAL TOLERANCES ARE COMPUTED AS A FUNCTION OF THE NUMBER C OF DIGITS CARRIED IN A WORD BY CALLS TO I1MACH AND D1MACH. C HOWEVER, THE COMPUTATIONAL ACCURACY IS LIMITED TO THE MAX- C IMUM OF UNIT ROUNDOFF (=D1MACH(4)) AND 1.0D-18 SINCE CRITICAL C CONSTANTS ARE GIVEN TO ONLY 18 DIGITS. C C DGAMRN CALLS D1MACH, I1MACH, XERROR. C C DESCRIPTION OF ARGUMENTS C C INPUT X IS DOUBLE PRECISION C X - ARGUMENT, X.GT.0.0D0 C C OUTPUT DGAMRN IS DOUBLE PRECISION C DGAMRN - RATIO GAMMA(X)/GAMMA(X+0.5) C C ERROR CONDITIONS C X.LE.0 IS A FATAL ERROR C INTEGER I, I1M11, K, MX, NX INTEGER I1MACH DOUBLE PRECISION FLN, GR, RLN, S, TOL, TRM, X, XDMY, XINC, XM, * XMIN, XP, XSQ DOUBLE PRECISION D1MACH DIMENSION GR(12) C DATA GR(1), GR(2), GR(3), GR(4), GR(5), GR(6), GR(7), GR(8), * GR(9), GR(10), GR(11), GR(12) /1.00000000000000000D+00, * -1.56250000000000000D-02,2.56347656250000000D-03, * -1.27983093261718750D-03,1.34351104497909546D-03, * -2.43289663922041655D-03,6.75423753364157164D-03, * -2.66369606131178216D-02,1.41527455519564332D-01, * -9.74384543032201613D-01,8.43686251229783675D+00, * -8.97258321640552515D+01/ C IF (X.LE.0.0D0) GO TO 60 NX = INT(SNGL(X)) TOL = DMAX1(D1MACH(4),1.0D-18) I1M11 = I1MACH(14) RLN = D1MACH(5)*DBLE(FLOAT(I1M11)) FLN = DMIN1(RLN,20.0D0) FLN = DMAX1(FLN,3.0D0) FLN = FLN - 3.0D0 XM = 2.0D0 + FLN*(0.2366D0+0.01723D0*FLN) MX = INT(SNGL(XM)) + 1 XMIN = DBLE(FLOAT(MX)) XDMY = X - 0.25D0 XINC = 0.0D0 IF (X.GE.XMIN) GO TO 10 XINC = XMIN - DBLE(FLOAT(NX)) XDMY = XDMY + XINC 10 CONTINUE S = 1.0D0 IF (XDMY*TOL.GT.1.0D0) GO TO 30 XSQ = 1.0D0/(XDMY*XDMY) XP = XSQ DO 20 K=2,12 TRM = GR(K)*XP IF (DABS(TRM).LT.TOL) GO TO 30 S = S + TRM XP = XP*XSQ 20 CONTINUE 30 CONTINUE S = S/DSQRT(XDMY) IF (XINC.NE.0.0D0) GO TO 40 DGAMRN = S RETURN 40 CONTINUE NX = INT(SNGL(XINC)) XP = 0.0D0 DO 50 I=1,NX S = S*(1.0D0+0.5D0/(X+XP)) XP = XP + 1.0D0 50 CONTINUE DGAMRN = S RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 60 CALL XERROR(35H DGAMRN, X IS NOT GREATER THAN XERO, 35, 2, 1) RETURN END SUBROUTINE DHKSEQ(X, M, H) DHK 10 C C DHKSEQ IS AN ADAPTATION OF SUBROUTINE PSIFN, ACM TRANS. MATH. C SOFTWARE, 1983 TO GENERATE H(K,X)=(-X)**(K+1)*(PSI(K,X)- C -PSI(K,X+0.5))/GAMMA(K+1) FOR K=0,...,M. C INTEGER I, J, K, M, MX, NX INTEGER I1MACH DOUBLE PRECISION B, FK, FLN, FN, FNP, H, HRX, RLN, RXSQ, R1M5, S, * SLOPE, T, TK, TRM, TRMH, TRMR, TST, U, V, WDTOL, X, XDMY, XH, * XINC, XM, XMIN, YINT DOUBLE PRECISION D1MACH DIMENSION B(22), TRM(22), TRMR(25), TRMH(25), U(25), V(25), H(1) C----------------------------------------------------------------------- C SCALED BERNOULLI NUMBERS 2.0*B(2K)*(1-2**(-2K)) C----------------------------------------------------------------------- DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19), * B(20), B(21), B(22) /1.00000000000000000D+00, * -5.00000000000000000D-01,2.50000000000000000D-01, * -6.25000000000000000D-02,4.68750000000000000D-02, * -6.64062500000000000D-02,1.51367187500000000D-01, * -5.06103515625000000D-01,2.33319091796875000D+00, * -1.41840972900390625D+01,1.09941936492919922D+02, * -1.05824747562408447D+03,1.23842434241771698D+04, * -1.73160495905935764D+05,2.85103429084961116D+06, * -5.45964619322445132D+07,1.20316174668075304D+09, * -3.02326315271452307D+10,8.59229286072319606D+11, * -2.74233104097776039D+13,9.76664637943633248D+14, * -3.85931586838450360D+16/ C C WDTOL = DMAX1(D1MACH(4),1.0D-18) FN = DBLE(FLOAT(M-1)) FNP = FN + 1.0D0 C----------------------------------------------------------------------- C COMPUTE XMIN C----------------------------------------------------------------------- R1M5 = D1MACH(5) RLN = R1M5*DBLE(FLOAT(I1MACH(14))) RLN = DMIN1(RLN,18.06D0) FLN = DMAX1(RLN,3.0D0) - 3.0D0 YINT = 3.50D0 + 0.40D0*FLN SLOPE = 0.21D0 + FLN*(0.0006038D0*FLN+0.008677D0) XM = YINT + SLOPE*FN MX = INT(SNGL(XM)) + 1 XMIN = DBLE(FLOAT(MX)) C----------------------------------------------------------------------- C GENERATE H(M-1,XDMY)*XDMY**(M) BY THE ASYMPTOTIC EXPANSION C----------------------------------------------------------------------- XDMY = X XINC = 0.0D0 IF (X.GE.XMIN) GO TO 10 NX = INT(SNGL(X)) XINC = XMIN - DBLE(FLOAT(NX)) XDMY = X + XINC 10 CONTINUE RXSQ = 1.0D0/(XDMY*XDMY) HRX = 0.5D0/XDMY TST = 0.5D0*WDTOL T = FNP*HRX C----------------------------------------------------------------------- C INITIALIZE COEFFICIENT ARRAY C----------------------------------------------------------------------- S = T*B(3) IF (DABS(S).LT.TST) GO TO 30 TK = 2.0D0 DO 20 K=4,22 T = T*((TK+FN+1.0D0)/(TK+1.0D0))*((TK+FN)/(TK+2.0D0))*RXSQ TRM(K) = T*B(K) IF (DABS(TRM(K)).LT.TST) GO TO 30 S = S + TRM(K) TK = TK + 2.0D0 20 CONTINUE GO TO 110 30 CONTINUE H(M) = S + 0.5D0 IF (M.EQ.1) GO TO 70 C----------------------------------------------------------------------- C GENERATE LOWER DERIVATIVES, I.LT.M-1 C----------------------------------------------------------------------- DO 60 I=2,M FNP = FN FN = FN - 1.0D0 S = FNP*HRX*B(3) IF (DABS(S).LT.TST) GO TO 50 FK = FNP + 3.0D0 DO 40 K=4,22 TRM(K) = TRM(K)*FNP/FK IF (DABS(TRM(K)).LT.TST) GO TO 50 S = S + TRM(K) FK = FK + 2.0D0 40 CONTINUE GO TO 110 50 CONTINUE MX = M - I + 1 H(MX) = S + 0.5D0 60 CONTINUE 70 CONTINUE IF (XINC.EQ.0.0D0) RETURN C----------------------------------------------------------------------- C RECUR BACKWARD FROM XDMY TO X C----------------------------------------------------------------------- XH = X + 0.5D0 S = 0.0D0 NX = INT(SNGL(XINC)) DO 80 I=1,NX TRMR(I) = X/(X+DBLE(FLOAT(NX-I))) U(I) = TRMR(I) TRMH(I) = X/(XH+DBLE(FLOAT(NX-I))) V(I) = TRMH(I) S = S + U(I) - V(I) 80 CONTINUE MX = NX + 1 TRMR(MX) = X/XDMY U(MX) = TRMR(MX) H(1) = H(1)*TRMR(MX) + S IF (M.EQ.1) RETURN DO 100 J=2,M S = 0.0D0 DO 90 I=1,NX TRMR(I) = TRMR(I)*U(I) TRMH(I) = TRMH(I)*V(I) S = S + TRMR(I) - TRMH(I) 90 CONTINUE TRMR(MX) = TRMR(MX)*U(MX) H(J) = H(J)*TRMR(MX) + S 100 CONTINUE RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 110 CALL XERROR(45H DHKSEQ, ERROR TEST IN ASYMPTOTIC SUM NOT MET, 45, * 2, 1) RETURN END SUBROUTINE DEXINT(X, N, KODE, M, TOL, EN, IERR) DEX 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, PP. 365-377, 420-428. C C ABSTRACT *** A DOUBLE PRECISION ROUTINE *** C DEXINT 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) = DEXP(-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) = DEXP(-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 THE MAXIMUM NUMBER OF SIGNIFICANT DIGITS OBTAINABLE C IS THE SMALLER OF 18 AND THE NUMBER OF DIGITS CARRIED IN C DOUBLE PRECISION ARITHMETIC. C C DEXINT CALLS I1MACH, D1MACH, DPSIXN, XERROR C C DESCRIPTION OF ARGUMENTS C C INPUT * X AND TOL ARE DOUBLE PRECISION * 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 DEXP(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 IS THE LARGER OF DOUBLE PRECISION UNIT C ROUNDOFF = D1MACH(4) AND 0.5D-18 C C OUTPUT * EN IS A DOUBLE PRECISION VECTOR * C EN A VECTOR OF DIMENSION AT LEAST M CONTAINING VALUES C EN(K) = E(N+K-1,X) OR DEXP(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.0D0 , K=1,M RETURNED ON KODE=1 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***END PROLOG C DOUBLE PRECISION A, AA, AAMS, AH, AK, AT, B, BK, BT, CC, CNORM, * CT, EM, EMX, EN, ETOL, FNM, FX, PT, P1, P2, S, TOL, TX, X, XCUT, * XLIM, XTOL, Y, YT, Y1, Y2 DOUBLE PRECISION D1MACH, DPSIXN INTEGER I, IC, ICASE, ICT, IERR, IK, IND, IX, I1M, JSET, K, KK, * KN, KODE, KS, M, ML, MU, N, ND, NM INTEGER I1MACH DIMENSION EN(1), A(99), B(99), Y(2) C DATA XCUT / 2.0D0 / C ETOL = DMAX1(D1MACH(4),0.5D-18) I1M = -I1MACH(15) PT = 2.3026D0*DBLE(FLOAT(I1M))*D1MACH(5) XLIM = PT - 6.907755D0 BT = PT + DBLE(FLOAT(N+M-1)) IF (BT.GT.1000.0D0) XLIM = PT - DLOG(BT) 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.1D0) GO TO 290 C IERR = 0 IF (X.GT.XCUT) GO TO 100 IF (X.LT.0.0D0) GO TO 300 IF (X.EQ.0.0D0 .AND. N.EQ.1) GO TO 310 IF (X.EQ.0.0D0 .AND. N.GT.1) GO TO 80 C----------------------------------------------------------------------- C SERIES FOR E(N,X) FOR X.LE.XCUT C----------------------------------------------------------------------- TX = X + 0.5D0 IX = INT(SNGL(TX)) C----------------------------------------------------------------------- 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 C----------------------------------------------------------------------- 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 = DBLE(FLOAT(NM)) S = 0.0D0 XTOL = 3.0D0*TOL IF (ND.EQ.1) GO TO 10 XTOL = 0.3333D0*TOL S = 1.0D0/FNM 10 CONTINUE AA = 1.0D0 AK = 1.0D0 IC = 35 IF (X.LT.ETOL) IC = 1 DO 50 I=1,IC AA = -AA*X/AK IF (I.EQ.NM) GO TO 30 S = S - AA/(AK-FNM) IF (DABS(AA).LE.XTOL*DABS(S)) GO TO 20 AK = AK + 1.0D0 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.0D0 GO TO 50 30 S = S + AA*(-DLOG(X)+DPSIXN(ND)) XTOL = 3.0D0*TOL 40 AK = AK + 1.0D0 50 CONTINUE IF (IC.NE.1) GO TO 320 60 IF (ND.EQ.1) S = S + (-DLOG(X)+DPSIXN(1)) IF (KODE.EQ.2) S = S*DEXP(X) EN(1) = S EMX = 1.0D0 IF (M.EQ.1) GO TO 70 EN(IND) = S AA = DBLE(FLOAT(KS)) IF (KODE.EQ.1) EMX = DEXP(-X) GO TO (220, 240), ICASE 70 IF (ICASE.EQ.2) RETURN IF (KODE.EQ.1) EMX = DEXP(-X) EN(1) = (EMX-S)/X RETURN 80 CONTINUE DO 90 I=1,M EN(I) = 1.0D0/DBLE(FLOAT(N+I-2)) 90 CONTINUE RETURN C----------------------------------------------------------------------- C BACKWARD RECURSIVE MILLER ALGORITHM FOR C E(N,X)=DEXP(-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.0D0 IF (KODE.EQ.2) GO TO 130 IF (X.LE.XLIM) GO TO 120 IERR = 1 DO 110 I=1,M EN(I) = 0.0D0 110 CONTINUE RETURN 120 EMX = DEXP(-X) 130 CONTINUE TX = X + 0.5D0 IX = INT(SNGL(TX)) 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 = DBLE(FLOAT(IK)) JSET = 1 + KS - (IK+IK) C----------------------------------------------------------------------- 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. C----------------------------------------------------------------------- IC = 0 AA = AH + AH AAMS = AA - 1.0D0 AAMS = AAMS*AAMS TX = X + X FX = TX + TX AK = AH XTOL = TOL IF (TOL.LE.1.0D-3) XTOL = 20.0D0*TOL CT = AAMS + FX*AH EM = (AH+1.0D0)/((X+AA)*XTOL*DSQRT(CT)) BK = AA CC = AH*AH C----------------------------------------------------------------------- C FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD C RECURSION C----------------------------------------------------------------------- P1 = 0.0D0 P2 = 1.0D0 190 CONTINUE IF (IC.EQ.99) GO TO 330 IC = IC + 1 AK = AK + 1.0D0 AT = BK/(BK+AK+CC+DBLE(FLOAT(IC))) BK = BK + AK + AK A(IC) = AT BT = (AK+AK+X)/(AK+1.0D0) B(IC) = BT PT = P2 P2 = BT*P2 - AT*P1 P1 = PT CT = CT + FX EM = EM*AT*(1.0D0-TX/CT) IF (EM*(AK+1.0D0).GT.P1*P1) GO TO 190 ICT = IC KK = IC + 1 BT = TX/(CT+FX) Y2 = (BK/(BK+CC+DBLE(FLOAT(KK))))*(P1/P2)*(1.0D0-BT+0.375D0*BT*BT) Y1 = 1.0D0 C----------------------------------------------------------------------- C BACKWARD RECURRENCE FOR C Y1= C*U( A ,A,X) C Y2= C*(A/(1+A/2))*U(A+1,A,X) C----------------------------------------------------------------------- DO 200 K=1,ICT KK = KK - 1 YT = Y1 Y1 = (B(KK)*Y1-Y2)/A(KK) Y2 = YT 200 CONTINUE C----------------------------------------------------------------------- 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 C----------------------------------------------------------------------- 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 = DBLE(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.0D0 EN(K) = (EMX-AA*EN(K+1))/X K = K - 1 230 CONTINUE IF (MU.LE.0) RETURN AA = DBLE(FLOAT(KS)) 240 K = IND DO 250 I=1,MU EN(K+1) = (EMX-X*EN(K))/AA AA = AA + 1.0D0 K = K + 1 250 CONTINUE RETURN C----------------------------------------------------------------------- C ERROR MESSAGES C----------------------------------------------------------------------- 260 CALL XERROR(31HIN DEXINT, N NOT GREATER THAN 0, 31, 2, 1) RETURN 270 CALL XERROR(26HIN DEXINT, KODE NOT 1 OR 2, 26, 2, 1) RETURN 280 CALL XERROR(31HIN DEXINT, M NOT GREATER THAN 0, 31, 2, 1) RETURN 290 CALL XERROR(32HIN DEXINT, TOL NOT WITHIN LIMITS, 32, 2, 1) RETURN 300 CALL XERROR(36HIN DEXINT, X IS NOT ZERO OR POSITIVE, 36, 2, 1) RETURN 310 CALL XERROR(66HIN DEXINT, THE EXPONENTIAL INTEGRAL IS NOT DEFINED *FOR X=0 AND N=1, 66, 2, 1) RETURN 320 CALL XERROR(73HIN DEXINT, RELATIVE ERROR TEST FOR SERIES TERMINATI *ON NOT MET IN 36 TERMS, 73, 3, 1) RETURN 330 CALL XERROR(68HIN DEXINT, TERMINATION TEST FOR MILLER ALGORITHM NO *T MET IN 99 STEPS, 68, 3, 1) RETURN 340 CALL XERROR(92HIN DEXINT, AN ERROR IN PLACING INT(X+0.5) WITH RESP *ECT TO N AND N+M-1 OCCURRED FOR X.GT.XCUT, 92, 8, 1) RETURN END DOUBLE PRECISION FUNCTION DPSIXN(N) DPS 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 C INTEGER N, K DOUBLE PRECISION AX, B, C, FN, RFN2, TRM, S, WDTOL DOUBLE PRECISION D1MACH DIMENSION B(6), C(100) C C DPSIXN(N), N = 1,100 DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), 2 C(19), C(20), C(21), C(22), C(23), C(24)/ 3 -5.77215664901532861D-01, 4.22784335098467139D-01, 4 9.22784335098467139D-01, 1.25611766843180047D+00, 5 1.50611766843180047D+00, 1.70611766843180047D+00, 6 1.87278433509846714D+00, 2.01564147795561000D+00, 7 2.14064147795561000D+00, 2.25175258906672111D+00, 8 2.35175258906672111D+00, 2.44266167997581202D+00, 9 2.52599501330914535D+00, 2.60291809023222227D+00, A 2.67434666166079370D+00, 2.74101332832746037D+00, B 2.80351332832746037D+00, 2.86233685773922507D+00, C 2.91789241329478063D+00, 2.97052399224214905D+00, D 3.02052399224214905D+00, 3.06814303986119667D+00, E 3.11359758531574212D+00, 3.15707584618530734D+00/ DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ 3 3.19874251285197401D+00, 3.23874251285197401D+00, 4 3.27720405131351247D+00, 3.31424108835054951D+00, 5 3.34995537406483522D+00, 3.38443813268552488D+00, 6 3.41777146601885821D+00, 3.45002953053498724D+00, 7 3.48127953053498724D+00, 3.51158256083801755D+00, 8 3.54099432554389990D+00, 3.56956575411532847D+00, 9 3.59734353189310625D+00, 3.62437055892013327D+00, A 3.65068634839381748D+00, 3.67632737403484313D+00, B 3.70132737403484313D+00, 3.72571761793728215D+00, C 3.74952714174680596D+00, 3.77278295570029433D+00, D 3.79551022842756706D+00, 3.81773245064978928D+00, E 3.83947158108457189D+00, 3.86074817682925274D+00/ DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), 2 C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ 3 3.88158151016258607D+00, 3.90198967342789220D+00, 4 3.92198967342789220D+00, 3.94159751656514710D+00, 5 3.96082828579591633D+00, 3.97969621032421822D+00, 6 3.99821472884273674D+00, 4.01639654702455492D+00, 7 4.03425368988169777D+00, 4.05179754953082058D+00, 8 4.06903892884116541D+00, 4.08598808138353829D+00, 9 4.10265474805020496D+00, 4.11904819067315578D+00, A 4.13517722293122029D+00, 4.15105023880423617D+00, B 4.16667523880423617D+00, 4.18205985418885155D+00, C 4.19721136934036670D+00, 4.21213674247469506D+00, D 4.22684262482763624D+00, 4.24133537845082464D+00, E 4.25562109273653893D+00, 4.26970559977879245D+00/ DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ 3 4.28359448866768134D+00, 4.29729311880466764D+00, 4 4.31080663231818115D+00, 4.32413996565151449D+00, 5 4.33729786038835659D+00, 4.35028487337536958D+00, 6 4.36310538619588240D+00, 4.37576361404398366D+00, 7 4.38826361404398366D+00, 4.40060929305632934D+00, 8 4.41280441500754886D+00, 4.42485260777863319D+00, 9 4.43675736968339510D+00, 4.44852207556574804D+00, A 4.46014998254249223D+00, 4.47164423541605544D+00, B 4.48300787177969181D+00, 4.49424382683587158D+00, C 4.50535493794698269D+00, 4.51634394893599368D+00, D 4.52721351415338499D+00, 4.53796620232542800D+00, E 4.54860450019776842D+00, 4.55913081598724211D+00/ DATA C(97), C(98), C(99), C(100)/ 1 4.56954748265390877D+00, 4.57985676100442424D+00, 2 4.59006084263707730D+00, 4.60016185273808740D+00/ C COEFFICIENTS OF ASYMPTOTIC EXPANSION DATA B(1), B(2), B(3), B(4), B(5), B(6)/ 1 8.33333333333333333D-02, -8.33333333333333333D-03, 2 3.96825396825396825D-03, -4.16666666666666666D-03, 3 7.57575757575757576D-03, -2.10927960927960928D-02/ C IF (N.GT.100) GO TO 10 DPSIXN = C(N) RETURN 10 CONTINUE WDTOL = DMAX1(D1MACH(4),1.0D-18) FN = DBLE(FLOAT(N)) AX = 1.0D0 S = -0.5D0/FN IF (DABS(S).LE.WDTOL) GO TO 30 RFN2 = 1.0D0/(FN*FN) DO 20 K=1,6 AX = AX*RFN2 TRM = -B(K)*AX IF (DABS(TRM).LT.WDTOL) GO TO 30 S = S + TRM 20 CONTINUE 30 CONTINUE DPSIXN = S + DLOG(FN) RETURN END C=========================== subsidiary routines C=========================== used to be in the file 609.slatec SUBROUTINE FDUMP FDUMP 2 C***BEGIN PROLOGUE FDUMP C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q2,N2 C***KEYWORDS SYMBOLIC DUMP,DUMP C***DATE WRITTEN 8/79 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C SYMBOLIC DUMP (SHOULD BE LOCALLY WRITTEN) C***DESCRIPTION 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 LATEST REVISION --- 23 MAY 1979 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE FDUMP C C***FIRST EXECUTABLE STATEMENT FDUMP RETURN END FUNCTION J4SAVE(IWHICH,IVALUE,ISET) J4SAVE 2 C***BEGIN PROLOGUE J4SAVE C***REFER TO XERROR 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 WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C ADAPTED FROM BELL LABORATORIES PORT LIBRARY ERROR HANDLER C LATEST REVISION --- 23 MAY 1979 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE J4SAVE C LOGICAL ISET INTEGER IPARAM(9) 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 SUBROUTINE S88FMT(N,IVALUE,IFMT) S88FMT 2 C***BEGIN PROLOGUE S88FMT C***REFER TO XERROR C ABSTRACT C S88FMT REPLACES IFMT(1), ... ,IFMT(N) WITH THE C CHARACTERS CORRESPONDING TO THE N LEAST SIGNIFICANT C DIGITS OF IVALUE. C C TAKEN FROM THE BELL LABORATORIES PORT LIBRARY ERROR HANDLER C LATEST REVISION --- 7 JUNE 1978 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE S88FMT C DIMENSION IFMT(N),IDIGIT(10) DATA IDIGIT(1),IDIGIT(2),IDIGIT(3),IDIGIT(4),IDIGIT(5), 1 IDIGIT(6),IDIGIT(7),IDIGIT(8),IDIGIT(9),IDIGIT(10) 2 /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ C***FIRST EXECUTABLE STATEMENT S88FMT NT = N IT = IVALUE 10 IF (NT .EQ. 0) RETURN INDEX = MOD(IT,10) IFMT(NT) = IDIGIT(INDEX+1) IT = IT/10 NT = NT - 1 GO TO 10 END SUBROUTINE XERABT(MESSG,NMESSG) XERABT 2 C***BEGIN PROLOGUE XERABT C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,ABORT C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C ABORTS PROGRAM EXECUTION AND PRINTS ERROR MESSAGE C***DESCRIPTION C ABSTRACT C ***NOTE*** MACHINE DEPENDENT ROUTINE C XERABT 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 AND NMESSG ARE AS IN XERROR, EXCEPT THAT NMESSG MAY C BE ZERO, IN WHICH CASE NO MESSAGE IS BEING SUPPLIED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 7 JUNE 1978 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERABT C DIMENSION MESSG(NMESSG) C***FIRST EXECUTABLE STATEMENT XERABT STOP C RETURN END SUBROUTINE XERCLR XERCLR 2 C***BEGIN PROLOGUE XERCLR C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,CLEAR C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C RESETS CURRENT ERROR NUMBER TO ZERO C***DESCRIPTION C ABSTRACT C THIS ROUTINE SIMPLY RESETS THE CURRENT ERROR NUMBER TO ZERO. C THIS MAY BE NECESSARY TO DO IN ORDER TO DETERMINE THAT C A CERTAIN ERROR HAS OCCURRED AGAIN SINCE THE LAST TIME C NUMXER WAS REFERENCED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 7 JUNE 1978 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XERCLR C C***FIRST EXECUTABLE STATEMENT XERCLR JUNK = J4SAVE(1,0,.TRUE.) RETURN END SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL) XERCTL 2 C***BEGIN PROLOGUE XERCTL C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,CONTROL C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C ALLOWS USER CONTROL OVER HANDLING OF INDIVIDUAL ERRORS C***DESCRIPTION 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 XERCTL. C IF THE USER HAS PROVIDED HIS OWN VERSION OF XERCTL, 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 MESSG1 - THE FIRST WORD (ONLY) OF THE ERROR MESSAGE. C NMESSG - SAME AS IN THE CALL TO XERROR OR XERRWV. C NERR - SAME AS IN THE CALL TO XERROR OR XERRWV. C LEVEL - SAME AS IN THE CALL TO XERROR OR XERRWV. 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 C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERCTL C C***FIRST EXECUTABLE STATEMENT XERCTL RETURN END SUBROUTINE XERDMP XERDMP 2 C***BEGIN PROLOGUE XERDMP C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,TABLES C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C PRINTS THE ERROR TABLES AND THEN CLEARS THEM C***DESCRIPTION C ABSTRACT C XERDMP PRINTS THE ERROR TABLES, THEN CLEARS THEM. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 7 JUNE 1978 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED XERSAV C***END PROLOGUE XERDMP C C***FIRST EXECUTABLE STATEMENT XERDMP CALL XERSAV(1H ,0,0,0,KOUNT) RETURN END SUBROUTINE XERMAX(MAX) XERMAX 2 C***BEGIN PROLOGUE XERMAX C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,MAXIMUM ERROR COUNT C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C SETS MAXIMUM NUMBER OF TIMES ANY ERROR MESSAGE IS TO BE PRINTED C***DESCRIPTION C ABSTRACT C XERMAX SETS THE MAXIMUM NUMBER OF TIMES ANY MESSAGE C IS TO BE PRINTED. THAT IS, NON-FATAL MESSAGES ARE C NOT TO BE PRINTED AFTER THEY HAVE OCCURED MAX TIMES. C SUCH NON-FATAL MESSAGES MAY BE PRINTED LESS THAN C MAX TIMES EVEN IF THEY OCCUR MAX TIMES, IF ERROR C SUPPRESSION MODE (KONTRL=0) IS EVER IN EFFECT. C C DESCRIPTION OF PARAMETER C --INPUT-- C MAX - THE MAXIMUM NUMBER OF TIMES ANY ONE MESSAGE C IS TO BE PRINTED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 7 JUNE 1978 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XERMAX C C***FIRST EXECUTABLE STATEMENT XERMAX JUNK = J4SAVE(4,MAX,.TRUE.) RETURN END SUBROUTINE XERPRT(MESSG,NMESSG) XERPRT 2 C***BEGIN PROLOGUE XERPRT C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,PRINT C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C PRINTS ERROR MESSAGES C***DESCRIPTION C ABSTRACT C PRINT THE HOLLERITH MESSAGE IN MESSG, OF LENGTH NMESSG, C ON EACH FILE INDICATED BY XGETUA. C C PREPARE FORMAT C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED XGETUA,S88FMT,I1MACH C***END PROLOGUE XERPRT C INTEGER F(10),G(14),LUN(5) DIMENSION MESSG(NMESSG) DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10) 1 / 1H( ,1H1 ,1HX ,1H, ,1H ,1H ,1HA ,1H ,1H ,1H) / DATA G(1),G(2),G(3),G(4),G(5),G(6),G(7),G(8),G(9),G(10) 1 / 1H( ,1H1 ,1HX ,1H ,1H ,1H ,1H ,1H ,1H ,1H / DATA G(11),G(12),G(13),G(14) 1 / 1H ,1H ,1H ,1H) / DATA LA/1HA/,LCOM/1H,/,LBLANK/1H / C PREPARE FORMAT FOR WHOLE LINES C***FIRST EXECUTABLE STATEMENT XERPRT NCHAR = I1MACH(6) NFIELD = 72/NCHAR CALL S88FMT(2,NFIELD,F(5)) CALL S88FMT(2,NCHAR,F(8)) C PREPARE FORMAT FOR LAST, PARTIAL LINE, IF NEEDED NCHARL = NFIELD*NCHAR NLINES = NMESSG/NCHARL NWORD = NLINES*NFIELD NCHREM = NMESSG - NLINES*NCHARL IF (NCHREM.LE.0) GO TO 40 DO 10 I=4,13 10 G(I) = LBLANK NFIELD = NCHREM/NCHAR IF (NFIELD.LE.0) GO TO 20 C PREPARE WHOLE WORD FIELDS G(4) = LCOM CALL S88FMT(2,NFIELD,G(5)) G(7) = LA CALL S88FMT(2,NCHAR,G(8)) 20 CONTINUE NCHLST = MOD(NCHREM,NCHAR) IF (NCHLST.LE.0) GO TO 30 C PREPARE PARTIAL WORD FIELD G(10) = LCOM G(11) = LA CALL S88FMT(2,NCHLST,G(12)) 30 CONTINUE 40 CONTINUE C PRINT THE MESSAGE NWORD1 = NWORD+1 NWORD2 = (NMESSG+NCHAR-1)/NCHAR CALL XGETUA(LUN,NUNIT) DO 50 KUNIT = 1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) IF (NWORD.GT.0) WRITE (IUNIT,F) (MESSG(I),I=1,NWORD) IF (NCHREM.GT.0) WRITE (IUNIT,G) (MESSG(I),I=NWORD1,NWORD2) 50 CONTINUE RETURN END SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL) XERROR 2 C***BEGIN PROLOGUE XERROR C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C PROCESSES AN ERROR (DIAGNOSTIC) MESSAGE C***DESCRIPTION C ABSTRACT C XERROR PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER C DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE C OF THE LIBRARY ERROR CONTROL FLAG, KONTRL. C (SEE SUBROUTINE XSETF FOR DETAILS.) C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED, CONTAINING C NO MORE THAN 72 CHARACTERS. C NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG. C NERR - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE. C NERR MUST NOT BE ZERO. C LEVEL - ERROR CATEGORY. C =2 MEANS THIS IS AN UNCONDITIONALLY FATAL ERROR. C =1 MEANS THIS IS A RECOVERABLE ERROR. (I.E., IT IS C NON-FATAL IF XSETF HAS BEEN APPROPRIATELY CALLED.) C =0 MEANS THIS IS A WARNING MESSAGE ONLY. C =-1 MEANS THIS IS A WARNING MESSAGE WHICH IS TO BE C PRINTED AT MOST ONCE, REGARDLESS OF HOW MANY C TIMES THIS CALL IS EXECUTED. C C EXAMPLES C CALL XERROR(23HSMOOTH -- NUM WAS ZERO.,23,1,2) C CALL XERROR(43HINTEG -- LESS THAN FULL ACCURACY ACHIEVED., C 43,2,1) C CALL XERROR(65HROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL C 1 FULLY COLLAPSED.,65,3,0) C CALL XERROR(39HEXP -- UNDERFLOWS BEING SET TO ZERO.,39,1,-1) C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 7 FEB 1979 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED XERRWV C***END PROLOGUE XERROR C DIMENSION MESSG(NMESSG) C***FIRST EXECUTABLE STATEMENT XERROR CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.) RETURN END SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) XERRWV 2 C***BEGIN PROLOGUE XERRWV C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,PRINT WITH VALUES C***DATE WRITTEN MARCH 19, 1980 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C PROCESSES ERROR MESSAGE ALLOWING 2 INTEGER AND TWO REAL VALUES C TO BE INCLUDED IN THE MESSAGE. C***DESCRIPTION C ABSTRACT C XERRWV PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER C DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE C OF THE LIBRARY ERROR CONTROL FLAG, KONTRL. C (SEE SUBROUTINE XSETF FOR DETAILS.) C IN ADDITION, UP TO TWO INTEGER VALUES AND TWO REAL C VALUES MAY BE PRINTED ALONG WITH THE MESSAGE. C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED. C NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG. C NERR - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE. C NERR MUST NOT BE ZERO. C LEVEL - ERROR CATEGORY. C =2 MEANS THIS IS AN UNCONDITIONALLY FATAL ERROR. C =1 MEANS THIS IS A RECOVERABLE ERROR. (I.E., IT IS C NON-FATAL IF XSETF HAS BEEN APPROPRIATELY CALLED.) C =0 MEANS THIS IS A WARNING MESSAGE ONLY. C =-1 MEANS THIS IS A WARNING MESSAGE WHICH IS TO BE C PRINTED AT MOST ONCE, REGARDLESS OF HOW MANY C TIMES THIS CALL IS EXECUTED. C NI - NUMBER OF INTEGER VALUES TO BE PRINTED. (O TO 2) C I1 - FIRST INTEGER VALUE. C I2 - SECOND INTEGER VALUE. C NR - NUMBER OF REAL VALUES TO BE PRINTED. (0 TO 2) C R1 - FIRST REAL VALUE. C R2 - SECOND REAL VALUE. C C EXAMPLES C CALL XERROR(29HSMOOTH -- NUM (=I1) WAS ZERO.,29,1,2, C 1 1,NUM,0,0,0.,0.) C CALL XERRWV(54HQUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM C 1 (R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN) C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 19 MAR 1980 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT C J4SAVE C***END PROLOGUE XERRWV C DIMENSION MESSG(NMESSG),LUN(5) C***FIRST EXECUTABLE STATEMENT XERRWV LKNTRL = J4SAVE(2,0,.FALSE.) MAXMES = J4SAVE(4,0,.FALSE.) C CHECK FOR VALID INPUT IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND. 1 (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT(23HXERROR -- INVALID INPUT,23) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV(1H ,0,0,0,KDUMMY) CALL XERABT(23HXERROR -- INVALID INPUT,23) RETURN 10 CONTINUE C RECORD MESSAGE JUNK = J4SAVE(1,NERR,.TRUE.) CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT) C LET USER OVERRIDE LFIRST = MESSG(1) LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL) C RESET TO ORIGINAL VALUES LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL LKNTRL = MAX0(-2,MIN0(2,LKNTRL)) MKNTRL = IABS(LKNTRL) C DECIDE WHETHER TO PRINT MESSAGE IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100 IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES))) 1.OR.((LLEVEL.EQ.0) .AND.(KOUNT.GT.MAXMES)) 2.OR.((LLEVEL.EQ.1) .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1)) 3.OR.((LLEVEL.EQ.2) .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100 IF (LKNTRL.LE.0) GO TO 20 CALL XERPRT(1H ,1) C INTRODUCTION IF (LLEVEL.EQ.(-1)) CALL XERPRT 1(57HWARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.,57) IF (LLEVEL.EQ.0) CALL XERPRT(13HWARNING IN...,13) IF (LLEVEL.EQ.1) CALL XERPRT 1 (23HRECOVERABLE ERROR IN...,23) IF (LLEVEL.EQ.2) CALL XERPRT(17HFATAL ERROR IN...,17) 20 CONTINUE C MESSAGE CALL XERPRT(MESSG,LMESSG) CALL XGETUA(LUN,NUNIT) DO 50 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) IF (NI.GE.1) WRITE (IUNIT,22) I1 IF (NI.GE.2) WRITE (IUNIT,23) I2 IF (NR.GE.1) WRITE (IUNIT,24) R1 IF (NR.GE.2) WRITE (IUNIT,25) R2 22 FORMAT (11X,21HIN ABOVE MESSAGE, I1=,I10) 23 FORMAT (11X,21HIN ABOVE MESSAGE, I2=,I10) 24 FORMAT (11X,21HIN ABOVE MESSAGE, R1=,E20.10) 25 FORMAT (11X,21HIN ABOVE MESSAGE, R2=,E20.10) IF (LKNTRL.LE.0) GO TO 40 C ERROR NUMBER WRITE (IUNIT,30) LERR 30 FORMAT (15H ERROR NUMBER =,I10) 40 CONTINUE 50 CONTINUE C TRACE-BACK CALL FDUMP 100 CONTINUE IFATAL = 0 IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2))) 1IFATAL = 1 C QUIT HERE IF MESSAGE IS NOT FATAL IF (IFATAL.LE.0) RETURN IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120 C PRINT REASON FOR ABORT IF (LLEVEL.EQ.1) CALL XERPRT 1 (35HJOB ABORT DUE TO UNRECOVERED ERROR.,35) IF (LLEVEL.EQ.2) CALL XERPRT 1 (29HJOB ABORT DUE TO FATAL ERROR.,29) C PRINT ERROR SUMMARY CALL XERSAV(1H ,-1,0,0,KDUMMY) 120 CONTINUE C ABORT IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0 CALL XERABT(MESSG,LMESSG) RETURN END SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT) XERSAV 2 C***BEGIN PROLOGUE XERSAV C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,RECORD C***DATE WRITTEN MARCH 19, 1980 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C RECORDS THAT AN ERROR OCCURRED. C***DESCRIPTION C ABSTRACT C RECORD THAT THIS ERROR OCCURRED. C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG, NMESSG, NERR, LEVEL ARE AS IN XERROR, C EXCEPT THAT WHEN NMESSG=0 THE TABLES WILL BE C DUMPED AND CLEARED, AND WHEN NMESSG IS LESS THAN ZERO THE C TABLES WILL BE DUMPED AND NOT CLEARED. C --OUTPUT-- C ICOUNT WILL BE THE NUMBER OF TIMES THIS MESSAGE HAS C BEEN SEEN, OR ZERO IF THE TABLE HAS OVERFLOWED AND C DOES NOT CONTAIN THIS MESSAGE SPECIFICALLY. C WHEN NMESSG=0, ICOUNT WILL NOT BE ALTERED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 19 MAR 1980 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED XGETUA,S88FMT,I1MACH C***END PROLOGUE XERSAV C INTEGER F(17),LUN(5) DIMENSION MESSG(1) DIMENSION MESTAB(10),NERTAB(10),LEVTAB(10),KOUNT(10) DATA MESTAB(1),MESTAB(2),MESTAB(3),MESTAB(4),MESTAB(5), 1 MESTAB(6),MESTAB(7),MESTAB(8),MESTAB(9),MESTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA NERTAB(1),NERTAB(2),NERTAB(3),NERTAB(4),NERTAB(5), 1 NERTAB(6),NERTAB(7),NERTAB(8),NERTAB(9),NERTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA LEVTAB(1),LEVTAB(2),LEVTAB(3),LEVTAB(4),LEVTAB(5), 1 LEVTAB(6),LEVTAB(7),LEVTAB(8),LEVTAB(9),LEVTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5), 1 KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA KOUNTX/0/ DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10), 1 F(11),F(12),F(13),F(14),F(15),F(16),F(17) 2 /1H( ,1H1 ,1HX ,1H, ,1HA ,1H ,1H ,1H, ,1HI ,1H , 3 1H ,1H, ,1H2 ,1HI ,1H1 ,1H0 ,1H) / C***FIRST EXECUTABLE STATEMENT XERSAV IF (NMESSG.GT.0) GO TO 80 C DUMP THE TABLE IF (KOUNT(1).EQ.0) RETURN C PREPARE FORMAT NCHAR = I1MACH(6) CALL S88FMT(2,NCHAR,F(6)) NCOL = 20 - NCHAR CALL S88FMT(2,NCOL,F(10)) C PRINT TO EACH UNIT CALL XGETUA(LUN,NUNIT) DO 60 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) C PRINT TABLE HEADER WRITE (IUNIT,10) 10 FORMAT (32H0 ERROR MESSAGE SUMMARY/ 1 41H FIRST WORD NERR LEVEL COUNT) C PRINT BODY OF TABLE DO 20 I=1,10 IF (KOUNT(I).EQ.0) GO TO 30 WRITE (IUNIT,F) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I) 20 CONTINUE 30 CONTINUE C PRINT NUMBER OF OTHER ERRORS IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX 40 FORMAT (41H0OTHER ERRORS NOT INDIVIDUALLY TABULATED=,I10) WRITE (IUNIT,50) 50 FORMAT (1X) 60 CONTINUE IF (NMESSG.LT.0) RETURN C CLEAR THE ERROR TABLES DO 70 I=1,10 70 KOUNT(I) = 0 KOUNTX = 0 RETURN 80 CONTINUE 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. DO 90 I=1,10 II = I IF (KOUNT(I).EQ.0) GO TO 110 IF (MESSG(1).NE.MESTAB(I)) GO TO 90 IF (NERR.NE.NERTAB(I)) GO TO 90 IF (LEVEL.NE.LEVTAB(I)) GO TO 90 GO TO 100 90 CONTINUE C THREE POSSIBLE CASES... C TABLE IS FULL KOUNTX = KOUNTX+1 ICOUNT = 1 RETURN C MESSAGE FOUND IN TABLE 100 KOUNT(II) = KOUNT(II) + 1 ICOUNT = KOUNT(II) RETURN C EMPTY SLOT FOUND FOR NEW MESSAGE 110 MESTAB(II) = MESSG(1) NERTAB(II) = NERR LEVTAB(II) = LEVEL KOUNT(II) = 1 ICOUNT = 1 RETURN END SUBROUTINE XGETF(KONTRL) XGETF 2 C***BEGIN PROLOGUE XGETF C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,FLAG C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C RETURNS CURRENT VALUE OF ERROR CONTROL FLAG C***DESCRIPTION C ABSTRACT C XGETF RETURNS THE CURRENT VALUE OF THE ERROR CONTROL FLAG C IN KONTRL. SEE SUBROUTINE XSETF FOR FLAG VALUE MEANINGS. C (KONTRL IS AN OUTPUT PARAMETER ONLY.) C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 7 JUNE 1978 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XGETF C C***FIRST EXECUTABLE STATEMENT XGETF KONTRL = J4SAVE(2,0,.FALSE.) RETURN END SUBROUTINE XGETUA(IUNIT,N) XGETUA 2 C***BEGIN PROLOGUE XGETUA C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,UNIT NUMBER C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C RETURNS UNIT NUMBER(S) TO WHICH ERROR MESSAGES ARE BEING SENT C***DESCRIPTION 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 WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XGETUA C DIMENSION IUNIT(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 IUNIT(I) = J4SAVE(INDEX,0,.FALSE.) 30 CONTINUE RETURN END SUBROUTINE XGETUN(IUNIT) XGETUN 2 C***BEGIN PROLOGUE XGETUN C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,OUTPUT UNIT C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C RETURNS THE (FIRST) OUTPUT FILE TO WHICH MESSAGES ARE BEING SENT C***DESCRIPTION C ABSTRACT C XGETUN GETS THE (FIRST) OUTPUT FILE TO WHICH ERROR MESSAGES C ARE BEING SENT. TO FIND OUT IF MORE THAN ONE FILE IS BEING C USED, ONE MUST USE THE XGETUA ROUTINE. C C DESCRIPTION OF PARAMETER C --OUTPUT-- C IUNIT - THE LOGICAL UNIT NUMBER OF THE (FIRST) UNIT TO C WHICH ERROR MESSAGES ARE BEING SENT. C A VALUE OF ZERO MEANS THAT THE DEFAULT FILE, AS C DEFINED BY THE I1MACH ROUTINE, IS BEING USED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 23 MAY 1979 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XGETUN C C***FIRST EXECUTABLE STATEMENT XGETUN IUNIT = J4SAVE(3,0,.FALSE.) RETURN END SUBROUTINE XSETF(KONTRL) XSETF 2 C***BEGIN PROLOGUE XSETF C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,EROR,SET FLAG C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C SETS THE ERROR CONTROL FLAG C***DESCRIPTION C ABSTRACT C XSETF SETS THE ERROR CONTROL FLAG VALUE TO KONTRL. C (KONTRL IS AN INPUT PARAMETER ONLY.) C THE FOLLOWING TABLE SHOWS HOW EACH MESSAGE IS TREATED, C DEPENDING ON THE VALUES OF KONTRL AND LEVEL. (SEE XERROR C FOR DESCRIPTION OF LEVEL.) C C IF KONTRL IS ZERO OR NEGATIVE, NO INFORMATION OTHER THAN THE C MESSAGE ITSELF (INCLUDING NUMERIC VALUES, IF ANY) WILL BE C PRINTED. IF KONTRL IS POSITIVE, INTRODUCTORY MESSAGES, C TRACE-BACKS, ETC., WILL BE PRINTED IN ADDITION TO THE MESSAGE. C C IABS(KONTRL) C LEVEL 0 1 2 C VALUE C 2 FATAL FATAL FATAL C C 1 NOT PRINTED PRINTED FATAL C C 0 NOT PRINTED PRINTED PRINTED C C -1 NOT PRINTED PRINTED PRINTED C ONLY ONLY C ONCE ONCE C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 23 MAY 1979 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE,XERRWV C***END PROLOGUE XSETF C C***FIRST EXECUTABLE STATEMENT XSETF IF ((KONTRL.GE.(-2)).AND.(KONTRL.LE.2)) GO TO 10 CALL XERRWV(39HXSETF -- INVALID VALUE OF KONTRL (I1).,33,1,2, 1 1,KONTRL,0,0,0.,0.) RETURN 10 JUNK = J4SAVE(2,KONTRL,.TRUE.) RETURN END SUBROUTINE XSETUA(IUNIT,N) XSETUA 2 C***BEGIN PROLOGUE XSETUA C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,UNIT NUMBERS C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C SETS UP TO 5 UNIT NUMBERS TO WHICH MESSAGES ARE TO BE SENT C***DESCRIPTION C ABSTRACT C XSETUA MAY BE CALLED TO DECLARE A LIST OF UP TO FIVE C LOGICAL UNITS, EACH OF WHICH IS TO RECEIVE A COPY OF C EACH ERROR MESSAGE PROCESSED BY THIS PACKAGE. C THE PURPOSE OF XSETUA IS TO ALLOW SIMULTANEOUS PRINTING C OF EACH ERROR MESSAGE ON, SAY, A MAIN OUTPUT FILE, C AN INTERACTIVE TERMINAL, AND OTHER FILES SUCH AS GRAPHICS C COMMUNICATION FILES. C C DESCRIPTION OF PARAMETERS C --INPUT-- C IUNIT - AN ARRAY OF UP TO FIVE UNIT NUMBERS. C NORMALLY THESE NUMBERS SHOULD ALL BE DIFFERENT. C (BUT DUPLICATES ARE NOT PROHIBITED.) C N - THE NUMBER OF UNIT NUMBERS PROVIDED IN IUNIT. C MUST HAVE 1 .LE. N .LE. 5. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 23 MAY 1979 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE,XERRWV C***END PROLOGUE XSETUA C DIMENSION IUNIT(5) C***FIRST EXECUTABLE STATEMENT XSETUA IF ((N.GE.1).AND.(N.LE.5)) GO TO 10 CALL XERRWV(34HXSETUA -- INVALID VALUE OF N (I1).,34,1,2, 1 1,N,0,0,0.,0.) RETURN 10 CONTINUE DO 20 I=1,N INDEX = I+4 IF (I.EQ.1) INDEX = 3 JUNK = J4SAVE(INDEX,IUNIT(I),.TRUE.) 20 CONTINUE JUNK = J4SAVE(5,N,.TRUE.) RETURN END SUBROUTINE XSETUN(IUNIT) XSETUN 2 C***BEGIN PROLOGUE XSETUN C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS XERROR,ERROR,OUTPUT UNIT C***DATE WRITTEN AUGUST 1979 C***AUTHOR JONES R.E. (SLA) C***PURPOSE C SETS OUTPUT FILE TO WHICH ERROR MESSAGES ARE TO BE SENT C***DESCRIPTION C ABSTRACT C XSETUN SETS THE OUTPUT FILE TO WHICH ERROR MESSAGES ARE TO C BE SENT. ONLY ONE FILE WILL BE USED. SEE XSETUA FOR C HOW TO DECLARE MORE THAN ONE FILE. C C DESCRIPTION OF PARAMETER C --INPUT-- C IUNIT - AN INPUT PARAMETER GIVING THE LOGICAL UNIT NUMBER C TO WHICH ERROR MESSAGES ARE TO BE SENT. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C LATEST REVISION --- 7 JUNE 1978 C C***REFERENCES C JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR HANDLING C PACKAGE*, SAND78-1189, SANDIA LABORATORIES, 1978. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XSETUN C C***FIRST EXECUTABLE STATEMENT XSETUN JUNK = J4SAVE(3,IUNIT,.TRUE.) JUNK = J4SAVE(5,1,.TRUE.) RETURN END DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***REVISION DATE 820701 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS MACHINE CONSTANTS,DOUBLE PRECISION C***DATE WRITTEN 1975 C***AUTHOR FOX P.A., HALL A.D., SCHRYER N.L. (BELL LABS) C***PURPOSE C RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C C D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C D = D1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF D ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C DOUBLE-PRECISION MACHINE CONSTANTS 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 TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C C WHERE POSSIBLE, OCTAL OR HEXADECIMAL CONSTANTS HAVE BEEN USED C TO SPECIFY THE CONSTANTS EXACTLY WHICH HAS IN SOME CASES C REQUIRED THE USE OF EQUIVALENT INTEGER ARRAYS. C C***REFERENCES C FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A PORTABLE LIBRARY*, C ACM TRANSACTION ON MATHEMATICAL SOFTWARE, VOL. 4, NO. 2, C JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR 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) 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 BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C 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 C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C 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 C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00604000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37767777777777777777B / C DATA LARGE(2) / 37167777777777777777B / C C DATA RIGHT(1) / 15604000000000000000B / C DATA RIGHT(2) / 15000000000000000000B / C C DATA DIVER(1) / 15614000000000000000B / C DATA DIVER(2) / 15010000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA SMALL(1) / 200004000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 577777777777777777777B / C DATA LARGE(2) / 000007777777777777777B / C C DATA RIGHT(1) / 376414000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376424000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ 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 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 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 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 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 FORTRANS 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 FORTRANS 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 C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C 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 C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 / 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 MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR 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 C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** 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***FIRST EXECUTABLE STATEMENT D1MACH C IF (I .LT. 1 .OR. I .GT. 5) 1 CALL XERROR(25HD1MACH -- I OUT OF BOUNDS,25,1,2) C D1MACH = DMACH(I) RETURN C END INTEGER FUNCTION I1MACH(I) C***BEGIN PROLOGUE I1MACH C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS MACHINE CONSTANTS,INTEGER C***DATE WRITTEN 1975 C***AUTHOR FOX P.A.,HALL A.D.,SCHRYER N.L. (BELL LABS) C***PURPOSE C RETURNS INTEGER MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C THESE MACHINE CONSTANT ROUTINES MUST BE ACTIVATED FOR C A PARTICULAR ENVIRONMENT. C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C I1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C K = I1MACH(I) C C WHERE I=1,...,16. THE (OUTPUT) VALUE OF K ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE 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, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. ALSO, THE VALUES OF C I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY C WITH THE LOCAL OPERATING SYSTEM. C C***REFERENCES C FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A PORTABLE LIBRARY*, C ACM TRANSACTION ON MATHEMATICAL SOFTWARE, VOL. 4, NO. 2, C JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT C C EQUIVALENCE (IMACH(4),OUTPUT) 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 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) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CRAY 1 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) / 48 / C DATA IMACH(12) / -8192 / C DATA IMACH(13) / 8191 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -8192 / C DATA IMACH(16) / 8191 / 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 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 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 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 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 FORTRANS 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 FORTRANS 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 UNIVAC 1100 SERIES. FTN COMPILER C 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 UNIVAC 1100 SERIES. FOR COMPILER 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) / 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) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024/ C DATA IMACH(16) / 1023 / C C C MACHINE CONSTANTS FOR THE VAX 11/780 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***FIRST EXECUTABLE STATEMENT I1MACH C IF (I .LT. 1 .OR. I .GT. 16) GO TO 10 C I1MACH=IMACH(I) RETURN C 10 CONTINUE WRITE(OUTPUT,9000) 9000 FORMAT(39H1ERROR 1 IN I1MACH - I OUT OF BOUNDS ) C C CALL FDUMP C C STOP END REAL FUNCTION R1MACH(I) C***BEGIN PROLOGUE R1MACH C***REVISION DATE 820701 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS MACHINE CONSTANTS C***DATE WRITTEN 1979 C***AUTHOR FOX P.A.,.HALL A.D., SCHRYER N.L. (BELL LABS) C***PURPOSE C RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C C R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C A = R1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF A ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C SINGLE-PRECISION MACHINE CONSTANTS C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C R1MACH(5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C C WHERE POSSIBLE, OCTAL OR HEXADECIMAL CONSTANTS HAVE BEEN USED C TO SPECIFY THE CONSTANTS EXACTLY WHICH HAS IN SOME CASES C REQUIRED THE USE OF EQUIVALENT INTEGER ARRAYS. C C***REFERENCES C FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A PORTABLE LIBRARY*, C ACM TRANSACTION ON MATHEMATICAL SOFTWARE, VOL. 4, NO. 2, C JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE R1MACH C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C REAL RMACH(5) C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / 00014000000000000000B / C DATA RMACH(2) / 37767777777777777777B / C DATA RMACH(3) / 16404000000000000000B / C DATA RMACH(4) / 16414000000000000000B / C DATA RMACH(5) / 17164642023241175720B / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA RMACH(1) / 200004000000000000000B / C DATA RMACH(2) / 577777777777777777777B / C DATA RMACH(3) / 377214000000000000000B / C DATA RMACH(4) / 377224000000000000000B / C DATA RMACH(5) / 377774642023241175720B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC RMACH(5) C C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1),SMALL(2) / 20000000, 00000201 / C DATA LARGE(1),LARGE(2) / 37777777, 00000177 / C DATA RIGHT(1),RIGHT(2) / 20000000, 00000352 / C DATA DIVER(1),DIVER(2) / 20000000, 00000353 / C DATA LOG10(1),LOG10(2) / 23210115, 00000377 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 / C C MACHINE CONSTANTS FOR THE HP 2100 C C 3 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE91), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / 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 RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 / C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA RIGHT(1),RIGHT(2) / 13440, 0 / C DATA DIVER(1),DIVER(2) / 13568, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8347 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA RIGHT(1),RIGHT(2) / O032200, O000000 / C DATA DIVER(1),DIVER(2) / O032400, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020233 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 / C C MACHINE CONSTANTS FOR THE VAX 11/780 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) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 / C C DATA SMALL(1) / Z00000080 / C DATA LARGE(1) / ZFFFF7FFF / C DATA RIGHT(1) / Z00003480 / C DATA DIVER(1) / Z00003500 / C DATA LOG10(1) / Z209B3F9A / C C***FIRST EXECUTABLE STATEMENT R1MACH C IF (I .LT. 1 .OR. I .GT. 5) 1 CALL XERROR (25HR1MACH -- I OUT OF BOUNDS,25,1,2) C R1MACH = RMACH(I) RETURN C END