Logo Search packages:      
Sourcecode: octave-matcompat version File versions

expint.f

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

Generated by  Doxygen 1.6.0   Back to index