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
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)
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.
27
C Description of Arguments
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)
36
C C X can vary from A to B
37
C C FUN(X) should be finite for all X on interval.
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.
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
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.
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.
63
C***ROUTINES CALLED D1MACH, I1MACH, XERMSG
64
C***REVISION HISTORY (YYMMDD)
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
74
C***END PROLOGUE DQNC79
75
C .. Scalar Arguments ..
76
DOUBLE PRECISION A, ANS, B, ERR
78
C .. Function Arguments ..
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
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)
90
C .. External Functions ..
91
DOUBLE PRECISION D1MACH
93
EXTERNAL D1MACH, I1MACH
94
C .. External Subroutines ..
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/
103
C***FIRST EXECUTABLE STATEMENT DQNC79
109
NBITS = D1MACH(5)*I1MACH(14)/0.30102000D0
110
NLMX = MIN(99,(NBITS*4)/5)
117
IF (A .EQ. B) GO TO 260
120
IF (B .EQ. 0.0D0) GO TO 100
121
IF (SIGN(1.0D0,B)*A .LE. 0.0D0) GO TO 100
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
129
100 TOL = MAX(ABS(ERR),2.0D0**(5-NBITS))
130
IF (ERR .EQ. 0.0D0) TOL = SQRT(D1MACH(4))
136
F(I) = FUN(A+(I-1)*HH(1))
147
C Compute refined estimates, estimate the error, etc.
149
120 DO 130 I = 2,12,2
150
F(I) = FUN(AA(L)+(I-1)*HH(L))
154
C Compute left and right half estimates
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)))
161
C Update estimate of integral of absolute value
163
AREA = AREA + (ABS(Q7L)+ABS(Q7R(L))-ABS(Q7))
165
C Do not bother to test convergence before minimum refinement level
167
IF (L .LT. LMN) GO TO 180
169
C Estimate the error in new value for whole interval, Q13
174
C Compute nominal allowed error
178
C Borrow from bank account, but not too much
180
TEST = MIN(AE+0.8D0*BANK,10.0D0*AE)
182
C Don't ask for excessive accuracy
184
TEST = MAX(TEST,TOL*ABS(Q13),0.00003D0*TOL*AREA)
186
C Now, did this interval pass or not?
188
IF (EE-TEST) 150,150,170
190
C Have hit maximum refinement level -- penalize the cumulative error
192
140 CE = CE + (Q7-Q13)
195
C On good intervals accumulate the theoretical estimate
197
150 CE = CE + (Q7-Q13)/255.0D0
199
C Update the bank account. Don't go into debt.
201
160 BANK = BANK + (AE-EE)
202
IF (BANK .LT. 0.0D0) BANK = 0.0D0
204
C Did we just finish a left half or a right half?
206
IF (LR(L)) 190,190,210
208
C Consider the left half of next deeper level
210
170 IF (K .GT. KMX) LMX = MIN(KML,LMX)
211
IF (L .GE. LMX) GO TO 140
214
IF (L .LE. 17) EF = EF/SQ2
215
HH(L) = HH(L-1)*0.5D0
234
C Proceed to right half at this level
239
AA(L) = AA(L) + 12.0D0*HH(L)
249
C Left and right halves are done, so go back up a level
252
220 IF (L .LE. 1) GO TO 250
253
IF (L .LE. 17) EF = EF*SQ2
256
IF (LR(L)) 230,230,240
257
230 VL(L) = VL(L+1) + VR
259
240 VR = VL(L+1) + VR
265
IF (ABS(CE) .LE. 2.0D0*TOL*AREA) GO TO 270
267
CALL XERMSG ('SLATEC', 'DQNC79',
268
+ 'ANS is probably insufficiently accurate.', 2, 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)
added
removed
src/numerical/slatec/fortran/dqnc79.f
