2
SUBROUTINE QAWC (F, A, B, C, EPSABS, EPSREL, RESULT, ABSERR,
3
+ NEVAL, IER, LIMIT, LENW, LAST, IWORK, WORK)
4
C***BEGIN PROLOGUE QAWC
5
C***PURPOSE The routine calculates an approximation result to a
6
C Cauchy principal value I = INTEGRAL of F*W over (A,B)
7
C (W(X) = 1/((X-C), C.NE.A, C.NE.B), hopefully satisfying
8
C following claim for accuracy
9
C ABS(I-RESULT).LE.MAX(EPSABE,EPSREL*ABS(I)).
10
C***LIBRARY SLATEC (QUADPACK)
11
C***CATEGORY H2A2A1, J4
12
C***TYPE SINGLE PRECISION (QAWC-S, DQAWC-D)
13
C***KEYWORDS AUTOMATIC INTEGRATOR, CAUCHY PRINCIPAL VALUE,
14
C CLENSHAW-CURTIS METHOD, GLOBALLY ADAPTIVE, QUADPACK,
15
C QUADRATURE, SPECIAL-PURPOSE
16
C***AUTHOR Piessens, Robert
17
C Applied Mathematics and Programming Division
20
C Applied Mathematics and Programming Division
24
C Computation of a Cauchy principal value
25
C Standard fortran subroutine
32
C Function subprogram defining the integrand
33
C Function F(X). The actual name for F needs to be
34
C declared E X T E R N A L in the driver program.
37
C Under limit of integration
40
C Upper limit of integration
42
C C - Parameter in the weight function, C.NE.A, C.NE.B.
43
C If C = A or C = B, the routine will end with
47
C Absolute accuracy requested
49
C Relative accuracy requested
51
C and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
52
C the routine will end with IER = 6.
56
C Approximation to the integral
59
C Estimate or the modulus of the absolute error,
60
C Which should equal or exceed ABS(I-RESULT)
63
C Number of integrand evaluations
66
C IER = 0 Normal and reliable termination of the
67
C routine. It is assumed that the requested
68
C accuracy has been achieved.
69
C IER.GT.0 Abnormal termination of the routine
70
C the estimates for integral and error are
71
C less reliable. It is assumed that the
72
C requested accuracy has not been achieved.
74
C IER = 1 Maximum number of subdivisions allowed
75
C has been achieved. One can allow more sub-
76
C divisions by increasing the value of LIMIT
77
C (and taking the according dimension
78
C adjustments into account). However, if
79
C this yields no improvement it is advised
80
C to analyze the integrand in order to
81
C determine the integration difficulties.
82
C If the position of a local difficulty
83
C can be determined (e.g. SINGULARITY,
84
C DISCONTINUITY within the interval) one
85
C will probably gain from splitting up the
86
C interval at this point and calling
87
C appropriate integrators on the subranges.
88
C = 2 The occurrence of roundoff error is detec-
89
C ted, which prevents the requested
90
C tolerance from being achieved.
91
C = 3 Extremely bad integrand behaviour occurs
92
C at some points of the integration
94
C = 6 The input is invalid, because
97
C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
98
C or LIMIT.LT.1 or LENW.LT.LIMIT*4.
99
C RESULT, ABSERR, NEVAL, LAST are set to
100
C zero. Except when LENW or LIMIT is
101
C invalid, IWORK(1), WORK(LIMIT*2+1) and
102
C WORK(LIMIT*3+1) are set to zero, WORK(1)
103
C is set to A and WORK(LIMIT+1) to B.
105
C DIMENSIONING PARAMETERS
107
C Dimensioning parameter for IWORK
108
C LIMIT determines the maximum number of subintervals
109
C in the partition of the given integration interval
111
C If LIMIT.LT.1, the routine will end with IER = 6.
114
C Dimensioning parameter for WORK
115
C LENW must be at least LIMIT*4.
116
C If LENW.LT.LIMIT*4, the routine will end with
120
C On return, LAST equals the number of subintervals
121
C produced in the subdivision process, which
122
C determines the number of significant elements
123
C actually in the WORK ARRAYS.
127
C Vector of dimension at least LIMIT, the first K
128
C elements of which contain pointers
129
C to the error estimates over the subintervals,
130
C such that WORK(LIMIT*3+IWORK(1)), ... ,
131
C WORK(LIMIT*3+IWORK(K)) form a decreasing
132
C sequence, with K = LAST if LAST.LE.(LIMIT/2+2),
133
C and K = LIMIT+1-LAST otherwise
136
C Vector of dimension at least LENW
138
C WORK(1), ..., WORK(LAST) contain the left
139
C end points of the subintervals in the
140
C partition of (A,B),
141
C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) contain
142
C the right end points,
143
C WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) contain
144
C the integral approximations over the subintervals,
145
C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
146
C contain the error estimates.
148
C***REFERENCES (NONE)
149
C***ROUTINES CALLED QAWCE, XERMSG
150
C***REVISION HISTORY (YYMMDD)
151
C 800101 DATE WRITTEN
152
C 890831 Modified array declarations. (WRB)
153
C 890831 REVISION DATE from Version 3.2
154
C 891214 Prologue converted to Version 4.0 format. (BAB)
155
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
156
C***END PROLOGUE QAWC
158
REAL A,ABSERR,B,C,EPSABS,EPSREL,F,RESULT,WORK
159
INTEGER IER,IWORK,LENW,LIMIT,LVL,L1,L2,L3,NEVAL
161
DIMENSION IWORK(*),WORK(*)
165
C CHECK VALIDITY OF LIMIT AND LENW.
167
C***FIRST EXECUTABLE STATEMENT QAWC
173
IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10
175
C PREPARE CALL FOR QAWCE.
180
CALL QAWCE(F,A,B,C,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,NEVAL,IER,
181
1 WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST)
183
C CALL ERROR HANDLER IF NECESSARY.
186
10 IF(IER.EQ.6) LVL = 1
187
IF (IER .NE. 0) CALL XERMSG ('SLATEC', 'QAWC',
188
+ 'ABNORMAL RETURN', IER, LVL)