relibc_openlibm/slatec/dqnc79.f

*DECK DQNC79
      SUBROUTINE DQNC79 (FUN, A, B, ERR, ANS, IERR, K)
C***BEGIN PROLOGUE  DQNC79
C***PURPOSE  Integrate a function using a 7-point adaptive Newton-Cotes
C            quadrature rule.
C***LIBRARY   SLATEC
C***CATEGORY  H2A1A1
C***TYPE      DOUBLE PRECISION (QNC79-S, DQNC79-D)
C***KEYWORDS  ADAPTIVE QUADRATURE, INTEGRATION, NEWTON-COTES
C***AUTHOR  Kahaner, D. K., (NBS)
C           Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract  *** a DOUBLE PRECISION routine ***
C       DQNC79 is a general purpose program for evaluation of
C       one dimensional integrals of user defined functions.
C       DQNC79 will pick its own points for evaluation of the
C       integrand and these will vary from problem to problem.
C       Thus, DQNC79 is not designed to integrate over data sets.
C       Moderately smooth integrands will be integrated efficiently
C       and reliably.  For problems with strong singularities,
C       oscillations etc., the user may wish to use more sophis-
C       ticated routines such as those in QUADPACK.  One measure
C       of the reliability of DQNC79 is the output parameter K,
C       giving the number of integrand evaluations that were needed.
C
C     Description of Arguments
C
C     --Input--* FUN, A, B, ERR are DOUBLE PRECISION *
C       FUN  - name of external function to be integrated.  This name
C              must be in an EXTERNAL statement in your calling
C              program.  You must write a Fortran function to evaluate
C              FUN.  This should be of the form
C                    DOUBLE PRECISION FUNCTION FUN (X)
C              C
C              C     X can vary from A to B
C              C     FUN(X) should be finite for all X on interval.
C              C
C                    FUN = ...
C                    RETURN
C                    END
C       A    - lower limit of integration
C       B    - upper limit of integration (may be less than A)
C       ERR  - is a requested error tolerance.  Normally, pick a value
C              0 .LT. ERR .LT. 1.0D-8.
C
C     --Output--
C       ANS  - computed value of the integral.  Hopefully, ANS is
C              accurate to within ERR * integral of ABS(FUN(X)).
C       IERR - a status code
C            - Normal codes
C               1  ANS most likely meets requested error tolerance.
C              -1  A equals B, or A and B are too nearly equal to
C                  allow normal integration.  ANS is set to zero.
C            - Abnormal code
C               2  ANS probably does not meet requested error tolerance.
C       K    - the number of function evaluations actually used to do
C              the integration.  A value of K .GT. 1000 indicates a
C              difficult problem; other programs may be more efficient.
C              DQNC79 will gracefully give up if K exceeds 2000.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, I1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   790601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920218  Code redone to parallel QNC79.  (WRB)
C   930120  Increase array size 80->99, and KMX 2000->5000 for SUN -r8
C           wordlength.  (RWC)
C***END PROLOGUE  DQNC79
C     .. Scalar Arguments ..
      DOUBLE PRECISION A, ANS, B, ERR
      INTEGER IERR, K
C     .. Function Arguments ..
      DOUBLE PRECISION FUN
      EXTERNAL FUN
C     .. Local Scalars ..
      DOUBLE PRECISION AE, AREA, BANK, BLOCAL, C, CE, EE, EF, EPS, Q13,
     +                 Q7, Q7L, SQ2, TEST, TOL, VR, W1, W2, W3, W4
      INTEGER I, KML, KMX, L, LMN, LMX, NBITS, NIB, NLMN, NLMX
      LOGICAL FIRST
C     .. Local Arrays ..
      DOUBLE PRECISION AA(99), F(13), F1(99), F2(99), F3(99), F4(99),
     +                 F5(99), F6(99), F7(99), HH(99), Q7R(99), VL(99)
      INTEGER LR(99)
C     .. External Functions ..
      DOUBLE PRECISION D1MACH
      INTEGER I1MACH
      EXTERNAL D1MACH, I1MACH
C     .. External Subroutines ..
      EXTERNAL XERMSG
C     .. Intrinsic Functions ..
      INTRINSIC ABS, LOG, MAX, MIN, SIGN, SQRT
C     .. Save statement ..
      SAVE NBITS, NLMX, FIRST, SQ2, W1, W2, W3, W4
C     .. Data statements ..
      DATA KML /7/, KMX /5000/, NLMN /2/
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DQNC79
      IF (FIRST) THEN
        W1 = 41.0D0/140.0D0
        W2 = 216.0D0/140.0D0
        W3 = 27.0D0/140.0D0
        W4 = 272.0D0/140.0D0
        NBITS = D1MACH(5)*I1MACH(14)/0.30102000D0
        NLMX = MIN(99,(NBITS*4)/5)
        SQ2 = SQRT(2.0D0)
      ENDIF
      FIRST = .FALSE.
      ANS = 0.0D0
      IERR = 1
      CE = 0.0D0
      IF (A .EQ. B) GO TO 260
      LMX = NLMX
      LMN = NLMN
      IF (B .EQ. 0.0D0) GO TO 100
      IF (SIGN(1.0D0,B)*A .LE. 0.0D0) GO TO 100
      C = ABS(1.0D0-A/B)
      IF (C .GT. 0.1D0) GO TO 100
      IF (C .LE. 0.0D0) GO TO 260
      NIB = 0.5D0 - LOG(C)/LOG(2.0D0)
      LMX = MIN(NLMX,NBITS-NIB-4)
      IF (LMX .LT. 2) GO TO 260
      LMN = MIN(LMN,LMX)
  100 TOL = MAX(ABS(ERR),2.0D0**(5-NBITS))
      IF (ERR .EQ. 0.0D0) TOL = SQRT(D1MACH(4))
      EPS = TOL
      HH(1) = (B-A)/12.0D0
      AA(1) = A
      LR(1) = 1
      DO 110 I = 1,11,2
        F(I) = FUN(A+(I-1)*HH(1))
  110 CONTINUE
      BLOCAL = B
      F(13) = FUN(BLOCAL)
      K = 7
      L = 1
      AREA = 0.0D0
      Q7 = 0.0D0
      EF = 256.0D0/255.0D0
      BANK = 0.0D0
C
C     Compute refined estimates, estimate the error, etc.
C
  120 DO 130 I = 2,12,2
        F(I) = FUN(AA(L)+(I-1)*HH(L))
  130 CONTINUE
      K = K + 6
C
C     Compute left and right half estimates
C
      Q7L = HH(L)*((W1*(F(1)+F(7))+W2*(F(2)+F(6)))+
     +      (W3*(F(3)+F(5))+W4*F(4)))
      Q7R(L) = HH(L)*((W1*(F(7)+F(13))+W2*(F(8)+F(12)))+
     +         (W3*(F(9)+F(11))+W4*F(10)))
C
C     Update estimate of integral of absolute value
C
      AREA = AREA + (ABS(Q7L)+ABS(Q7R(L))-ABS(Q7))
C
C     Do not bother to test convergence before minimum refinement level
C
      IF (L .LT. LMN) GO TO 180
C
C     Estimate the error in new value for whole interval, Q13
C
      Q13 = Q7L + Q7R(L)
      EE = ABS(Q7-Q13)*EF
C
C     Compute nominal allowed error
C
      AE = EPS*AREA
C
C     Borrow from bank account, but not too much
C
      TEST = MIN(AE+0.8D0*BANK,10.0D0*AE)
C
C     Don't ask for excessive accuracy
C
      TEST = MAX(TEST,TOL*ABS(Q13),0.00003D0*TOL*AREA)
C
C     Now, did this interval pass or not?
C
      IF (EE-TEST) 150,150,170
C
C     Have hit maximum refinement level -- penalize the cumulative error
C
  140 CE = CE + (Q7-Q13)
      GO TO 160
C
C     On good intervals accumulate the theoretical estimate
C
  150 CE = CE + (Q7-Q13)/255.0D0
C
C     Update the bank account.  Don't go into debt.
C
  160 BANK = BANK + (AE-EE)
      IF (BANK .LT. 0.0D0) BANK = 0.0D0
C
C     Did we just finish a left half or a right half?
C
      IF (LR(L)) 190,190,210
C
C     Consider the left half of next deeper level
C
  170 IF (K .GT. KMX) LMX = MIN(KML,LMX)
      IF (L .GE. LMX) GO TO 140
  180 L = L + 1
      EPS = EPS*0.5D0
      IF (L .LE. 17) EF = EF/SQ2
      HH(L) = HH(L-1)*0.5D0
      LR(L) = -1
      AA(L) = AA(L-1)
      Q7 = Q7L
      F1(L) = F(7)
      F2(L) = F(8)
      F3(L) = F(9)
      F4(L) = F(10)
      F5(L) = F(11)
      F6(L) = F(12)
      F7(L) = F(13)
      F(13) = F(7)
      F(11) = F(6)
      F(9) = F(5)
      F(7) = F(4)
      F(5) = F(3)
      F(3) = F(2)
      GO TO 120
C
C     Proceed to right half at this level
C
  190 VL(L) = Q13
  200 Q7 = Q7R(L-1)
      LR(L) = 1
      AA(L) = AA(L) + 12.0D0*HH(L)
      F(1) = F1(L)
      F(3) = F2(L)
      F(5) = F3(L)
      F(7) = F4(L)
      F(9) = F5(L)
      F(11) = F6(L)
      F(13) = F7(L)
      GO TO 120
C
C     Left and right halves are done, so go back up a level
C
  210 VR = Q13
  220 IF (L .LE. 1) GO TO 250
      IF (L .LE. 17) EF = EF*SQ2
      EPS = EPS*2.0D0
      L = L - 1
      IF (LR(L)) 230,230,240
  230 VL(L) = VL(L+1) + VR
      GO TO 200
  240 VR = VL(L+1) + VR
      GO TO 220
C
C     Exit
C
  250 ANS = VR
      IF (ABS(CE) .LE. 2.0D0*TOL*AREA) GO TO 270
      IERR = 2
      CALL XERMSG ('SLATEC', 'DQNC79',
     +   'ANS is probably insufficiently accurate.', 2, 1)
      GO TO 270
  260 IERR = -1
      CALL XERMSG ('SLATEC', 'DQNC79',
     +   'A and B are too nearly equal to allow normal integration. $$'
     +   // 'ANS is set to zero and IERR to -1.', -1, -1)
  270 RETURN
      END