GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dqagp.f
Go to the documentation of this file.
1  SUBROUTINE dqagp(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR,
2  * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK)
3 C***BEGIN PROLOGUE DQAGP
4 C***DATE WRITTEN 800101 (YYMMDD)
5 C***REVISION DATE 830518 (YYMMDD)
6 C***CATEGORY NO. H2A2A1
7 C***KEYWORDS AUTOMATIC INTEGRATOR, GENERAL-PURPOSE,
8 C SINGULARITIES AT USER SPECIFIED POINTS,
9 C EXTRAPOLATION, GLOBALLY ADAPTIVE
10 C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV - K.U.LEUVEN
11 C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
12 C***PURPOSE THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
13 C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
14 C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
15 C BREAK POINTS OF THE INTEGRATION INTERVAL, WHERE LOCAL
16 C DIFFICULTIES OF THE INTEGRAND MAY OCCUR (E.G.
17 C SINGULARITIES, DISCONTINUITIES), ARE PROVIDED BY THE USER.
18 C***DESCRIPTION
19 C
20 C COMPUTATION OF A DEFINITE INTEGRAL
21 C STANDARD FORTRAN SUBROUTINE
22 C DOUBLE PRECISION VERSION
23 C
24 C PARAMETERS
25 C ON ENTRY
26 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND
27 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
28 C DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
29 C
30 C A - DOUBLE PRECISION
31 C LOWER LIMIT OF INTEGRATION
32 C
33 C B - DOUBLE PRECISION
34 C UPPER LIMIT OF INTEGRATION
35 C
36 C NPTS2 - INTEGER
37 C NUMBER EQUAL TO TWO MORE THAN THE NUMBER OF
38 C USER-SUPPLIED BREAK POINTS WITHIN THE INTEGRATION
39 C RANGE, NPTS.GE.2.
40 C IF NPTS2.LT.2, THE ROUTINE WILL END WITH IER = 6.
41 C
42 C POINTS - DOUBLE PRECISION
43 C VECTOR OF DIMENSION NPTS2, THE FIRST (NPTS2-2)
44 C ELEMENTS OF WHICH ARE THE USER PROVIDED BREAK
45 C POINTS. IF THESE POINTS DO NOT CONSTITUTE AN
46 C ASCENDING SEQUENCE THERE WILL BE AN AUTOMATIC
47 C SORTING.
48 C
49 C EPSABS - DOUBLE PRECISION
50 C ABSOLUTE ACCURACY REQUESTED
51 C EPSREL - DOUBLE PRECISION
52 C RELATIVE ACCURACY REQUESTED
53 C IF EPSABS.LE.0
54 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
55 C THE ROUTINE WILL END WITH IER = 6.
56 C
57 C ON RETURN
58 C RESULT - DOUBLE PRECISION
59 C APPROXIMATION TO THE INTEGRAL
60 C
61 C ABSERR - DOUBLE PRECISION
62 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
63 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
64 C
65 C NEVAL - INTEGER
66 C NUMBER OF INTEGRAND EVALUATIONS
67 C
68 C IER - INTEGER
69 C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
70 C ROUTINE. IT IS ASSUMED THAT THE REQUESTED
71 C ACCURACY HAS BEEN ACHIEVED.
72 C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE.
73 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE
74 C LESS RELIABLE. IT IS ASSUMED THAT THE
75 C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
76 C ERROR MESSAGES
77 C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
78 C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
79 C SUBDIVISIONS BY INCREASING THE VALUE OF
80 C LIMIT (AND TAKING THE ACCORDING DIMENSION
81 C ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
82 C THIS YIELDS NO IMPROVEMENT IT IS ADVISED
83 C TO ANALYZE THE INTEGRAND IN ORDER TO
84 C DETERMINE THE INTEGRATION DIFFICULTIES. IF
85 C THE POSITION OF A LOCAL DIFFICULTY CAN BE
86 C DETERMINED (I.E. SINGULARITY,
87 C DISCONTINUITY WITHIN THE INTERVAL), IT
88 C SHOULD BE SUPPLIED TO THE ROUTINE AS AN
89 C ELEMENT OF THE VECTOR POINTS. IF NECESSARY
90 C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
91 C MUST BE USED, WHICH IS DESIGNED FOR
92 C HANDLING THE TYPE OF DIFFICULTY INVOLVED.
93 C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
94 C DETECTED, WHICH PREVENTS THE REQUESTED
95 C TOLERANCE FROM BEING ACHIEVED.
96 C THE ERROR MAY BE UNDER-ESTIMATED.
97 C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
98 C AT SOME POINTS OF THE INTEGRATION
99 C INTERVAL.
100 C = 4 THE ALGORITHM DOES NOT CONVERGE.
101 C ROUNDOFF ERROR IS DETECTED IN THE
102 C EXTRAPOLATION TABLE.
103 C IT IS PRESUMED THAT THE REQUESTED
104 C TOLERANCE CANNOT BE ACHIEVED, AND THAT
105 C THE RETURNED RESULT IS THE BEST WHICH
106 C CAN BE OBTAINED.
107 C = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
108 C SLOWLY CONVERGENT. IT MUST BE NOTED THAT
109 C DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
110 C OF IER.GT.0.
111 C = 6 THE INPUT IS INVALID BECAUSE
112 C NPTS2.LT.2 OR
113 C BREAK POINTS ARE SPECIFIED OUTSIDE
114 C THE INTEGRATION RANGE OR
115 C (EPSABS.LE.0 AND
116 C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
117 C RESULT, ABSERR, NEVAL, LAST ARE SET TO
118 C ZERO. EXEPT WHEN LENIW OR LENW OR NPTS2 IS
119 C INVALID, IWORK(1), IWORK(LIMIT+1),
120 C WORK(LIMIT*2+1) AND WORK(LIMIT*3+1)
121 C ARE SET TO ZERO.
122 C WORK(1) IS SET TO A AND WORK(LIMIT+1)
123 C TO B (WHERE LIMIT = (LENIW-NPTS2)/2).
124 C
125 C DIMENSIONING PARAMETERS
126 C LENIW - INTEGER
127 C DIMENSIONING PARAMETER FOR IWORK
128 C LENIW DETERMINES LIMIT = (LENIW-NPTS2)/2,
129 C WHICH IS THE MAXIMUM NUMBER OF SUBINTERVALS IN THE
130 C PARTITION OF THE GIVEN INTEGRATION INTERVAL (A,B),
131 C LENIW.GE.(3*NPTS2-2).
132 C IF LENIW.LT.(3*NPTS2-2), THE ROUTINE WILL END WITH
133 C IER = 6.
134 C
135 C LENW - INTEGER
136 C DIMENSIONING PARAMETER FOR WORK
137 C LENW MUST BE AT LEAST LENIW*2-NPTS2.
138 C IF LENW.LT.LENIW*2-NPTS2, THE ROUTINE WILL END
139 C WITH IER = 6.
140 C
141 C LAST - INTEGER
142 C ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
143 C PRODUCED IN THE SUBDIVISION PROCESS, WHICH
144 C DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS
145 C ACTUALLY IN THE WORK ARRAYS.
146 C
147 C WORK ARRAYS
148 C IWORK - INTEGER
149 C VECTOR OF DIMENSION AT LEAST LENIW. ON RETURN,
150 C THE FIRST K ELEMENTS OF WHICH CONTAIN
151 C POINTERS TO THE ERROR ESTIMATES OVER THE
152 C SUBINTERVALS, SUCH THAT WORK(LIMIT*3+IWORK(1)),...,
153 C WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
154 C SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
155 C K = LIMIT+1-LAST OTHERWISE
156 C IWORK(LIMIT+1), ...,IWORK(LIMIT+LAST) CONTAIN THE
157 C SUBDIVISION LEVELS OF THE SUBINTERVALS, I.E.
158 C IF (AA,BB) IS A SUBINTERVAL OF (P1,P2)
159 C WHERE P1 AS WELL AS P2 IS A USER-PROVIDED
160 C BREAK POINT OR INTEGRATION LIMIT, THEN (AA,BB) HAS
161 C LEVEL L IF ABS(BB-AA) = ABS(P2-P1)*2**(-L),
162 C IWORK(LIMIT*2+1), ..., IWORK(LIMIT*2+NPTS2) HAVE
163 C NO SIGNIFICANCE FOR THE USER,
164 C NOTE THAT LIMIT = (LENIW-NPTS2)/2.
165 C
166 C WORK - DOUBLE PRECISION
167 C VECTOR OF DIMENSION AT LEAST LENW
168 C ON RETURN
169 C WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
170 C END POINTS OF THE SUBINTERVALS IN THE
171 C PARTITION OF (A,B),
172 C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
173 C THE RIGHT END POINTS,
174 C WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
175 C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
176 C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
177 C CONTAIN THE CORRESPONDING ERROR ESTIMATES,
178 C WORK(LIMIT*4+1), ..., WORK(LIMIT*4+NPTS2)
179 C CONTAIN THE INTEGRATION LIMITS AND THE
180 C BREAK POINTS SORTED IN AN ASCENDING SEQUENCE.
181 C NOTE THAT LIMIT = (LENIW-NPTS2)/2.
182 C
183 C***REFERENCES (NONE)
184 C***ROUTINES CALLED DQAGPE,XERROR
185 C***END PROLOGUE DQAGP
186 C
187  DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,POINTS,RESULT,WORK
188  INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,NEVAL,
189  * npts2
190 C
191  dimension iwork(leniw),points(npts2),work(lenw)
192 C
193  EXTERNAL f
194 C
195 C CHECK VALIDITY OF LIMIT AND LENW.
196 C
197 C***FIRST EXECUTABLE STATEMENT DQAGP
198  ier = 6
199  neval = 0
200  last = 0
201  result = 0.0d+00
202  abserr = 0.0d+00
203  IF(leniw.LT.(3*npts2-2).OR.lenw.LT.(leniw*2-npts2).OR.npts2.LT.2)
204  * GO TO 10
205 C
206 C PREPARE CALL FOR DQAGPE.
207 C
208  limit = (leniw-npts2)/2
209  l1 = limit+1
210  l2 = limit+l1
211  l3 = limit+l2
212  l4 = limit+l3
213 C
214  CALL dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr,
215  * neval,ier,work(1),work(l1),work(l2),work(l3),work(l4),
216  * iwork(1),iwork(l1),iwork(l2),last)
217 C
218 C CALL ERROR HANDLER IF NECESSARY.
219 C
220  lvl = 0
221 10 IF(ier.EQ.6) lvl = 1
222  IF(ier.GT.0) CALL xerror('ABNORMAL RETURN FROM DQAGP',26,ier,lvl)
223  RETURN
224  END
subroutine dqagp(F, A, B, NPTS2, POINTS, EPSABS, EPSREL, RESULT, ABSERR, NEVAL, IER, LENIW, LENW, LAST, IWORK, WORK)
Definition: dqagp.f:3
subroutine dqagpe(F, A, B, NPTS2, POINTS, EPSABS, EPSREL, LIMIT, RESULT, ABSERR, NEVAL, IER, ALIST, BLIST, RLIST, ELIST, PTS, IORD, LEVEL, NDIN, LAST)
Definition: dqagpe.f:4
subroutine xerror(MESSG, NMESSG, NERR, LEVEL)
Definition: xerror.f:2