src/numerical/slatec/fortran/dqnc79.f

   1 *DECK DQNC79
   2       SUBROUTINE DQNC79 (FUN, A, B, ERR, ANS, IERR, K)
   3 C***BEGIN PROLOGUE  DQNC79
   4 C***PURPOSE  Integrate a function using a 7-point adaptive Newton-Cotes
   5 C            quadrature rule.
   6 C***LIBRARY   SLATEC
   7 C***CATEGORY  H2A1A1
   8 C***TYPE      DOUBLE PRECISION (QNC79-S, DQNC79-D)
   9 C***KEYWORDS  ADAPTIVE QUADRATURE, INTEGRATION, NEWTON-COTES
  10 C***AUTHOR  Kahaner, D. K., (NBS)
  11 C           Jones, R. E., (SNLA)
  12 C***DESCRIPTION
  13 C
  14 C     Abstract  *** a DOUBLE PRECISION routine ***
  15 C       DQNC79 is a general purpose program for evaluation of
  16 C       one dimensional integrals of user defined functions.
  17 C       DQNC79 will pick its own points for evaluation of the
  18 C       integrand and these will vary from problem to problem.
  19 C       Thus, DQNC79 is not designed to integrate over data sets.
  20 C       Moderately smooth integrands will be integrated efficiently
  21 C       and reliably.  For problems with strong singularities,
  22 C       oscillations etc., the user may wish to use more sophis-
  23 C       ticated routines such as those in QUADPACK.  One measure
  24 C       of the reliability of DQNC79 is the output parameter K,
  25 C       giving the number of integrand evaluations that were needed.
  26 C
  27 C     Description of Arguments
  28 C
  29 C     --Input--* FUN, A, B, ERR are DOUBLE PRECISION *
  30 C       FUN  - name of external function to be integrated.  This name
  31 C              must be in an EXTERNAL statement in your calling
  32 C              program.  You must write a Fortran function to evaluate
  33 C              FUN.  This should be of the form
  34 C                    DOUBLE PRECISION FUNCTION FUN (X)
  35 C              C
  36 C              C     X can vary from A to B
  37 C              C     FUN(X) should be finite for all X on interval.
  38 C              C
  39 C                    FUN = ...
  40 C                    RETURN
  41 C                    END
  42 C       A    - lower limit of integration
  43 C       B    - upper limit of integration (may be less than A)
  44 C       ERR  - is a requested error tolerance.  Normally, pick a value
  45 C              0 .LT. ERR .LT. 1.0D-8.
  46 C
  47 C     --Output--
  48 C       ANS  - computed value of the integral.  Hopefully, ANS is
  49 C              accurate to within ERR * integral of ABS(FUN(X)).
  50 C       IERR - a status code
  51 C            - Normal codes
  52 C               1  ANS most likely meets requested error tolerance.
  53 C              -1  A equals B, or A and B are too nearly equal to
  54 C                  allow normal integration.  ANS is set to zero.
  55 C            - Abnormal code
  56 C               2  ANS probably does not meet requested error tolerance.
  57 C       K    - the number of function evaluations actually used to do
  58 C              the integration.  A value of K .GT. 1000 indicates a
  59 C              difficult problem; other programs may be more efficient.
  60 C              DQNC79 will gracefully give up if K exceeds 2000.
  61 C
  62 C***REFERENCES  (NONE)
  63 C***ROUTINES CALLED  D1MACH, I1MACH, XERMSG
  64 C***REVISION HISTORY  (YYMMDD)
  65 C   790601  DATE WRITTEN
  66 C   890531  Changed all specific intrinsics to generic.  (WRB)
  67 C   890911  Removed unnecessary intrinsics.  (WRB)
  68 C   890911  REVISION DATE from Version 3.2
  69 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  70 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  71 C   920218  Code redone to parallel QNC79.  (WRB)
  72 C   930120  Increase array size 80->99, and KMX 2000->5000 for SUN -r8
  73 C           wordlength.  (RWC)
  74 C***END PROLOGUE  DQNC79
  75 C     .. Scalar Arguments ..
  76       DOUBLE PRECISION A, ANS, B, ERR
  77       INTEGER IERR, K
  78 C     .. Function Arguments ..
  79       DOUBLE PRECISION FUN
  80       EXTERNAL FUN
  81 C     .. Local Scalars ..
  82       DOUBLE PRECISION AE, AREA, BANK, BLOCAL, C, CE, EE, EF, EPS, Q13,
  83      +                 Q7, Q7L, SQ2, TEST, TOL, VR, W1, W2, W3, W4
  84       INTEGER I, KML, KMX, L, LMN, LMX, NBITS, NIB, NLMN, NLMX
  85       LOGICAL FIRST
  86 C     .. Local Arrays ..
  87       DOUBLE PRECISION AA(99), F(13), F1(99), F2(99), F3(99), F4(99),
  88      +                 F5(99), F6(99), F7(99), HH(99), Q7R(99), VL(99)
  89       INTEGER LR(99)
  90 C     .. External Functions ..
  91       DOUBLE PRECISION D1MACH
  92       INTEGER I1MACH
  93       EXTERNAL D1MACH, I1MACH
  94 C     .. External Subroutines ..
  95       EXTERNAL XERMSG
  96 C     .. Intrinsic Functions ..
  97       INTRINSIC ABS, LOG, MAX, MIN, SIGN, SQRT
  98 C     .. Save statement ..
  99       SAVE NBITS, NLMX, FIRST, SQ2, W1, W2, W3, W4
 100 C     .. Data statements ..
 101       DATA KML /7/, KMX /5000/, NLMN /2/
 102       DATA FIRST /.TRUE./
 103 C***FIRST EXECUTABLE STATEMENT  DQNC79
 104       IF (FIRST) THEN
 105         W1 = 41.0D0/140.0D0
 106         W2 = 216.0D0/140.0D0
 107         W3 = 27.0D0/140.0D0
 108         W4 = 272.0D0/140.0D0
 109         NBITS = D1MACH(5)*I1MACH(14)/0.30102000D0
 110         NLMX = MIN(99,(NBITS*4)/5)
 111         SQ2 = SQRT(2.0D0)
 112       ENDIF
 113       FIRST = .FALSE.
 114       ANS = 0.0D0
 115       IERR = 1
 116       CE = 0.0D0
 117       IF (A .EQ. B) GO TO 260
 118       LMX = NLMX
 119       LMN = NLMN
 120       IF (B .EQ. 0.0D0) GO TO 100
 121       IF (SIGN(1.0D0,B)*A .LE. 0.0D0) GO TO 100
 122       C = ABS(1.0D0-A/B)
 123       IF (C .GT. 0.1D0) GO TO 100
 124       IF (C .LE. 0.0D0) GO TO 260
 125       NIB = 0.5D0 - LOG(C)/LOG(2.0D0)
 126       LMX = MIN(NLMX,NBITS-NIB-4)
 127       IF (LMX .LT. 2) GO TO 260
 128       LMN = MIN(LMN,LMX)
 129   100 TOL = MAX(ABS(ERR),2.0D0**(5-NBITS))
 130       IF (ERR .EQ. 0.0D0) TOL = SQRT(D1MACH(4))
 131       EPS = TOL
 132       HH(1) = (B-A)/12.0D0
 133       AA(1) = A
 134       LR(1) = 1
 135       DO 110 I = 1,11,2
 136         F(I) = FUN(A+(I-1)*HH(1))
 137   110 CONTINUE
 138       BLOCAL = B
 139       F(13) = FUN(BLOCAL)
 140       K = 7
 141       L = 1
 142       AREA = 0.0D0
 143       Q7 = 0.0D0
 144       EF = 256.0D0/255.0D0
 145       BANK = 0.0D0
 146 C
 147 C     Compute refined estimates, estimate the error, etc.
 148 C
 149   120 DO 130 I = 2,12,2
 150         F(I) = FUN(AA(L)+(I-1)*HH(L))
 151   130 CONTINUE
 152       K = K + 6
 153 C
 154 C     Compute left and right half estimates
 155 C
 156       Q7L = HH(L)*((W1*(F(1)+F(7))+W2*(F(2)+F(6)))+
 157      +      (W3*(F(3)+F(5))+W4*F(4)))
 158       Q7R(L) = HH(L)*((W1*(F(7)+F(13))+W2*(F(8)+F(12)))+
 159      +         (W3*(F(9)+F(11))+W4*F(10)))
 160 C
 161 C     Update estimate of integral of absolute value
 162 C
 163       AREA = AREA + (ABS(Q7L)+ABS(Q7R(L))-ABS(Q7))
 164 C
 165 C     Do not bother to test convergence before minimum refinement level
 166 C
 167       IF (L .LT. LMN) GO TO 180
 168 C
 169 C     Estimate the error in new value for whole interval, Q13
 170 C
 171       Q13 = Q7L + Q7R(L)
 172       EE = ABS(Q7-Q13)*EF
 173 C
 174 C     Compute nominal allowed error
 175 C
 176       AE = EPS*AREA
 177 C
 178 C     Borrow from bank account, but not too much
 179 C
 180       TEST = MIN(AE+0.8D0*BANK,10.0D0*AE)
 181 C
 182 C     Don't ask for excessive accuracy
 183 C
 184       TEST = MAX(TEST,TOL*ABS(Q13),0.00003D0*TOL*AREA)
 185 C
 186 C     Now, did this interval pass or not?
 187 C
 188       IF (EE-TEST) 150,150,170
 189 C
 190 C     Have hit maximum refinement level -- penalize the cumulative error
 191 C
 192   140 CE = CE + (Q7-Q13)
 193       GO TO 160
 194 C
 195 C     On good intervals accumulate the theoretical estimate
 196 C
 197   150 CE = CE + (Q7-Q13)/255.0D0
 198 C
 199 C     Update the bank account.  Don't go into debt.
 200 C
 201   160 BANK = BANK + (AE-EE)
 202       IF (BANK .LT. 0.0D0) BANK = 0.0D0
 203 C
 204 C     Did we just finish a left half or a right half?
 205 C
 206       IF (LR(L)) 190,190,210
 207 C
 208 C     Consider the left half of next deeper level
 209 C
 210   170 IF (K .GT. KMX) LMX = MIN(KML,LMX)
 211       IF (L .GE. LMX) GO TO 140
 212   180 L = L + 1
 213       EPS = EPS*0.5D0
 214       IF (L .LE. 17) EF = EF/SQ2
 215       HH(L) = HH(L-1)*0.5D0
 216       LR(L) = -1
 217       AA(L) = AA(L-1)
 218       Q7 = Q7L
 219       F1(L) = F(7)
 220       F2(L) = F(8)
 221       F3(L) = F(9)
 222       F4(L) = F(10)
 223       F5(L) = F(11)
 224       F6(L) = F(12)
 225       F7(L) = F(13)
 226       F(13) = F(7)
 227       F(11) = F(6)
 228       F(9) = F(5)
 229       F(7) = F(4)
 230       F(5) = F(3)
 231       F(3) = F(2)
 232       GO TO 120
 233 C
 234 C     Proceed to right half at this level
 235 C
 236   190 VL(L) = Q13
 237   200 Q7 = Q7R(L-1)
 238       LR(L) = 1
 239       AA(L) = AA(L) + 12.0D0*HH(L)
 240       F(1) = F1(L)
 241       F(3) = F2(L)
 242       F(5) = F3(L)
 243       F(7) = F4(L)
 244       F(9) = F5(L)
 245       F(11) = F6(L)
 246       F(13) = F7(L)
 247       GO TO 120
 248 C
 249 C     Left and right halves are done, so go back up a level
 250 C
 251   210 VR = Q13
 252   220 IF (L .LE. 1) GO TO 250
 253       IF (L .LE. 17) EF = EF*SQ2
 254       EPS = EPS*2.0D0
 255       L = L - 1
 256       IF (LR(L)) 230,230,240
 257   230 VL(L) = VL(L+1) + VR
 258       GO TO 200
 259   240 VR = VL(L+1) + VR
 260       GO TO 220
 261 C
 262 C     Exit
 263 C
 264   250 ANS = VR
 265       IF (ABS(CE) .LE. 2.0D0*TOL*AREA) GO TO 270
 266       IERR = 2
 267       CALL XERMSG ('SLATEC', 'DQNC79',
 268      +   'ANS is probably insufficiently accurate.', 2, 1)
 269       GO TO 270
 270   260 IERR = -1
 271       CALL XERMSG ('SLATEC', 'DQNC79',
 272      +   'A and B are too nearly equal to allow normal integration. $$'
 273      +   // 'ANS is set to zero and IERR to -1.', -1, -1)
 274   270 RETURN
 275       END