V3FIT
dcuhre.f
1  SUBROUTINE dcuhre(NDIM,NUMFUN,A,B,MINPTS,MAXPTS,FUNSUB,EPSABS,
2  + EPSREL,KEY,NW,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,
3  + WORK)
4 C***BEGIN PROLOGUE DCUHRE
5 C***DATE WRITTEN 900116 (YYMMDD)
6 C***REVISION DATE 900116 (YYMMDD)
7 C***CATEGORY NO. H2B1A1
8 C***AUTHOR
9 C Jarle Berntsen, The Computing Centre,
10 C University of Bergen, Thormohlens gt. 55,
11 C N-5008 Bergen, Norway
12 C Phone.. 47-5-544055
13 C Email.. jarle@eik.ii.uib.no
14 C Terje O. Espelid, Department of Informatics,
15 C University of Bergen, Thormohlens gt. 55,
16 C N-5008 Bergen, Norway
17 C Phone.. 47-5-544180
18 C Email.. terje@eik.ii.uib.no
19 C Alan Genz, Computer Science Department, Washington State
20 C University, Pullman, WA 99163-2752, USA
21 C Email.. acg@eecs.wsu.edu
22 C***KEYWORDS automatic multidimensional integrator,
23 C n-dimensional hyper-rectangles,
24 C general purpose, global adaptive
25 C***PURPOSE The routine calculates an approximation to a given
26 C vector of definite integrals
27 C
28 C B(1) B(2) B(NDIM)
29 C I I ... I (F ,F ,...,F ) DX(NDIM)...DX(2)DX(1),
30 C A(1) A(2) A(NDIM) 1 2 NUMFUN
31 C
32 C where F = F (X ,X ,...,X ), I = 1,2,...,NUMFUN.
33 C I I 1 2 NDIM
34 C
35 C hopefully satisfying for each component of I the following
36 C claim for accuracy:
37 C ABS(I(K)-RESULT(K)).LE.MAX(EPSABS,EPSREL*ABS(I(K)))
38 C***DESCRIPTION Computation of integrals over hyper-rectangular
39 C regions.
40 C DCUHRE is a driver for the integration routine
41 C DADHRE, which repeatedly subdivides the region
42 C of integration and estimates the integrals and the
43 C errors over the subregions with greatest
44 C estimated errors until the error request
45 C is met or MAXPTS function evaluations have been used.
46 C
47 C For NDIM = 2 the default integration rule is of
48 C degree 13 and uses 65 evaluation points.
49 C For NDIM = 3 the default integration rule is of
50 C degree 11 and uses 127 evaluation points.
51 C For NDIM greater then 3 the default integration rule
52 C is of degree 9 and uses NUM evaluation points where
53 C NUM = 1 + 4*2*NDIM + 2*NDIM*(NDIM-1) + 4*NDIM*(NDIM-1) +
54 C 4*NDIM*(NDIM-1)*(NDIM-2)/3 + 2**NDIM
55 C The degree 9 rule may also be applied for NDIM = 2
56 C and NDIM = 3.
57 C A rule of degree 7 is available in all dimensions.
58 C The number of evaluation
59 C points used by the degree 7 rule is
60 C NUM = 1 + 3*2*NDIM + 2*NDIM*(NDIM-1) + 2**NDIM
61 C
62 C When DCUHRE computes estimates to a vector of
63 C integrals, all components of the vector are given
64 C the same treatment. That is, I(F ) and I(F ) for
65 C J K
66 C J not equal to K, are estimated with the same
67 C subdivision of the region of integration.
68 C For integrals with enough similarity, we may save
69 C time by applying DCUHRE to all integrands in one call.
70 C For integrals that varies continuously as functions of
71 C some parameter, the estimates produced by DCUHRE will
72 C also vary continuously when the same subdivision is
73 C applied to all components. This will generally not be
74 C the case when the different components are given
75 C separate treatment.
76 C
77 C On the other hand this feature should be used with
78 C caution when the different components of the integrals
79 C require clearly different subdivisions.
80 C
81 C ON ENTRY
82 C
83 C NDIM Integer.
84 C Number of variables. 1 < NDIM <= 15.
85 C NUMFUN Integer.
86 C Number of components of the integral.
87 C A Real array of dimension NDIM.
88 C Lower limits of integration.
89 C B Real array of dimension NDIM.
90 C Upper limits of integration.
91 C MINPTS Integer.
92 C Minimum number of function evaluations.
93 C MAXPTS Integer.
94 C Maximum number of function evaluations.
95 C The number of function evaluations over each subregion
96 C is NUM.
97 C If (KEY = 0 or KEY = 1) and (NDIM = 2) Then
98 C NUM = 65
99 C Elseif (KEY = 0 or KEY = 2) and (NDIM = 3) Then
100 C NUM = 127
101 C Elseif (KEY = 0 and NDIM > 3) or (KEY = 3) Then
102 C NUM = 1 + 4*2*NDIM + 2*NDIM*(NDIM-1) + 4*NDIM*(NDIM-1) +
103 C 4*NDIM*(NDIM-1)*(NDIM-2)/3 + 2**NDIM
104 C Elseif (KEY = 4) Then
105 C NUM = 1 + 3*2*NDIM + 2*NDIM*(NDIM-1) + 2**NDIM
106 C MAXPTS >= 3*NUM and MAXPTS >= MINPTS
107 C For 3 < NDIM < 13 the minimum values for MAXPTS are:
108 C NDIM = 4 5 6 7 8 9 10 11 12
109 C KEY = 3: 459 819 1359 2151 3315 5067 7815 12351 20235
110 C KEY = 4: 195 309 483 765 1251 2133 3795 7005 13299
111 C FUNSUB Externally declared subroutine for computing
112 C all components of the integrand at the given
113 C evaluation point.
114 C It must have parameters (NDIM,X,NUMFUN,FUNVLS)
115 C Input parameters:
116 C NDIM Integer that defines the dimension of the
117 C integral.
118 C X Real array of dimension NDIM
119 C that defines the evaluation point.
120 C NUMFUN Integer that defines the number of
121 C components of I.
122 C Output parameter:
123 C FUNVLS Real array of dimension NUMFUN
124 C that defines NUMFUN components of the integrand.
125 C
126 C EPSABS Real.
127 C Requested absolute error.
128 C EPSREL Real.
129 C Requested relative error.
130 C KEY Integer.
131 C Key to selected local integration rule.
132 C KEY = 0 is the default value.
133 C For NDIM = 2 the degree 13 rule is selected.
134 C For NDIM = 3 the degree 11 rule is selected.
135 C For NDIM > 3 the degree 9 rule is selected.
136 C KEY = 1 gives the user the 2 dimensional degree 13
137 C integration rule that uses 65 evaluation points.
138 C KEY = 2 gives the user the 3 dimensional degree 11
139 C integration rule that uses 127 evaluation points.
140 C KEY = 3 gives the user the degree 9 integration rule.
141 C KEY = 4 gives the user the degree 7 integration rule.
142 C This is the recommended rule for problems that
143 C require great adaptivity.
144 C NW Integer.
145 C Defines the length of the working array WORK.
146 C Let MAXSUB denote the maximum allowed number of subregions
147 C for the given values of MAXPTS, KEY and NDIM.
148 C MAXSUB = (MAXPTS-NUM)/(2*NUM) + 1
149 C NW should be greater or equal to
150 C MAXSUB*(2*NDIM+2*NUMFUN+2) + 17*NUMFUN + 1
151 C For efficient execution on parallel computers
152 C NW should be chosen greater or equal to
153 C MAXSUB*(2*NDIM+2*NUMFUN+2) + 17*NUMFUN*MDIV + 1
154 C where MDIV is the number of subregions that are divided in
155 C each subdivision step.
156 C MDIV is default set internally in DCUHRE equal to 1.
157 C For efficient execution on parallel computers
158 C with NPROC processors MDIV should be set equal to
159 C the smallest integer such that MOD(2*MDIV,NPROC) = 0.
160 C
161 C RESTAR Integer.
162 C If RESTAR = 0, this is the first attempt to compute
163 C the integral.
164 C If RESTAR = 1, then we restart a previous attempt.
165 C In this case the only parameters for DCUHRE that may
166 C be changed (with respect to the previous call of DCUHRE)
167 C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR.
168 C
169 C ON RETURN
170 C
171 C RESULT Real array of dimension NUMFUN.
172 C Approximations to all components of the integral.
173 C ABSERR Real array of dimension NUMFUN.
174 C Estimates of absolute errors.
175 C NEVAL Integer.
176 C Number of function evaluations used by DCUHRE.
177 C IFAIL Integer.
178 C IFAIL = 0 for normal exit, when ABSERR(K) <= EPSABS or
179 C ABSERR(K) <= ABS(RESULT(K))*EPSREL with MAXPTS or less
180 C function evaluations for all values of K,
181 C 1 <= K <= NUMFUN .
182 C IFAIL = 1 if MAXPTS was too small for DCUHRE
183 C to obtain the required accuracy. In this case DCUHRE
184 C returns values of RESULT with estimated absolute
185 C errors ABSERR.
186 C IFAIL = 2 if KEY is less than 0 or KEY greater than 4.
187 C IFAIL = 3 if NDIM is less than 2 or NDIM greater than 15.
188 C IFAIL = 4 if KEY = 1 and NDIM not equal to 2.
189 C IFAIL = 5 if KEY = 2 and NDIM not equal to 3.
190 C IFAIL = 6 if NUMFUN is less than 1.
191 C IFAIL = 7 if volume of region of integration is zero.
192 C IFAIL = 8 if MAXPTS is less than 3*NUM.
193 C IFAIL = 9 if MAXPTS is less than MINPTS.
194 C IFAIL = 10 if EPSABS < 0 and EPSREL < 0.
195 C IFAIL = 11 if NW is too small.
196 C IFAIL = 12 if unlegal RESTAR.
197 C WORK Real array of dimension NW.
198 C Used as working storage.
199 C WORK(NW) = NSUB, the number of subregions in the data
200 C structure.
201 C Let WRKSUB=(NW-1-17*NUMFUN*MDIV)/(2*NDIM+2*NUMFUN+2)
202 C WORK(1),...,WORK(NUMFUN*WRKSUB) contain
203 C the estimated components of the integrals over the
204 C subregions.
205 C WORK(NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*WRKSUB) contain
206 C the estimated errors over the subregions.
207 C WORK(2*NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*WRKSUB+NDIM*
208 C WRKSUB) contain the centers of the subregions.
209 C WORK(2*NUMFUN*WRKSUB+NDIM*WRKSUB+1),...,WORK((2*NUMFUN+
210 C NDIM)*WRKSUB+NDIM*WRKSUB) contain subregion half widths.
211 C WORK(2*NUMFUN*WRKSUB+2*NDIM*WRKSUB+1),...,WORK(2*NUMFUN*
212 C WRKSUB+2*NDIM*WRKSUB+WRKSUB) contain the greatest errors
213 C in each subregion.
214 C WORK((2*NUMFUN+2*NDIM+1)*WRKSUB+1),...,WORK((2*NUMFUN+
215 C 2*NDIM+1)*WRKSUB+WRKSUB) contain the direction of
216 C subdivision in each subregion.
217 C WORK(2*(NDIM+NUMFUN+1)*WRKSUB),...,WORK(2*(NDIM+NUMFUN+1)*
218 C WRKSUB+ 17*MDIV*NUMFUN) is used as temporary
219 C storage in DADHRE.
220 C
221 C
222 C DCUHRE Example Test Program
223 C
224 C
225 C DTEST1 is a simple test driver for DCUHRE.
226 C
227 C Output produced on a SUN 3/50.
228 c
229 C DCUHRE TEST RESULTS
230 C
231 C FTEST CALLS = 3549, IFAIL = 0
232 C N ESTIMATED ERROR INTEGRAL
233 C 1 0.000000 0.138508
234 C 2 0.000000 0.063695
235 C 3 0.000009 0.058618
236 C 4 0.000000 0.054070
237 C 5 0.000000 0.050056
238 C 6 0.000000 0.046546
239 C
240 C PROGRAM DTEST1
241 C EXTERNAL FTEST
242 C INTEGER KEY, N, NF, NDIM, MINCLS, MAXCLS, IFAIL, NEVAL, NW
243 C PARAMETER (NDIM = 5, NW = 5000, NF = NDIM+1)
244 C DOUBLE PRECISION A(NDIM), B(NDIM), WRKSTR(NW)
245 C DOUBLE PRECISION ABSEST(NF), FINEST(NF), ABSREQ, RELREQ
246 C DO 10 N = 1,NDIM
247 C A(N) = 0
248 C B(N) = 1
249 C 10 CONTINUE
250 C MINCLS = 0
251 C MAXCLS = 10000
252 C KEY = 0
253 C ABSREQ = 0
254 C RELREQ = 1E-3
255 C CALL DCUHRE(NDIM, NF, A, B, MINCLS, MAXCLS, FTEST, ABSREQ, RELREQ,
256 C * KEY, NW, 0, FINEST, ABSEST, NEVAL, IFAIL, WRKSTR)
257 C PRINT 9999, NEVAL, IFAIL
258 C9999 FORMAT (8X, 'DCUHRE TEST RESULTS', //' FTEST CALLS = ', I4,
259 C * ', IFAIL = ', I2, /' N ESTIMATED ERROR INTEGRAL')
260 C DO 20 N = 1,NF
261 C PRINT 9998, N, ABSEST(N), FINEST(N)
262 C9998 FORMAT (3X, I2, 2F15.6)
263 C 20 CONTINUE
264 C END
265 C SUBROUTINE FTEST(NDIM, Z, NFUN, F)
266 C INTEGER N, NDIM, NFUN
267 C DOUBLE PRECISION Z(NDIM), F(NFUN), SUM
268 C SUM = 0
269 C DO 10 N = 1,NDIM
270 C SUM = SUM + N*Z(N)**2
271 C 10 CONTINUE
272 C F(1) = EXP(-SUM/2)
273 C DO 20 N = 1,NDIM
274 C F(N+1) = Z(N)*F(1)
275 C 20 CONTINUE
276 C END
277 C
278 C***LONG DESCRIPTION
279 C
280 C The information for each subregion is contained in the
281 C data structure WORK.
282 C When passed on to DADHRE, WORK is split into eight
283 C arrays VALUES, ERRORS, CENTRS, HWIDTS, GREATE, DIR,
284 C OLDRES and WORK.
285 C VALUES contains the estimated values of the integrals.
286 C ERRORS contains the estimated errors.
287 C CENTRS contains the centers of the subregions.
288 C HWIDTS contains the half widths of the subregions.
289 C GREATE contains the greatest estimated error for each subregion.
290 C DIR contains the directions for further subdivision.
291 C OLDRES and WORK are used as work arrays in DADHRE.
292 C
293 C The data structures for the subregions are in DADHRE organized
294 C as a heap, and the size of GREATE(I) defines the position of
295 C region I in the heap. The heap is maintained by the program
296 C DTRHRE.
297 C
298 C The subroutine DADHRE is written for efficient execution on shared
299 C memory parallel computer. On a computer with NPROC processors we will
300 C in each subdivision step divide MDIV regions, where MDIV is
301 C chosen such that MOD(2*MDIV,NPROC) = 0, in totally 2*MDIV new regions.
302 C Each processor will then compute estimates of the integrals and errors
303 C over 2*MDIV/NPROC subregions in each subdivision step.
304 C The subroutine for estimating the integral and the error over
305 C each subregion, DRLHRE, uses WORK2 as a work array.
306 C We must make sure that each processor writes its results to
307 C separate parts of the memory, and therefore the sizes of WORK and
308 C WORK2 are functions of MDIV.
309 C In order to achieve parallel processing of subregions, compiler
310 C directives should be placed in front of the DO 200
311 C loop in DADHRE on machines like Alliant and CRAY.
312 C
313 C***REFERENCES
314 C J.Berntsen, T.O.Espelid and A.Genz, An Adaptive Algorithm
315 C for the Approximate Calculation of Multiple Integrals,
316 C To be published.
317 C
318 C J.Berntsen, T.O.Espelid and A.Genz, DCUHRE: An Adaptive
319 C Multidimensional Integration Routine for a Vector of
320 C Integrals, To be published.
321 C
322 C***ROUTINES CALLED DCHHRE,DADHRE
323 C***END PROLOGUE DCUHRE
324 C
325 C Global variables.
326 C
327  EXTERNAL funsub
328  INTEGER NDIM,NUMFUN,MINPTS,MAXPTS,KEY,NW,RESTAR
329  INTEGER NEVAL,IFAIL
330  DOUBLE PRECISION A(NDIM),B(NDIM),EPSABS,EPSREL
331  DOUBLE PRECISION RESULT(NUMFUN),ABSERR(NUMFUN),WORK(NW)
332 C
333 C Local variables.
334 C
335 C MDIV Integer.
336 C MDIV is the number of subregions that are divided in
337 C each subdivision step in DADHRE.
338 C MDIV is chosen default to 1.
339 C For efficient execution on parallel computers
340 C with NPROC processors MDIV should be set equal to
341 C the smallest integer such that MOD(2*MDIV,NPROC) = 0.
342 C MAXDIM Integer.
343 C The maximum allowed value of NDIM.
344 C MAXWT Integer. The maximum number of weights used by the
345 C integration rule.
346 C WTLENG Integer.
347 C The number of generators used by the selected rule.
348 C WORK2 Real work space. The length
349 C depends on the parameters MDIV,MAXDIM and MAXWT.
350 C MAXSUB Integer.
351 C The maximum allowed number of subdivisions
352 C for the given values of KEY, NDIM and MAXPTS.
353 C MINSUB Integer.
354 C The minimum allowed number of subregions for the given
355 C values of MINPTS, KEY and NDIM.
356 C WRKSUB Integer.
357 C The maximum allowed number of subregions as a function
358 C of NW, NUMFUN, NDIM and MDIV. This determines the length
359 C of the main work arrays.
360 C NUM Integer. The number of integrand evaluations needed
361 C over each subregion.
362 C
363  INTEGER MDIV,MAXWT,WTLENG,MAXDIM,LENW2,MAXSUB,MINSUB
364  INTEGER NUM,NSUB,LENW,KEYF
365  parameter(mdiv=1)
366  parameter(maxdim=15)
367  parameter(maxwt=14)
368  parameter(lenw2=2*mdiv*maxdim* (maxwt+1)+12*maxwt+2*maxdim)
369  INTEGER WRKSUB,I1,I2,I3,I4,I5,I6,I7,I8,K1,K2,K3,K4,K5,K6,K7,K8
370  DOUBLE PRECISION WORK2(LENW2)
371 C
372 C***FIRST EXECUTABLE STATEMENT DCUHRE
373 C
374 C Compute NUM, WTLENG, MAXSUB and MINSUB,
375 C and check the input parameters.
376 C
377  CALL dchhre(maxdim,ndim,numfun,mdiv,a,b,minpts,maxpts,epsabs,
378  + epsrel,key,nw,restar,num,maxsub,minsub,keyf,
379  + ifail,wtleng)
380  wrksub = (nw - 1 - 17*mdiv*numfun)/(2*ndim + 2*numfun + 2)
381  IF (ifail.NE.0) THEN
382  GO TO 999
383  END IF
384 C
385 C Split up the work space.
386 C
387  i1 = 1
388  i2 = i1 + wrksub*numfun
389  i3 = i2 + wrksub*numfun
390  i4 = i3 + wrksub*ndim
391  i5 = i4 + wrksub*ndim
392  i6 = i5 + wrksub
393  i7 = i6 + wrksub
394  i8 = i7 + numfun*mdiv
395  k1 = 1
396  k2 = k1 + 2*mdiv*wtleng*ndim
397  k3 = k2 + wtleng*5
398  k4 = k3 + wtleng
399  k5 = k4 + ndim
400  k6 = k5 + ndim
401  k7 = k6 + 2*mdiv*ndim
402  k8 = k7 + 3*wtleng
403 C
404 C On restart runs the number of subregions from the
405 C previous call is assigned to NSUB.
406 C
407  IF (restar.EQ.1) THEN
408  nsub = work(nw)
409  END IF
410 C
411 C Compute the size of the temporary work space needed in DADHRE.
412 C
413  lenw = 16*mdiv*numfun
414  CALL dadhre(ndim,numfun,mdiv,a,b,minsub,maxsub,funsub,epsabs,
415  + epsrel,keyf,restar,num,lenw,wtleng,
416  + result,abserr,neval,nsub,ifail,work(i1),work(i2),
417  + work(i3),work(i4),work(i5),work(i6),work(i7),work(i8),
418  + work2(k1),work2(k2),work2(k3),work2(k4),work2(k5),
419  + work2(k6),work2(k7),work2(k8))
420  work(nw) = nsub
421 999 RETURN
422 C
423 C***END DCUHRE
424 C
425  END
426  SUBROUTINE dchhre(MAXDIM,NDIM,NUMFUN,MDIV,A,B,MINPTS,MAXPTS,
427  + EPSABS,EPSREL,KEY,NW,RESTAR,NUM,MAXSUB,MINSUB,
428  + KEYF,IFAIL,WTLENG)
429 C***BEGIN PROLOGUE DCHHRE
430 C***PURPOSE DCHHRE checks the validity of the
431 C input parameters to DCUHRE.
432 C***DESCRIPTION
433 C DCHHRE computes NUM, MAXSUB, MINSUB, KEYF, WTLENG and
434 C IFAIL as functions of the input parameters to DCUHRE,
435 C and checks the validity of the input parameters to DCUHRE.
436 C
437 C ON ENTRY
438 C
439 C MAXDIM Integer.
440 C The maximum allowed number of dimensions.
441 C NDIM Integer.
442 C Number of variables. 1 < NDIM <= MAXDIM.
443 C NUMFUN Integer.
444 C Number of components of the integral.
445 C MDIV Integer.
446 C MDIV is the number of subregions that are divided in
447 C each subdivision step in DADHRE.
448 C MDIV is chosen default to 1.
449 C For efficient execution on parallel computers
450 C with NPROC processors MDIV should be set equal to
451 C the smallest integer such that MOD(2*MDIV,NPROC) = 0.
452 C A Real array of dimension NDIM.
453 C Lower limits of integration.
454 C B Real array of dimension NDIM.
455 C Upper limits of integration.
456 C MINPTS Integer.
457 C Minimum number of function evaluations.
458 C MAXPTS Integer.
459 C Maximum number of function evaluations.
460 C The number of function evaluations over each subregion
461 C is NUM.
462 C If (KEY = 0 or KEY = 1) and (NDIM = 2) Then
463 C NUM = 65
464 C Elseif (KEY = 0 or KEY = 2) and (NDIM = 3) Then
465 C NUM = 127
466 C Elseif (KEY = 0 and NDIM > 3) or (KEY = 3) Then
467 C NUM = 1 + 4*2*NDIM + 2*NDIM*(NDIM-1) + 4*NDIM*(NDIM-1) +
468 C 4*NDIM*(NDIM-1)*(NDIM-2)/3 + 2**NDIM
469 C Elseif (KEY = 4) Then
470 C NUM = 1 + 3*2*NDIM + 2*NDIM*(NDIM-1) + 2**NDIM
471 C MAXPTS >= 3*NUM and MAXPTS >= MINPTS
472 C EPSABS Real.
473 C Requested absolute error.
474 C EPSREL Real.
475 C Requested relative error.
476 C KEY Integer.
477 C Key to selected local integration rule.
478 C KEY = 0 is the default value.
479 C For NDIM = 2 the degree 13 rule is selected.
480 C For NDIM = 3 the degree 11 rule is selected.
481 C For NDIM > 3 the degree 9 rule is selected.
482 C KEY = 1 gives the user the 2 dimensional degree 13
483 C integration rule that uses 65 evaluation points.
484 C KEY = 2 gives the user the 3 dimensional degree 11
485 C integration rule that uses 127 evaluation points.
486 C KEY = 3 gives the user the degree 9 integration rule.
487 C KEY = 4 gives the user the degree 7 integration rule.
488 C This is the recommended rule for problems that
489 C require great adaptivity.
490 C NW Integer.
491 C Defines the length of the working array WORK.
492 C Let MAXSUB denote the maximum allowed number of subregions
493 C for the given values of MAXPTS, KEY and NDIM.
494 C MAXSUB = (MAXPTS-NUM)/(2*NUM) + 1
495 C NW should be greater or equal to
496 C MAXSUB*(2*NDIM+2*NUMFUN+2) + 17*NUMFUN + 1
497 C For efficient execution on parallel computers
498 C NW should be chosen greater or equal to
499 C MAXSUB*(2*NDIM+2*NUMFUN+2) + 17*NUMFUN*MDIV + 1
500 C where MDIV is the number of subregions that are divided in
501 C each subdivision step.
502 C MDIV is default set internally in DCUHRE equal to 1.
503 C For efficient execution on parallel computers
504 C with NPROC processors MDIV should be set equal to
505 C the smallest integer such that MOD(2*MDIV,NPROC) = 0.
506 C RESTAR Integer.
507 C If RESTAR = 0, this is the first attempt to compute
508 C the integral.
509 C If RESTAR = 1, then we restart a previous attempt.
510 C
511 C ON RETURN
512 C
513 C NUM Integer.
514 C The number of function evaluations over each subregion.
515 C MAXSUB Integer.
516 C The maximum allowed number of subregions for the
517 C given values of MAXPTS, KEY and NDIM.
518 C MINSUB Integer.
519 C The minimum allowed number of subregions for the given
520 C values of MINPTS, KEY and NDIM.
521 C IFAIL Integer.
522 C IFAIL = 0 for normal exit.
523 C IFAIL = 2 if KEY is less than 0 or KEY greater than 4.
524 C IFAIL = 3 if NDIM is less than 2 or NDIM greater than
525 C MAXDIM.
526 C IFAIL = 4 if KEY = 1 and NDIM not equal to 2.
527 C IFAIL = 5 if KEY = 2 and NDIM not equal to 3.
528 C IFAIL = 6 if NUMFUN less than 1.
529 C IFAIL = 7 if volume of region of integration is zero.
530 C IFAIL = 8 if MAXPTS is less than 3*NUM.
531 C IFAIL = 9 if MAXPTS is less than MINPTS.
532 C IFAIL = 10 if EPSABS < 0 and EPSREL < 0.
533 C IFAIL = 11 if NW is too small.
534 C IFAIL = 12 if unlegal RESTAR.
535 C KEYF Integer.
536 C Key to selected integration rule.
537 C WTLENG Integer.
538 C The number of generators of the chosen integration rule.
539 C
540 C***ROUTINES CALLED-NONE
541 C***END PROLOGUE DCHHRE
542 C
543 C Global variables.
544 C
545  INTEGER NDIM,NUMFUN,MDIV,MINPTS,MAXPTS,KEY,NW,MINSUB,MAXSUB
546  INTEGER RESTAR,NUM,KEYF,IFAIL,MAXDIM,WTLENG
547  DOUBLE PRECISION A(NDIM),B(NDIM),EPSABS,EPSREL
548 C
549 C Local variables.
550 C
551  INTEGER LIMIT,J
552 C
553 C***FIRST EXECUTABLE STATEMENT DCHHRE
554 C
555  ifail = 0
556 C
557 C Check on legal KEY.
558 C
559  IF (key.LT.0 .OR. key.GT.4) THEN
560  ifail = 2
561  GO TO 999
562  END IF
563 C
564 C Check on legal NDIM.
565 C
566  IF (ndim.LT.2 .OR. ndim.GT.maxdim) THEN
567  ifail = 3
568  GO TO 999
569  END IF
570 C
571 C For KEY = 1, NDIM must be equal to 2.
572 C
573  IF (key.EQ.1 .AND. ndim.NE.2) THEN
574  ifail = 4
575  GO TO 999
576  END IF
577 C
578 C For KEY = 2, NDIM must be equal to 3.
579 C
580  IF (key.EQ.2 .AND. ndim.NE.3) THEN
581  ifail = 5
582  GO TO 999
583  END IF
584 C
585 C For KEY = 0, we point at the selected integration rule.
586 C
587  IF (key.EQ.0) THEN
588  IF (ndim.EQ.2) THEN
589  keyf = 1
590  ELSE IF (ndim.EQ.3) THEN
591  keyf = 2
592  ELSE
593  keyf = 3
594  ENDIF
595  ELSE
596  keyf = key
597  ENDIF
598 C
599 C Compute NUM and WTLENG as a function of KEYF and NDIM.
600 C
601  IF (keyf.EQ.1) THEN
602  num = 65
603  wtleng = 14
604  ELSE IF (keyf.EQ.2) THEN
605  num = 127
606  wtleng = 13
607  ELSE IF (keyf.EQ.3) THEN
608  num = 1 + 4*2*ndim + 2*ndim* (ndim-1) + 4*ndim* (ndim-1) +
609  + 4*ndim* (ndim-1)* (ndim-2)/3 + 2**ndim
610  wtleng = 9
611  IF (ndim.EQ.2) wtleng = 8
612  ELSE IF (keyf.EQ.4) THEN
613  num = 1 + 3*2*ndim + 2*ndim* (ndim-1) + 2**ndim
614  wtleng = 6
615  END IF
616 C
617 C Compute MAXSUB.
618 C
619  maxsub = (maxpts-num)/ (2*num) + 1
620 C
621 C Compute MINSUB.
622 C
623  minsub = (minpts-num)/ (2*num) + 1
624  IF (mod(minpts-num,2*num).NE.0) THEN
625  minsub = minsub + 1
626  END IF
627  minsub = max(2,minsub)
628 C
629 C Check on positive NUMFUN.
630 C
631  IF (numfun.LT.1) THEN
632  ifail = 6
633  GO TO 999
634  END IF
635 C
636 C Check on legal upper and lower limits of integration.
637 C
638  DO 10 j = 1,ndim
639  IF (a(j)-b(j).EQ.0) THEN
640  ifail = 7
641  GO TO 999
642  END IF
643 10 CONTINUE
644 C
645 C Check on MAXPTS < 3*NUM.
646 C
647  IF (maxpts.LT.3*num) THEN
648  ifail = 8
649  GO TO 999
650  END IF
651 C
652 C Check on MAXPTS >= MINPTS.
653 C
654  IF (maxpts.LT.minpts) THEN
655  ifail = 9
656  GO TO 999
657  END IF
658 C
659 C Check on legal accuracy requests.
660 C
661  IF (epsabs.LT.0 .AND. epsrel.LT.0) THEN
662  ifail = 10
663  GO TO 999
664  END IF
665 C
666 C Check on big enough double precision workspace.
667 C
668  limit = maxsub* (2*ndim+2*numfun+2) + 17*mdiv*numfun + 1
669  IF (nw.LT.limit) THEN
670  ifail = 11
671  GO TO 999
672  END IF
673 C
674 C Check on legal RESTAR.
675 C
676  IF (restar.NE.0 .AND. restar.NE.1) THEN
677  ifail = 12
678  GO TO 999
679  END IF
680 999 RETURN
681 C
682 C***END DCHHRE
683 C
684  END
685  SUBROUTINE dadhre(NDIM,NUMFUN,MDIV,A,B,MINSUB,MAXSUB,FUNSUB,
686  + EPSABS,EPSREL,KEY,RESTAR,NUM,LENW,WTLENG,
687  + RESULT,ABSERR,NEVAL,NSUB,IFAIL,VALUES,
688  + ERRORS,CENTRS,HWIDTS,GREATE,DIR,OLDRES,WORK,G,W,
689  + RULPTS,CENTER,HWIDTH,X,SCALES,NORMS)
690 C***BEGIN PROLOGUE DADHRE
691 C***KEYWORDS automatic multidimensional integrator,
692 C n-dimensional hyper-rectangles,
693 C general purpose, global adaptive
694 C***PURPOSE The routine calculates an approximation to a given
695 C vector of definite integrals, I, over a hyper-rectangular
696 C region hopefully satisfying for each component of I the
697 C following claim for accuracy:
698 C ABS(I(K)-RESULT(K)).LE.MAX(EPSABS,EPSREL*ABS(I(K)))
699 C***DESCRIPTION Computation of integrals over hyper-rectangular
700 C regions.
701 C DADHRE repeatedly subdivides the region
702 C of integration and estimates the integrals and the
703 C errors over the subregions with greatest
704 C estimated errors until the error request
705 C is met or MAXSUB subregions are stored.
706 C The regions are divided in two equally sized parts along
707 C the direction with greatest absolute fourth divided
708 C difference.
709 C
710 C ON ENTRY
711 C
712 C NDIM Integer.
713 C Number of variables. 1 < NDIM <= MAXDIM.
714 C NUMFUN Integer.
715 C Number of components of the integral.
716 C MDIV Integer.
717 C Defines the number of new subregions that are divided
718 C in each subdivision step.
719 C A Real array of dimension NDIM.
720 C Lower limits of integration.
721 C B Real array of dimension NDIM.
722 C Upper limits of integration.
723 C MINSUB Integer.
724 C The computations proceed until there are at least
725 C MINSUB subregions in the data structure.
726 C MAXSUB Integer.
727 C The computations proceed until there are at most
728 C MAXSUB subregions in the data structure.
729 C
730 C FUNSUB Externally declared subroutine for computing
731 C all components of the integrand in the given
732 C evaluation point.
733 C It must have parameters (NDIM,X,NUMFUN,FUNVLS)
734 C Input parameters:
735 C NDIM Integer that defines the dimension of the
736 C integral.
737 C X Real array of dimension NDIM
738 C that defines the evaluation point.
739 C NUMFUN Integer that defines the number of
740 C components of I.
741 C Output parameter:
742 C FUNVLS Real array of dimension NUMFUN
743 C that defines NUMFUN components of the integrand.
744 C
745 C EPSABS Real.
746 C Requested absolute error.
747 C EPSREL Real.
748 C Requested relative error.
749 C KEY Integer.
750 C Key to selected local integration rule.
751 C KEY = 0 is the default value.
752 C For NDIM = 2 the degree 13 rule is selected.
753 C For NDIM = 3 the degree 11 rule is selected.
754 C For NDIM > 3 the degree 9 rule is selected.
755 C KEY = 1 gives the user the 2 dimensional degree 13
756 C integration rule that uses 65 evaluation points.
757 C KEY = 2 gives the user the 3 dimensional degree 11
758 C integration rule that uses 127 evaluation points.
759 C KEY = 3 gives the user the degree 9 integration rule.
760 C KEY = 4 gives the user the degree 7 integration rule.
761 C This is the recommended rule for problems that
762 C require great adaptivity.
763 C RESTAR Integer.
764 C If RESTAR = 0, this is the first attempt to compute
765 C the integral.
766 C If RESTAR = 1, then we restart a previous attempt.
767 C (In this case the output parameters from DADHRE
768 C must not be changed since the last
769 C exit from DADHRE.)
770 C NUM Integer.
771 C The number of function evaluations over each subregion.
772 C LENW Integer.
773 C Defines the length of the working array WORK.
774 C LENW should be greater or equal to
775 C 16*MDIV*NUMFUN.
776 C WTLENG Integer.
777 C The number of weights in the basic integration rule.
778 C NSUB Integer.
779 C If RESTAR = 1, then NSUB must specify the number
780 C of subregions stored in the previous call to DADHRE.
781 C
782 C ON RETURN
783 C
784 C RESULT Real array of dimension NUMFUN.
785 C Approximations to all components of the integral.
786 C ABSERR Real array of dimension NUMFUN.
787 C Estimates of absolute accuracies.
788 C NEVAL Integer.
789 C Number of function evaluations used by DADHRE.
790 C NSUB Integer.
791 C Number of stored subregions.
792 C IFAIL Integer.
793 C IFAIL = 0 for normal exit, when ABSERR(K) <= EPSABS or
794 C ABSERR(K) <= ABS(RESULT(K))*EPSREL with MAXSUB or less
795 C subregions processed for all values of K,
796 C 1 <= K <= NUMFUN.
797 C IFAIL = 1 if MAXSUB was too small for DADHRE
798 C to obtain the required accuracy. In this case DADHRE
799 C returns values of RESULT with estimated absolute
800 C accuracies ABSERR.
801 C VALUES Real array of dimension (NUMFUN,MAXSUB).
802 C Used to store estimated values of the integrals
803 C over the subregions.
804 C ERRORS Real array of dimension (NUMFUN,MAXSUB).
805 C Used to store the corresponding estimated errors.
806 C CENTRS Real array of dimension (NDIM,MAXSUB).
807 C Used to store the centers of the stored subregions.
808 C HWIDTS Real array of dimension (NDIM,MAXSUB).
809 C Used to store the half widths of the stored subregions.
810 C GREATE Real array of dimension MAXSUB.
811 C Used to store the greatest estimated errors in
812 C all subregions.
813 C DIR Real array of dimension MAXSUB.
814 C DIR is used to store the directions for
815 C further subdivision.
816 C OLDRES Real array of dimension (NUMFUN,MDIV).
817 C Used to store old estimates of the integrals over the
818 C subregions.
819 C WORK Real array of dimension LENW.
820 C Used in DRLHRE and DTRHRE.
821 C G Real array of dimension (NDIM,WTLENG,2*MDIV).
822 C The fully symmetric sum generators for the rules.
823 C G(1,J,1),...,G(NDIM,J,1) are the generators for the
824 C points associated with the Jth weights.
825 C When MDIV subregions are divided in 2*MDIV
826 C subregions, the subregions may be processed on different
827 C processors and we must make a copy of the generators
828 C for each processor.
829 C W Real array of dimension (5,WTLENG).
830 C The weights for the basic and null rules.
831 C W(1,1), ..., W(1,WTLENG) are weights for the basic rule.
832 C W(I,1), ..., W(I,WTLENG) , for I > 1 are null rule weights.
833 C RULPTS Real array of dimension WTLENG.
834 C Work array used in DINHRE.
835 C CENTER Real array of dimension NDIM.
836 C Work array used in DTRHRE.
837 C HWIDTH Real array of dimension NDIM.
838 C Work array used in DTRHRE.
839 C X Real array of dimension (NDIM,2*MDIV).
840 C Work array used in DRLHRE.
841 C SCALES Real array of dimension (3,WTLENG).
842 C Work array used by DINHRE and DRLHRE.
843 C NORMS Real array of dimension (3,WTLENG).
844 C Work array used by DINHRE and DRLHRE.
845 C
846 C***REFERENCES
847 C
848 C P. van Dooren and L. de Ridder, Algorithm 6, An adaptive algorithm
849 C for numerical integration over an n-dimensional cube, J.Comput.Appl.
850 C Math. 2(1976)207-217.
851 C
852 C A.C.Genz and A.A.Malik, Algorithm 019. Remarks on algorithm 006:
853 C An adaptive algorithm for numerical integration over an
854 C N-dimensional rectangular region,J.Comput.Appl.Math. 6(1980)295-302.
855 C
856 C***ROUTINES CALLED DTRHRE,DINHRE,DRLHRE
857 C***END PROLOGUE DADHRE
858 C
859 C Global variables.
860 C
861  EXTERNAL funsub
862  INTEGER NDIM,NUMFUN,MDIV,MINSUB,MAXSUB,KEY,LENW,RESTAR
863  INTEGER NUM,NEVAL,NSUB,IFAIL,WTLENG
864  DOUBLE PRECISION A(NDIM),B(NDIM),EPSABS,EPSREL
865  DOUBLE PRECISION RESULT(NUMFUN),ABSERR(NUMFUN)
866  DOUBLE PRECISION VALUES(NUMFUN,MAXSUB),ERRORS(NUMFUN,MAXSUB)
867  DOUBLE PRECISION CENTRS(NDIM,MAXSUB)
868  DOUBLE PRECISION HWIDTS(NDIM,MAXSUB)
869  DOUBLE PRECISION GREATE(MAXSUB),DIR(MAXSUB)
870  DOUBLE PRECISION OLDRES(NUMFUN,MDIV)
871  DOUBLE PRECISION WORK(LENW),RULPTS(WTLENG)
872  DOUBLE PRECISION G(NDIM,WTLENG,2*MDIV),W(5,WTLENG)
873  DOUBLE PRECISION CENTER(NDIM),HWIDTH(NDIM),X(NDIM,2*MDIV)
874  DOUBLE PRECISION SCALES(3,WTLENG),NORMS(3,WTLENG)
875 C
876 C Local variables.
877 C
878 C INTSGN is used to get correct sign on the integral.
879 C SBRGNS is the number of stored subregions.
880 C NDIV The number of subregions to be divided in each main step.
881 C POINTR Pointer to the position in the data structure where
882 C the new subregions are to be stored.
883 C DIRECT Direction of subdivision.
884 C ERRCOF Heuristic error coeff. defined in DINHRE and used by DRLHRE
885 C and DADHRE.
886 C
887  INTEGER I,J,K
888  INTEGER INTSGN,SBRGNS
889  INTEGER L1
890  INTEGER NDIV,POINTR,DIRECT,INDEX
891  DOUBLE PRECISION OLDCEN,EST1,EST2,ERRCOF(6)
892 C
893 C***FIRST EXECUTABLE STATEMENT DADHRE
894 C
895 C Get the correct sign on the integral.
896 C
897  intsgn = 1
898  DO 10 j = 1,ndim
899  IF (b(j).LT.a(j)) THEN
900  intsgn = - intsgn
901  END IF
902 10 CONTINUE
903 C
904 C Call DINHRE to compute the weights and abscissas of
905 C the function evaluation points.
906 C
907  CALL dinhre(ndim,key,wtleng,w,g,errcof,rulpts,scales,norms)
908 C
909 C If RESTAR = 1, then this is a restart run.
910 C
911  IF (restar.EQ.1) THEN
912  sbrgns = nsub
913  GO TO 110
914  END IF
915 C
916 C Initialize the SBRGNS, CENTRS and HWIDTS.
917 C
918  sbrgns = 1
919  DO 15 j = 1,ndim
920  centrs(j,1) = (a(j)+b(j))/2
921  hwidts(j,1) = abs(b(j)-a(j))/2
922 15 CONTINUE
923 C
924 C Initialize RESULT, ABSERR and NEVAL.
925 C
926  DO 20 j = 1,numfun
927  result(j) = 0
928  abserr(j) = 0
929 20 CONTINUE
930  neval = 0
931 C
932 C Apply DRLHRE over the whole region.
933 C
934  CALL drlhre(ndim,centrs(1,1),hwidts(1,1),wtleng,g,w,errcof,numfun,
935  + funsub,scales,norms,x,work,values(1,1),errors(1,1),
936  + dir(1))
937  neval = neval + num
938 C
939 C Add the computed values to RESULT and ABSERR.
940 C
941  DO 55 j = 1,numfun
942  result(j) = result(j) + values(j,1)
943 55 CONTINUE
944  DO 65 j = 1,numfun
945  abserr(j) = abserr(j) + errors(j,1)
946 65 CONTINUE
947 C
948 C Store results in heap.
949 C
950  index = 1
951  CALL dtrhre(2,ndim,numfun,index,values,errors,centrs,hwidts,
952  + greate,work(1),work(numfun+1),center,hwidth,dir)
953 C
954 C***End initialisation.
955 C
956 C***Begin loop while the error is too great,
957 C and SBRGNS+1 is less than MAXSUB.
958 C
959 110 IF (sbrgns+1.LE.maxsub) THEN
960 C
961 C If we are allowed to divide further,
962 C prepare to apply basic rule over each half of the
963 C NDIV subregions with greatest errors.
964 C If MAXSUB is great enough, NDIV = MDIV
965 C
966  IF (mdiv.GT.1) THEN
967  ndiv = maxsub - sbrgns
968  ndiv = min(ndiv,mdiv,sbrgns)
969  ELSE
970  ndiv = 1
971  END IF
972 C
973 C Divide the NDIV subregions in two halves, and compute
974 C integral and error over each half.
975 C
976  DO 150 i = 1,ndiv
977  pointr = sbrgns + ndiv + 1 - i
978 C
979 C Adjust RESULT and ABSERR.
980 C
981  DO 115 j = 1,numfun
982  result(j) = result(j) - values(j,1)
983  abserr(j) = abserr(j) - errors(j,1)
984 115 CONTINUE
985 C
986 C Compute first half region.
987 C
988  DO 120 j = 1,ndim
989  centrs(j,pointr) = centrs(j,1)
990  hwidts(j,pointr) = hwidts(j,1)
991 120 CONTINUE
992  direct = dir(1)
993  dir(pointr) = direct
994  hwidts(direct,pointr) = hwidts(direct,1)/2
995  oldcen = centrs(direct,1)
996  centrs(direct,pointr) = oldcen - hwidts(direct,pointr)
997 C
998 C Save the computed values of the integrals.
999 C
1000  DO 125 j = 1,numfun
1001  oldres(j,ndiv-i+1) = values(j,1)
1002 125 CONTINUE
1003 C
1004 C Adjust the heap.
1005 C
1006  CALL dtrhre(1,ndim,numfun,sbrgns,values,errors,centrs,
1007  + hwidts,greate,work(1),work(numfun+1),center,
1008  + hwidth,dir)
1009 C
1010 C Compute second half region.
1011 C
1012  DO 130 j = 1,ndim
1013  centrs(j,pointr-1) = centrs(j,pointr)
1014  hwidts(j,pointr-1) = hwidts(j,pointr)
1015 130 CONTINUE
1016  centrs(direct,pointr-1) = oldcen + hwidts(direct,pointr)
1017  hwidts(direct,pointr-1) = hwidts(direct,pointr)
1018  dir(pointr-1) = direct
1019 150 CONTINUE
1020 C
1021 C Make copies of the generators for each processor.
1022 C
1023  DO 190 i = 2,2*ndiv
1024  DO 190 j = 1,ndim
1025  DO 190 k = 1,wtleng
1026  g(j,k,i) = g(j,k,1)
1027 190 CONTINUE
1028 C
1029 C Apply basic rule.
1030 C
1031 Cvd$l cncall
1032  DO 200 i = 1,2*ndiv
1033  index = sbrgns + i
1034  l1 = 1 + (i-1)*8*numfun
1035  CALL drlhre(ndim,centrs(1,index),hwidts(1,index),wtleng,
1036  + g(1,1,i),w,errcof,numfun,funsub,scales,norms,
1037  + x(1,i),work(l1),values(1,index),
1038  + errors(1,index),dir(index))
1039 200 CONTINUE
1040  neval = neval + 2*ndiv*num
1041 C
1042 C Add new contributions to RESULT.
1043 C
1044  DO 220 i = 1,2*ndiv
1045  DO 210 j = 1,numfun
1046  result(j) = result(j) + values(j,sbrgns+i)
1047 210 CONTINUE
1048 220 CONTINUE
1049 C
1050 C Check consistency of results and if necessary adjust
1051 C the estimated errors.
1052 C
1053  DO 240 i = 1,ndiv
1054  greate(sbrgns+2*i-1) = 0
1055  greate(sbrgns+2*i) = 0
1056  DO 230 j = 1,numfun
1057  est1 = abs(oldres(j,i)- (values(j,
1058  + sbrgns+2*i-1)+values(j,sbrgns+2*i)))
1059  est2 = errors(j,sbrgns+2*i-1) + errors(j,sbrgns+2*i)
1060  IF (est2.GT.0) THEN
1061  errors(j,sbrgns+2*i-1) = errors(j,sbrgns+2*i-1)*
1062  + (1+errcof(5)*est1/est2)
1063  errors(j,sbrgns+2*i) = errors(j,sbrgns+2*i)*
1064  + (1+errcof(5)*est1/est2)
1065  END IF
1066  errors(j,sbrgns+2*i-1) = errors(j,sbrgns+2*i-1) +
1067  + errcof(6)*est1
1068  errors(j,sbrgns+2*i) = errors(j,sbrgns+2*i) +
1069  + errcof(6)*est1
1070  IF (errors(j,sbrgns+2*i-1).GT.
1071  + greate(sbrgns+2*i-1)) THEN
1072  greate(sbrgns+2*i-1) = errors(j,sbrgns+2*i-1)
1073  END IF
1074  IF (errors(j,sbrgns+2*i).GT.greate(sbrgns+2*i)) THEN
1075  greate(sbrgns+2*i) = errors(j,sbrgns+2*i)
1076  END IF
1077  abserr(j) = abserr(j) + errors(j,sbrgns+2*i-1) +
1078  + errors(j,sbrgns+2*i)
1079 230 CONTINUE
1080 240 CONTINUE
1081 C
1082 C Store results in heap.
1083 C
1084  DO 250 i = 1,2*ndiv
1085  index = sbrgns + i
1086  CALL dtrhre(2,ndim,numfun,index,values,errors,centrs,
1087  + hwidts,greate,work(1),work(numfun+1),center,
1088  + hwidth,dir)
1089 250 CONTINUE
1090  sbrgns = sbrgns + 2*ndiv
1091 C
1092 C Check for termination.
1093 C
1094  IF (sbrgns.LT.minsub) THEN
1095  GO TO 110
1096  END IF
1097  DO 255 j = 1,numfun
1098  IF (abserr(j).GT.epsrel*abs(result(j)) .AND.
1099  + abserr(j).GT.epsabs) THEN
1100  GO TO 110
1101  END IF
1102 255 CONTINUE
1103  ifail = 0
1104  GO TO 499
1105 C
1106 C Else we did not succeed with the
1107 C given value of MAXSUB.
1108 C
1109  ELSE
1110  ifail = 1
1111  END IF
1112 C
1113 C Compute more accurate values of RESULT and ABSERR.
1114 C
1115 499 CONTINUE
1116  DO 500 j = 1,numfun
1117  result(j) = 0
1118  abserr(j) = 0
1119 500 CONTINUE
1120  DO 510 i = 1,sbrgns
1121  DO 505 j = 1,numfun
1122  result(j) = result(j) + values(j,i)
1123  abserr(j) = abserr(j) + errors(j,i)
1124 505 CONTINUE
1125 510 CONTINUE
1126 C
1127 C Compute correct sign on the integral.
1128 C
1129  DO 600 j = 1,numfun
1130  result(j) = result(j)*intsgn
1131 600 CONTINUE
1132  nsub = sbrgns
1133  RETURN
1134 C
1135 C***END DADHRE
1136 C
1137  END
1138  SUBROUTINE dinhre(NDIM,KEY,WTLENG,W,G,ERRCOF,RULPTS,SCALES,NORMS)
1139 C***BEGIN PROLOGUE DINHRE
1140 C***PURPOSE DINHRE computes abscissas and weights of the integration
1141 C rule and the null rules to be used in error estimation.
1142 C These are computed as functions of NDIM and KEY.
1143 C***DESCRIPTION DINHRE will for given values of NDIM and KEY compute or
1144 C select the correct values of the abscissas and
1145 C corresponding weights for different
1146 C integration rules and null rules and assign them to
1147 C G and W.
1148 C The heuristic error coefficients ERRCOF
1149 C will be computed as a function of KEY.
1150 C Scaling factors SCALES and normalization factors NORMS
1151 C used in the error estimation are computed.
1152 C
1153 C
1154 C ON ENTRY
1155 C
1156 C NDIM Integer.
1157 C Number of variables.
1158 C KEY Integer.
1159 C Key to selected local integration rule.
1160 C WTLENG Integer.
1161 C The number of weights in each of the rules.
1162 C
1163 C ON RETURN
1164 C
1165 C W Real array of dimension (5,WTLENG).
1166 C The weights for the basic and null rules.
1167 C W(1,1), ...,W(1,WTLENG) are weights for the basic rule.
1168 C W(I,1), ...,W(I,WTLENG), for I > 1 are null rule weights.
1169 C G Real array of dimension (NDIM,WTLENG).
1170 C The fully symmetric sum generators for the rules.
1171 C G(1,J),...,G(NDIM,J) are the generators for the points
1172 C associated with the the Jth weights.
1173 C ERRCOF Real array of dimension 6.
1174 C Heuristic error coefficients that are used in the
1175 C error estimation in BASRUL.
1176 C It is assumed that the error is computed using:
1177 C IF (N1*ERRCOF(1) < N2 and N2*ERRCOF(2) < N3)
1178 C THEN ERROR = ERRCOF(3)*N1
1179 C ELSE ERROR = ERRCOF(4)*MAX(N1, N2, N3)
1180 C ERROR = ERROR + EP*(ERRCOF(5)*ERROR/(ES+ERROR)+ERRCOF(6))
1181 C where N1-N3 are the null rules, EP is the error for
1182 C the parent
1183 C subregion and ES is the error for the sibling subregion.
1184 C RULPTS Real array of dimension WTLENG.
1185 C A work array containing the number of points produced by
1186 C each generator of the selected rule.
1187 C SCALES Real array of dimension (3,WTLENG).
1188 C Scaling factors used to construct new null rules,
1189 C N1, N2 and N3,
1190 C based on a linear combination of two successive null rules
1191 C in the sequence of null rules.
1192 C NORMS Real array of dimension (3,WTLENG).
1193 C 2**NDIM/(1-norm of the null rule constructed by each of
1194 C the scaling factors.)
1195 C
1196 C***ROUTINES CALLED D132RE,D113RE,D07HRE,D09HRE
1197 C***END PROLOGUE DINHRE
1198 C
1199 C Global variables.
1200 C
1201  INTEGER NDIM,KEY,WTLENG
1202  DOUBLE PRECISION G(NDIM,WTLENG),W(5,WTLENG),ERRCOF(6)
1203  DOUBLE PRECISION RULPTS(WTLENG),SCALES(3,WTLENG)
1204  DOUBLE PRECISION NORMS(3,WTLENG)
1205 C
1206 C Local variables.
1207 C
1208  INTEGER I,J,K
1209  DOUBLE PRECISION WE(14)
1210 C
1211 C***FIRST EXECUTABLE STATEMENT DINHRE
1212 C
1213 C Compute W, G and ERRCOF.
1214 C
1215  IF (key.EQ.1) THEN
1216  CALL d132re(wtleng,w,g,errcof,rulpts)
1217  ELSE IF (key.EQ.2) THEN
1218  CALL d113re(wtleng,w,g,errcof,rulpts)
1219  ELSE IF (key.EQ.3) THEN
1220  CALL d09hre(ndim,wtleng,w,g,errcof,rulpts)
1221  ELSE IF (key.EQ.4) THEN
1222  CALL d07hre(ndim,wtleng,w,g,errcof,rulpts)
1223  END IF
1224 C
1225 C Compute SCALES and NORMS.
1226 C
1227  DO 100 k = 1,3
1228  DO 50 i = 1,wtleng
1229  IF (w(k+1,i).NE.0) THEN
1230  scales(k,i) = - w(k+2,i)/w(k+1,i)
1231  ELSE
1232  scales(k,i) = 100
1233  END IF
1234  DO 30 j = 1,wtleng
1235  we(j) = w(k+2,j) + scales(k,i)*w(k+1,j)
1236 30 CONTINUE
1237  norms(k,i) = 0
1238  DO 40 j = 1,wtleng
1239  norms(k,i) = norms(k,i) + rulpts(j)*abs(we(j))
1240 40 CONTINUE
1241  norms(k,i) = 2**ndim/norms(k,i)
1242 50 CONTINUE
1243 100 CONTINUE
1244  RETURN
1245 C
1246 C***END DINHRE
1247 C
1248  END
1249  SUBROUTINE d132re(WTLENG,W,G,ERRCOF,RULPTS)
1250 C***BEGIN PROLOGUE D132RE
1251 C***AUTHOR Jarle Berntsen, EDB-senteret,
1252 C University of Bergen, Thormohlens gt. 55,
1253 C N-5008 Bergen, NORWAY
1254 C***PURPOSE D132RE computes abscissas and weights of a 2 dimensional
1255 C integration rule of degree 13.
1256 C Two null rules of degree 11, one null rule of degree 9
1257 C and one null rule of degree 7 to be used in error
1258 C estimation are also computed.
1259 C ***DESCRIPTION D132RE will select the correct values of the abscissas
1260 C and corresponding weights for different
1261 C integration rules and null rules and assign them to
1262 C G and W. The heuristic error coefficients ERRCOF
1263 C will also be assigned.
1264 C
1265 C
1266 C ON ENTRY
1267 C
1268 C WTLENG Integer.
1269 C The number of weights in each of the rules.
1270 C
1271 C ON RETURN
1272 C
1273 C W Real array of dimension (5,WTLENG).
1274 C The weights for the basic and null rules.
1275 C W(1,1),...,W(1,WTLENG) are weights for the basic rule.
1276 C W(I,1),...,W(I,WTLENG), for I > 1 are null rule weights.
1277 C G Real array of dimension (NDIM,WTLENG).
1278 C The fully symmetric sum generators for the rules.
1279 C G(1,J),...,G(NDIM,J) are the generators for the points
1280 C associated with the the Jth weights.
1281 C ERRCOF Real array of dimension 6.
1282 C Heuristic error coefficients that are used in the
1283 C error estimation in BASRUL.
1284 C RULPTS Real array of dimension WTLENG.
1285 C The number of points produced by each generator.
1286 C***REFERENCES S.Eriksen,
1287 C Thesis of the degree cand.scient, Dept. of Informatics,
1288 C Univ. of Bergen,Norway, 1984.
1289 C
1290 C***ROUTINES CALLED-NONE
1291 C***END PROLOGUE D132RE
1292 C
1293 C Global variables
1294 C
1295  INTEGER WTLENG
1296  DOUBLE PRECISION W(5,WTLENG),G(2,WTLENG),ERRCOF(6)
1297  DOUBLE PRECISION RULPTS(WTLENG)
1298 C
1299 C Local variables.
1300 C
1301  INTEGER I,J
1302  DOUBLE PRECISION DIM2G(16)
1303  DOUBLE PRECISION DIM2W(14,5)
1304 C
1305  DATA (dim2g(i),i=1,16)/0.2517129343453109d+00,
1306  + 0.7013933644534266d+00,0.9590960631619962d+00,
1307  + 0.9956010478552127d+00,0.5000000000000000d+00,
1308  + 0.1594544658297559d+00,0.3808991135940188d+00,
1309  + 0.6582769255267192d+00,0.8761473165029315d+00,
1310  + 0.9982431840531980d+00,0.9790222658168462d+00,
1311  + 0.6492284325645389d+00,0.8727421201131239d+00,
1312  + 0.3582614645881228d+00,0.5666666666666666d+00,
1313  + 0.2077777777777778d+00/
1314 C
1315  DATA (dim2w(i,1),i=1,14)/0.3379692360134460d-01,
1316  + 0.9508589607597761d-01,0.1176006468056962d+00,
1317  + 0.2657774586326950d-01,0.1701441770200640d-01,
1318  + 0.0000000000000000d+00,0.1626593098637410d-01,
1319  + 0.1344892658526199d+00,0.1328032165460149d+00,
1320  + 0.5637474769991870d-01,0.3908279081310500d-02,
1321  + 0.3012798777432150d-01,0.1030873234689166d+00,
1322  + 0.6250000000000000d-01/
1323 C
1324  DATA (dim2w(i,2),i=1,14)/0.3213775489050763d+00,
1325  + - .1767341636743844d+00,0.7347600537466072d-01,
1326  + - .3638022004364754d-01,0.2125297922098712d-01,
1327  + 0.1460984204026913d+00,0.1747613286152099d-01,
1328  + 0.1444954045641582d+00,0.1307687976001325d-03,
1329  + 0.5380992313941161d-03,0.1042259576889814d-03,
1330  + - .1401152865045733d-02,0.8041788181514763d-02,
1331  + - .1420416552759383d+00/
1332 C
1333  DATA (dim2w(i,3),i=1,14)/0.3372900883288987d+00,
1334  + - .1644903060344491d+00,0.7707849911634622d-01,
1335  + - .3804478358506310d-01,0.2223559940380806d-01,
1336  + 0.1480693879765931d+00,0.4467143702185814d-05,
1337  + 0.1508944767074130d+00,0.3647200107516215d-04,
1338  + 0.5777198999013880d-03,0.1041757313688177d-03,
1339  + - .1452822267047819d-02,0.8338339968783705d-02,
1340  + - .1472796329231960d+00/
1341 C
1342  DATA (dim2w(i,4),i=1,14)/ - .8264123822525677d+00,
1343  + 0.3065838614094360d+00,0.2389292538329435d-02,
1344  + - .1343024157997222d+00,0.8833366840533900d-01,
1345  + 0.0000000000000000d+00,0.9786283074168292d-03,
1346  + - .1319227889147519d+00,0.7990012200150630d-02,
1347  + 0.3391747079760626d-02,0.2294915718283264d-02,
1348  + - .1358584986119197d-01,0.4025866859057809d-01,
1349  + 0.3760268580063992d-02/
1350 C
1351  DATA (dim2w(i,5),i=1,14)/0.6539094339575232d+00,
1352  + - .2041614154424632d+00, - .1746981515794990d+00,
1353  + 0.3937939671417803d-01,0.6974520545933992d-02,
1354  + 0.0000000000000000d+00,0.6667702171778258d-02,
1355  + 0.5512960621544304d-01,0.5443846381278607d-01,
1356  + 0.2310903863953934d-01,0.1506937747477189d-01,
1357  + - .6057021648901890d-01,0.4225737654686337d-01,
1358  + 0.2561989142123099d-01/
1359 C
1360 C***FIRST EXECUTABLE STATEMENT D132RE
1361 C
1362 C Assign values to W.
1363 C
1364  DO 10 i = 1,14
1365  DO 10 j = 1,5
1366  w(j,i) = dim2w(i,j)
1367 10 CONTINUE
1368 C
1369 C Assign values to G.
1370 C
1371  DO 20 i = 1,2
1372  DO 20 j = 1,14
1373  g(i,j) = 0
1374 20 CONTINUE
1375  g(1,2) = dim2g(1)
1376  g(1,3) = dim2g(2)
1377  g(1,4) = dim2g(3)
1378  g(1,5) = dim2g(4)
1379  g(1,6) = dim2g(5)
1380  g(1,7) = dim2g(6)
1381  g(2,7) = g(1,7)
1382  g(1,8) = dim2g(7)
1383  g(2,8) = g(1,8)
1384  g(1,9) = dim2g(8)
1385  g(2,9) = g(1,9)
1386  g(1,10) = dim2g(9)
1387  g(2,10) = g(1,10)
1388  g(1,11) = dim2g(10)
1389  g(2,11) = g(1,11)
1390  g(1,12) = dim2g(11)
1391  g(2,12) = dim2g(12)
1392  g(1,13) = dim2g(13)
1393  g(2,13) = dim2g(14)
1394  g(1,14) = dim2g(15)
1395  g(2,14) = dim2g(16)
1396 C
1397 C Assign values to RULPTS.
1398 C
1399  rulpts(1) = 1
1400  DO 30 i = 2,11
1401  rulpts(i) = 4
1402 30 CONTINUE
1403  rulpts(12) = 8
1404  rulpts(13) = 8
1405  rulpts(14) = 8
1406 C
1407 C Assign values to ERRCOF.
1408 C
1409  errcof(1) = 10
1410  errcof(2) = 10
1411  errcof(3) = 1.
1412  errcof(4) = 5.
1413  errcof(5) = 0.5
1414  errcof(6) = 0.25
1415 C
1416 C***END D132RE
1417 C
1418  RETURN
1419  END
1420  SUBROUTINE d113re(WTLENG,W,G,ERRCOF,RULPTS)
1421 C***BEGIN PROLOGUE D113RE
1422 C***AUTHOR Jarle Berntsen, EDB-senteret,
1423 C University of Bergen, Thormohlens gt. 55,
1424 C N-5008 Bergen, NORWAY
1425 C***PURPOSE D113RE computes abscissas and weights of a 3 dimensional
1426 C integration rule of degree 11.
1427 C Two null rules of degree 9, one null rule of degree 7
1428 C and one null rule of degree 5 to be used in error
1429 C estimation are also computed.
1430 C***DESCRIPTION D113RE will select the correct values of the abscissas
1431 C and corresponding weights for different
1432 C integration rules and null rules and assign them to G
1433 C and W.
1434 C The heuristic error coefficients ERRCOF
1435 C will also be computed.
1436 C
1437 C
1438 C ON ENTRY
1439 C
1440 C WTLENG Integer.
1441 C The number of weights in each of the rules.
1442 C
1443 C ON RETURN
1444 C
1445 C W Real array of dimension (5,WTLENG).
1446 C The weights for the basic and null rules.
1447 C W(1,1),...,W(1,WTLENG) are weights for the basic rule.
1448 C W(I,1),...,W(I,WTLENG), for I > 1 are null rule weights.
1449 C G Real array of dimension (NDIM,WTLENG).
1450 C The fully symmetric sum generators for the rules.
1451 C G(1,J),...,G(NDIM,J) are the generators for the points
1452 C associated with the the Jth weights.
1453 C ERRCOF Real array of dimension 6.
1454 C Heuristic error coefficients that are used in the
1455 C error estimation in BASRUL.
1456 C RULPTS Real array of dimension WTLENG.
1457 C The number of points used by each generator.
1458 C
1459 C***REFERENCES J.Berntsen, Cautious adaptive numerical integration
1460 C over the 3-cube, Reports in Informatics 17, Dept. of
1461 C Inf.,Univ. of Bergen, Norway, 1985.
1462 C J.Berntsen and T.O.Espelid, On the construction of
1463 C higher degree three-dimensional embedded integration
1464 C rules, SIAM J. Numer. Anal.,Vol. 25,No. 1, pp.222-234,
1465 C 1988.
1466 C***ROUTINES CALLED-NONE
1467 C***END PROLOGUE D113RE
1468 C
1469 C Global variables.
1470 C
1471  INTEGER WTLENG
1472  DOUBLE PRECISION W(5,WTLENG),G(3,WTLENG),ERRCOF(6)
1473  DOUBLE PRECISION RULPTS(WTLENG)
1474 C
1475 C Local variables.
1476 C
1477  INTEGER I,J
1478  DOUBLE PRECISION DIM3G(14)
1479  DOUBLE PRECISION DIM3W(13,5)
1480 C
1481  DATA (dim3g(i),i=1,14)/0.1900000000000000d+00,
1482  + 0.5000000000000000d+00,0.7500000000000000d+00,
1483  + 0.8000000000000000d+00,0.9949999999999999d+00,
1484  + 0.9987344998351400d+00,0.7793703685672423d+00,
1485  + 0.9999698993088767d+00,0.7902637224771788d+00,
1486  + 0.4403396687650737d+00,0.4378478459006862d+00,
1487  + 0.9549373822794593d+00,0.9661093133630748d+00,
1488  + 0.4577105877763134d+00/
1489 C
1490  DATA (dim3w(i,1),i=1,13)/0.7923078151105734d-02,
1491  + 0.6797177392788080d-01,0.1086986538805825d-02,
1492  + 0.1838633662212829d+00,0.3362119777829031d-01,
1493  + 0.1013751123334062d-01,0.1687648683985235d-02,
1494  + 0.1346468564512807d+00,0.1750145884600386d-02,
1495  + 0.7752336383837454d-01,0.2461864902770251d+00,
1496  + 0.6797944868483039d-01,0.1419962823300713d-01/
1497 C
1498  DATA (dim3w(i,2),i=1,13)/0.1715006248224684d+01,
1499  + - .3755893815889209d+00,0.1488632145140549d+00,
1500  + - .2497046640620823d+00,0.1792501419135204d+00,
1501  + 0.3446126758973890d-02, - .5140483185555825d-02,
1502  + 0.6536017839876425d-02, - .6513454939229700d-03,
1503  + - .6304672433547204d-02,0.1266959399788263d-01,
1504  + - .5454241018647931d-02,0.4826995274768427d-02/
1505 C
1506  DATA (dim3w(i,3),i=1,13)/0.1936014978949526d+01,
1507  + - .3673449403754268d+00,0.2929778657898176d-01,
1508  + - .1151883520260315d+00,0.5086658220872218d-01,
1509  + 0.4453911087786469d-01, - .2287828257125900d-01,
1510  + 0.2908926216345833d-01, - .2898884350669207d-02,
1511  + - .2805963413307495d-01,0.5638741361145884d-01,
1512  + - .2427469611942451d-01,0.2148307034182882d-01/
1513 C
1514  DATA (dim3w(i,4),i=1,13)/0.5170828195605760d+00,
1515  + 0.1445269144914044d-01, - .3601489663995932d+00,
1516  + 0.3628307003418485d+00,0.7148802650872729d-02,
1517  + - .9222852896022966d-01,0.1719339732471725d-01,
1518  + - .1021416537460350d+00, - .7504397861080493d-02,
1519  + 0.1648362537726711d-01,0.5234610158469334d-01,
1520  + 0.1445432331613066d-01,0.3019236275367777d-02/
1521 C
1522  DATA (dim3w(i,5),i=1,13)/0.2054404503818520d+01,
1523  + 0.1377759988490120d-01, - .5768062917904410d+00,
1524  + 0.3726835047700328d-01,0.6814878939777219d-02,
1525  + 0.5723169733851849d-01, - .4493018743811285d-01,
1526  + 0.2729236573866348d-01,0.3547473950556990d-03,
1527  + 0.1571366799739551d-01,0.4990099219278567d-01,
1528  + 0.1377915552666770d-01,0.2878206423099872d-02/
1529 C
1530 C***FIRST EXECUTABLE STATEMENT D113RE
1531 C
1532 C Assign values to W.
1533 C
1534  DO 10 i = 1,13
1535  DO 10 j = 1,5
1536  w(j,i) = dim3w(i,j)
1537 10 CONTINUE
1538 C
1539 C Assign values to G.
1540 C
1541  DO 20 i = 1,3
1542  DO 20 j = 1,13
1543  g(i,j) = 0
1544 20 CONTINUE
1545  g(1,2) = dim3g(1)
1546  g(1,3) = dim3g(2)
1547  g(1,4) = dim3g(3)
1548  g(1,5) = dim3g(4)
1549  g(1,6) = dim3g(5)
1550  g(1,7) = dim3g(6)
1551  g(2,7) = g(1,7)
1552  g(1,8) = dim3g(7)
1553  g(2,8) = g(1,8)
1554  g(1,9) = dim3g(8)
1555  g(2,9) = g(1,9)
1556  g(3,9) = g(1,9)
1557  g(1,10) = dim3g(9)
1558  g(2,10) = g(1,10)
1559  g(3,10) = g(1,10)
1560  g(1,11) = dim3g(10)
1561  g(2,11) = g(1,11)
1562  g(3,11) = g(1,11)
1563  g(1,12) = dim3g(12)
1564  g(2,12) = dim3g(11)
1565  g(3,12) = g(2,12)
1566  g(1,13) = dim3g(13)
1567  g(2,13) = g(1,13)
1568  g(3,13) = dim3g(14)
1569 C
1570 C Assign values to RULPTS.
1571 C
1572  rulpts(1) = 1
1573  rulpts(2) = 6
1574  rulpts(3) = 6
1575  rulpts(4) = 6
1576  rulpts(5) = 6
1577  rulpts(6) = 6
1578  rulpts(7) = 12
1579  rulpts(8) = 12
1580  rulpts(9) = 8
1581  rulpts(10) = 8
1582  rulpts(11) = 8
1583  rulpts(12) = 24
1584  rulpts(13) = 24
1585 C
1586 C Assign values to ERRCOF.
1587 C
1588  errcof(1) = 4
1589  errcof(2) = 4.
1590  errcof(3) = 0.5
1591  errcof(4) = 3.
1592  errcof(5) = 0.5
1593  errcof(6) = 0.25
1594 C
1595 C***END D113RE
1596 C
1597  RETURN
1598  END
1599  SUBROUTINE d09hre(NDIM,WTLENG,W,G,ERRCOF,RULPTS)
1600 C***BEGIN PROLOGUE D09HRE
1601 C***KEYWORDS basic integration rule, degree 9
1602 C***PURPOSE To initialize a degree 9 basic rule and null rules.
1603 C***AUTHOR Alan Genz, Computer Science Department, Washington
1604 C State University, Pullman, WA 99163-1210 USA
1605 C***LAST MODIFICATION 88-05-20
1606 C***DESCRIPTION D09HRE initializes a degree 9 integration rule,
1607 C two degree 7 null rules, one degree 5 null rule and one
1608 C degree 3 null rule for the hypercube [-1,1]**NDIM.
1609 C
1610 C ON ENTRY
1611 C
1612 C NDIM Integer.
1613 C Number of variables.
1614 C WTLENG Integer.
1615 C The number of weights in each of the rules.
1616 C
1617 C ON RETURN
1618 C W Real array of dimension (5,WTLENG).
1619 C The weights for the basic and null rules.
1620 C W(1,1),...,W(1,WTLENG) are weights for the basic rule.
1621 C W(I,1),...,W(I,WTLENG), for I > 1 are null rule weights.
1622 C G Real array of dimension (NDIM, WTLENG).
1623 C The fully symmetric sum generators for the rules.
1624 C G(1, J), ..., G(NDIM, J) are the are the generators for the
1625 C points associated with the Jth weights.
1626 C ERRCOF Real array of dimension 6.
1627 C Heuristic error coefficients that are used in the
1628 C error estimation in BASRUL.
1629 C RULPTS Real array of dimension WTLENG.
1630 C A work array.
1631 C
1632 C***REFERENCES A. Genz and A. Malik,
1633 C "An Imbedded Family of Fully Symmetric Numerical
1634 C Integration Rules",
1635 C SIAM J Numer. Anal. 20 (1983), pp. 580-588.
1636 C***ROUTINES CALLED-NONE
1637 C***END PROLOGUE D09HRE
1638 C
1639 C Global variables
1640 C
1641  INTEGER NDIM,WTLENG
1642  DOUBLE PRECISION W(5,WTLENG),G(NDIM,WTLENG),ERRCOF(6)
1643  DOUBLE PRECISION RULPTS(WTLENG)
1644 C
1645 C Local Variables
1646 C
1647  DOUBLE PRECISION RATIO,LAM0,LAM1,LAM2,LAM3,LAMP,TWONDM
1648  INTEGER I,J
1649 C
1650 C***FIRST EXECUTABLE STATEMENT D09HRE
1651 C
1652 C
1653 C Initialize generators, weights and RULPTS
1654 C
1655  DO 30 j = 1,wtleng
1656  DO 10 i = 1,ndim
1657  g(i,j) = 0
1658 10 CONTINUE
1659  DO 20 i = 1,5
1660  w(i,j) = 0
1661 20 CONTINUE
1662  rulpts(j) = 2*ndim
1663 30 CONTINUE
1664  twondm = 2**ndim
1665  rulpts(wtleng) = twondm
1666  IF (ndim.GT.2) rulpts(8) = (4*ndim* (ndim-1)* (ndim-2))/3
1667  rulpts(7) = 4*ndim* (ndim-1)
1668  rulpts(6) = 2*ndim* (ndim-1)
1669  rulpts(1) = 1
1670 C
1671 C Compute squared generator parameters
1672 C
1673  lam0 = 0.4707
1674  lam1 = 4/ (15-5/lam0)
1675  ratio = (1-lam1/lam0)/27
1676  lam2 = (5-7*lam1-35*ratio)/ (7-35*lam1/3-35*ratio/lam0)
1677  ratio = ratio* (1-lam2/lam0)/3
1678  lam3 = (7-9* (lam2+lam1)+63*lam2*lam1/5-63*ratio)/
1679  + (9-63* (lam2+lam1)/5+21*lam2*lam1-63*ratio/lam0)
1680  lamp = 0.0625
1681 C
1682 C Compute degree 9 rule weights
1683 C
1684  w(1,wtleng) = 1/ (3*lam0)**4/twondm
1685  IF (ndim.GT.2) w(1,8) = (1-1/ (3*lam0))/ (6*lam1)**3
1686  w(1,7) = (1-7* (lam0+lam1)/5+7*lam0*lam1/3)/
1687  + (84*lam1*lam2* (lam2-lam0)* (lam2-lam1))
1688  w(1,6) = (1-7* (lam0+lam2)/5+7*lam0*lam2/3)/
1689  + (84*lam1*lam1* (lam1-lam0)* (lam1-lam2)) -
1690  + w(1,7)*lam2/lam1 - 2* (ndim-2)*w(1,8)
1691  w(1,4) = (1-9* ((lam0+lam1+lam2)/7- (lam0*lam1+lam0*lam2+
1692  + lam1*lam2)/5)-3*lam0*lam1*lam2)/
1693  + (18*lam3* (lam3-lam0)* (lam3-lam1)* (lam3-lam2))
1694  w(1,3) = (1-9* ((lam0+lam1+lam3)/7- (lam0*lam1+lam0*lam3+
1695  + lam1*lam3)/5)-3*lam0*lam1*lam3)/
1696  + (18*lam2* (lam2-lam0)* (lam2-lam1)* (lam2-lam3)) -
1697  + 2* (ndim-1)*w(1,7)
1698  w(1,2) = (1-9* ((lam0+lam2+lam3)/7- (lam0*lam2+lam0*lam3+
1699  + lam2*lam3)/5)-3*lam0*lam2*lam3)/
1700  + (18*lam1* (lam1-lam0)* (lam1-lam2)* (lam1-lam3)) -
1701  + 2* (ndim-1)* (w(1,7)+w(1,6)+ (ndim-2)*w(1,8))
1702 C
1703 C Compute weights for 2 degree 7, 1 degree 5 and 1 degree 3 rules
1704 C
1705  w(2,wtleng) = 1/ (108*lam0**4)/twondm
1706  IF (ndim.GT.2) w(2,8) = (1-27*twondm*w(2,9)*lam0**3)/ (6*lam1)**3
1707  w(2,7) = (1-5*lam1/3-15*twondm*w(2,wtleng)*lam0**2* (lam0-lam1))/
1708  + (60*lam1*lam2* (lam2-lam1))
1709  w(2,6) = (1-9* (8*lam1*lam2*w(2,7)+twondm*w(2,wtleng)*lam0**2))/
1710  + (36*lam1*lam1) - 2*w(2,8)* (ndim-2)
1711  w(2,4) = (1-7* ((lam1+lam2)/5-lam1*lam2/3+twondm*w(2,
1712  + wtleng)*lam0* (lam0-lam1)* (lam0-lam2)))/
1713  + (14*lam3* (lam3-lam1)* (lam3-lam2))
1714  w(2,3) = (1-7* ((lam1+lam3)/5-lam1*lam3/3+twondm*w(2,
1715  + wtleng)*lam0* (lam0-lam1)* (lam0-lam3)))/
1716  + (14*lam2* (lam2-lam1)* (lam2-lam3)) - 2* (ndim-1)*w(2,7)
1717  w(2,2) = (1-7* ((lam2+lam3)/5-lam2*lam3/3+twondm*w(2,
1718  + wtleng)*lam0* (lam0-lam2)* (lam0-lam3)))/
1719  + (14*lam1* (lam1-lam2)* (lam1-lam3)) -
1720  + 2* (ndim-1)* (w(2,7)+w(2,6)+ (ndim-2)*w(2,8))
1721  w(3,wtleng) = 5/ (324*lam0**4)/twondm
1722  IF (ndim.GT.2) w(3,8) = (1-27*twondm*w(3,9)*lam0**3)/ (6*lam1)**3
1723  w(3,7) = (1-5*lam1/3-15*twondm*w(3,wtleng)*lam0**2* (lam0-lam1))/
1724  + (60*lam1*lam2* (lam2-lam1))
1725  w(3,6) = (1-9* (8*lam1*lam2*w(3,7)+twondm*w(3,wtleng)*lam0**2))/
1726  + (36*lam1*lam1) - 2*w(3,8)* (ndim-2)
1727  w(3,5) = (1-7* ((lam1+lam2)/5-lam1*lam2/3+twondm*w(3,
1728  + wtleng)*lam0* (lam0-lam1)* (lam0-lam2)))/
1729  + (14*lamp* (lamp-lam1)* (lamp-lam2))
1730  w(3,3) = (1-7* ((lam1+lamp)/5-lam1*lamp/3+twondm*w(3,
1731  + wtleng)*lam0* (lam0-lam1)* (lam0-lamp)))/
1732  + (14*lam2* (lam2-lam1)* (lam2-lamp)) - 2* (ndim-1)*w(3,7)
1733  w(3,2) = (1-7* ((lam2+lamp)/5-lam2*lamp/3+twondm*w(3,
1734  + wtleng)*lam0* (lam0-lam2)* (lam0-lamp)))/
1735  + (14*lam1* (lam1-lam2)* (lam1-lamp)) -
1736  + 2* (ndim-1)* (w(3,7)+w(3,6)+ (ndim-2)*w(3,8))
1737  w(4,wtleng) = 2/ (81*lam0**4)/twondm
1738  IF (ndim.GT.2) w(4,8) = (2-27*twondm*w(4,9)*lam0**3)/ (6*lam1)**3
1739  w(4,7) = (2-15*lam1/9-15*twondm*w(4,wtleng)*lam0* (lam0-lam1))/
1740  + (60*lam1*lam2* (lam2-lam1))
1741  w(4,6) = (1-9* (8*lam1*lam2*w(4,7)+twondm*w(4,wtleng)*lam0**2))/
1742  + (36*lam1*lam1) - 2*w(4,8)* (ndim-2)
1743  w(4,4) = (2-7* ((lam1+lam2)/5-lam1*lam2/3+twondm*w(4,
1744  + wtleng)*lam0* (lam0-lam1)* (lam0-lam2)))/
1745  + (14*lam3* (lam3-lam1)* (lam3-lam2))
1746  w(4,3) = (2-7* ((lam1+lam3)/5-lam1*lam3/3+twondm*w(4,
1747  + wtleng)*lam0* (lam0-lam1)* (lam0-lam3)))/
1748  + (14*lam2* (lam2-lam1)* (lam2-lam3)) - 2* (ndim-1)*w(4,7)
1749  w(4,2) = (2-7* ((lam2+lam3)/5-lam2*lam3/3+twondm*w(4,
1750  + wtleng)*lam0* (lam0-lam2)* (lam0-lam3)))/
1751  + (14*lam1* (lam1-lam2)* (lam1-lam3)) -
1752  + 2* (ndim-1)* (w(4,7)+w(4,6)+ (ndim-2)*w(4,8))
1753  w(5,2) = 1/ (6*lam1)
1754 C
1755 C Set generator values
1756 C
1757  lam0 = sqrt(lam0)
1758  lam1 = sqrt(lam1)
1759  lam2 = sqrt(lam2)
1760  lam3 = sqrt(lam3)
1761  lamp = sqrt(lamp)
1762  DO 40 i = 1,ndim
1763  g(i,wtleng) = lam0
1764 40 CONTINUE
1765  IF (ndim.GT.2) THEN
1766  g(1,8) = lam1
1767  g(2,8) = lam1
1768  g(3,8) = lam1
1769  END IF
1770  g(1,7) = lam1
1771  g(2,7) = lam2
1772  g(1,6) = lam1
1773  g(2,6) = lam1
1774  g(1,5) = lamp
1775  g(1,4) = lam3
1776  g(1,3) = lam2
1777  g(1,2) = lam1
1778 C
1779 C Compute final weight values.
1780 C The null rule weights are computed from differences between
1781 C the degree 9 rule weights and lower degree rule weights.
1782 C
1783  w(1,1) = twondm
1784  DO 70 j = 2,5
1785  DO 50 i = 2,wtleng
1786  w(j,i) = w(j,i) - w(1,i)
1787  w(j,1) = w(j,1) - rulpts(i)*w(j,i)
1788 50 CONTINUE
1789 70 CONTINUE
1790  DO 80 i = 2,wtleng
1791  w(1,i) = twondm*w(1,i)
1792  w(1,1) = w(1,1) - rulpts(i)*w(1,i)
1793 80 CONTINUE
1794 C
1795 C Set error coefficients
1796 C
1797  errcof(1) = 5
1798  errcof(2) = 5
1799  errcof(3) = 1.
1800  errcof(4) = 5
1801  errcof(5) = 0.5
1802  errcof(6) = 0.25
1803 C
1804 C***END D09HRE
1805 C
1806  END
1807  SUBROUTINE d07hre(NDIM,WTLENG,W,G,ERRCOF,RULPTS)
1808 C***BEGIN PROLOGUE D07HRE
1809 C***KEYWORDS basic integration rule, degree 7
1810 C***PURPOSE To initialize a degree 7 basic rule, and null rules.
1811 C***AUTHOR Alan Genz, Computer Science Department, Washington
1812 C State University, Pullman, WA 99163-1210 USA
1813 C***LAST MODIFICATION 88-05-31
1814 C***DESCRIPTION D07HRE initializes a degree 7 integration rule,
1815 C two degree 5 null rules, one degree 3 null rule and one
1816 C degree 1 null rule for the hypercube [-1,1]**NDIM.
1817 C
1818 C ON ENTRY
1819 C
1820 C NDIM Integer.
1821 C Number of variables.
1822 C WTLENG Integer.
1823 C The number of weights in each of the rules.
1824 C WTLENG MUST be set equal to 6.
1825 C
1826 C ON RETURN
1827 C W Real array of dimension (5,WTLENG).
1828 C The weights for the basic and null rules.
1829 C W(1,1),...,W(1,WTLENG) are weights for the basic rule.
1830 C W(I,1),...,W(I,WTLENG), for I > 1 are null rule weights.
1831 C G Real array of dimension (NDIM, WTLENG).
1832 C The fully symmetric sum generators for the rules.
1833 C G(1, J), ..., G(NDIM, J) are the are the generators for the
1834 C points associated with the Jth weights.
1835 C ERRCOF Real array of dimension 6.
1836 C Heuristic error coefficients that are used in the
1837 C error estimation in BASRUL.
1838 C RULPTS Real array of dimension WTLENG.
1839 C A work array.
1840 C
1841 C***REFERENCES A. Genz and A. Malik,
1842 C "An Imbedded Family of Fully Symmetric Numerical
1843 C Integration Rules",
1844 C SIAM J Numer. Anal. 20 (1983), pp. 580-588.
1845 C***ROUTINES CALLED-NONE
1846 C***END PROLOGUE D07HRE
1847 C
1848 C Global variables
1849 C
1850  INTEGER NDIM,WTLENG
1851  DOUBLE PRECISION W(5,WTLENG),G(NDIM,WTLENG),ERRCOF(6)
1852  DOUBLE PRECISION RULPTS(WTLENG)
1853 C
1854 C Local Variables
1855 C
1856  DOUBLE PRECISION RATIO,LAM0,LAM1,LAM2,LAMP,TWONDM
1857  INTEGER I,J
1858 C
1859 C***FIRST EXECUTABLE STATEMENT D07HRE
1860 C
1861 C
1862 C Initialize generators, weights and RULPTS
1863 C
1864  DO 30 j = 1,wtleng
1865  DO 10 i = 1,ndim
1866  g(i,j) = 0
1867 10 CONTINUE
1868  DO 20 i = 1,5
1869  w(i,j) = 0
1870 20 CONTINUE
1871  rulpts(j) = 2*ndim
1872 30 CONTINUE
1873  twondm = 2**ndim
1874  rulpts(wtleng) = twondm
1875  rulpts(wtleng-1) = 2*ndim* (ndim-1)
1876  rulpts(1) = 1
1877 C
1878 C Compute squared generator parameters
1879 C
1880  lam0 = 0.4707
1881  lamp = 0.5625
1882  lam1 = 4/ (15-5/lam0)
1883  ratio = (1-lam1/lam0)/27
1884  lam2 = (5-7*lam1-35*ratio)/ (7-35*lam1/3-35*ratio/lam0)
1885 C
1886 C Compute degree 7 rule weights
1887 C
1888  w(1,6) = 1/ (3*lam0)**3/twondm
1889  w(1,5) = (1-5*lam0/3)/ (60* (lam1-lam0)*lam1**2)
1890  w(1,3) = (1-5*lam2/3-5*twondm*w(1,6)*lam0* (lam0-lam2))/
1891  + (10*lam1* (lam1-lam2)) - 2* (ndim-1)*w(1,5)
1892  w(1,2) = (1-5*lam1/3-5*twondm*w(1,6)*lam0* (lam0-lam1))/
1893  + (10*lam2* (lam2-lam1))
1894 C
1895 C Compute weights for 2 degree 5, 1 degree 3 and 1 degree 1 rules
1896 C
1897  w(2,6) = 1/ (36*lam0**3)/twondm
1898  w(2,5) = (1-9*twondm*w(2,6)*lam0**2)/ (36*lam1**2)
1899  w(2,3) = (1-5*lam2/3-5*twondm*w(2,6)*lam0* (lam0-lam2))/
1900  + (10*lam1* (lam1-lam2)) - 2* (ndim-1)*w(2,5)
1901  w(2,2) = (1-5*lam1/3-5*twondm*w(2,6)*lam0* (lam0-lam1))/
1902  + (10*lam2* (lam2-lam1))
1903  w(3,6) = 5/ (108*lam0**3)/twondm
1904  w(3,5) = (1-9*twondm*w(3,6)*lam0**2)/ (36*lam1**2)
1905  w(3,3) = (1-5*lamp/3-5*twondm*w(3,6)*lam0* (lam0-lamp))/
1906  + (10*lam1* (lam1-lamp)) - 2* (ndim-1)*w(3,5)
1907  w(3,4) = (1-5*lam1/3-5*twondm*w(3,6)*lam0* (lam0-lam1))/
1908  + (10*lamp* (lamp-lam1))
1909  w(4,6) = 1/ (54*lam0**3)/twondm
1910  w(4,5) = (1-18*twondm*w(4,6)*lam0**2)/ (72*lam1**2)
1911  w(4,3) = (1-10*lam2/3-10*twondm*w(4,6)*lam0* (lam0-lam2))/
1912  + (20*lam1* (lam1-lam2)) - 2* (ndim-1)*w(4,5)
1913  w(4,2) = (1-10*lam1/3-10*twondm*w(4,6)*lam0* (lam0-lam1))/
1914  + (20*lam2* (lam2-lam1))
1915 C
1916 C Set generator values
1917 C
1918  lam0 = sqrt(lam0)
1919  lam1 = sqrt(lam1)
1920  lam2 = sqrt(lam2)
1921  lamp = sqrt(lamp)
1922  DO 40 i = 1,ndim
1923  g(i,wtleng) = lam0
1924 40 CONTINUE
1925  g(1,wtleng-1) = lam1
1926  g(2,wtleng-1) = lam1
1927  g(1,wtleng-4) = lam2
1928  g(1,wtleng-3) = lam1
1929  g(1,wtleng-2) = lamp
1930 C
1931 C Compute final weight values.
1932 C The null rule weights are computed from differences between
1933 C the degree 7 rule weights and lower degree rule weights.
1934 C
1935  w(1,1) = twondm
1936  DO 70 j = 2,5
1937  DO 50 i = 2,wtleng
1938  w(j,i) = w(j,i) - w(1,i)
1939  w(j,1) = w(j,1) - rulpts(i)*w(j,i)
1940 50 CONTINUE
1941 70 CONTINUE
1942  DO 80 i = 2,wtleng
1943  w(1,i) = twondm*w(1,i)
1944  w(1,1) = w(1,1) - rulpts(i)*w(1,i)
1945 80 CONTINUE
1946 C
1947 C Set error coefficients
1948 C
1949  errcof(1) = 5
1950  errcof(2) = 5
1951  errcof(3) = 1
1952  errcof(4) = 5
1953  errcof(5) = 0.5
1954  errcof(6) = 0.25
1955 C
1956 C***END D07HRE
1957 C
1958  END
1959  SUBROUTINE drlhre(NDIM,CENTER,HWIDTH,WTLENG,G,W,ERRCOF,NUMFUN,
1960  + FUNSUB,SCALES,NORMS,X,NULL,BASVAL,RGNERR,DIRECT)
1961 C***BEGIN PROLOGUE DRLHRE
1962 C***KEYWORDS basic numerical integration rule
1963 C***PURPOSE To compute basic integration rule values.
1964 C***AUTHOR Alan Genz, Computer Science Department, Washington
1965 C State University, Pullman, WA 99163-1210 USA
1966 C***LAST MODIFICATION 90-02-06
1967 C***DESCRIPTION DRLHRE computes basic integration rule values for a
1968 C vector of integrands over a hyper-rectangular region.
1969 C These are estimates for the integrals. DRLHRE also computes
1970 C estimates for the errors and determines the coordinate axis
1971 C where the fourth difference for the integrands is largest.
1972 C
1973 C ON ENTRY
1974 C
1975 C NDIM Integer.
1976 C Number of variables.
1977 C CENTER Real array of dimension NDIM.
1978 C The coordinates for the center of the region.
1979 C HWIDTH Real Array of dimension NDIM.
1980 C HWIDTH(I) is half of the width of dimension I of the region.
1981 C WTLENG Integer.
1982 C The number of weights in the basic integration rule.
1983 C G Real array of dimension (NDIM,WTLENG).
1984 C The fully symmetric sum generators for the rules.
1985 C G(1,J), ..., G(NDIM,J) are the are the generators for the
1986 C points associated with the Jth weights.
1987 C W Real array of dimension (5,WTLENG).
1988 C The weights for the basic and null rules.
1989 C W(1,1),...,W(1,WTLENG) are weights for the basic rule.
1990 C W(I,1),...,W(I,WTLENG), for I > 1 are null rule weights.
1991 C ERRCOF Real array of dimension 6.
1992 C The error coefficients for the rules.
1993 C It is assumed that the error is computed using:
1994 C IF (N1*ERRCOF(1) < N2 and N2*ERRCOF(2) < N3)
1995 C THEN ERROR = ERRCOF(3)*N1
1996 C ELSE ERROR = ERRCOF(4)*MAX(N1, N2, N3)
1997 C ERROR = ERROR + EP*(ERRCOF(5)*ERROR/(ES+ERROR)+ERRCOF(6))
1998 C where N1-N4 are the null rules, EP is the error
1999 C for the parent
2000 C subregion and ES is the error for the sibling subregion.
2001 C NUMFUN Integer.
2002 C Number of components for the vector integrand.
2003 C FUNSUB Externally declared subroutine.
2004 C For computing the components of the integrand at a point X.
2005 C It must have parameters (NDIM,X,NUMFUN,FUNVLS).
2006 C Input Parameters:
2007 C X Real array of dimension NDIM.
2008 C Defines the evaluation point.
2009 C NDIM Integer.
2010 C Number of variables for the integrand.
2011 C NUMFUN Integer.
2012 C Number of components for the vector integrand.
2013 C Output Parameters:
2014 C FUNVLS Real array of dimension NUMFUN.
2015 C The components of the integrand at the point X.
2016 C SCALES Real array of dimension (3,WTLENG).
2017 C Scaling factors used to construct new null rules based
2018 C on a linear combination of two successive null rules
2019 C in the sequence of null rules.
2020 C NORMS Real array of dimension (3,WTLENG).
2021 C 2**NDIM/(1-norm of the null rule constructed by each of the
2022 C scaling factors.)
2023 C X Real Array of dimension NDIM.
2024 C A work array.
2025 C NULL Real array of dimension (NUMFUN, 8)
2026 C A work array.
2027 C
2028 C ON RETURN
2029 C
2030 C BASVAL Real array of dimension NUMFUN.
2031 C The values for the basic rule for each component
2032 C of the integrand.
2033 C RGNERR Real array of dimension NUMFUN.
2034 C The error estimates for each component of the integrand.
2035 C DIRECT Real.
2036 C The coordinate axis where the fourth difference of the
2037 C integrand values is largest.
2038 C
2039 C***REFERENCES
2040 C A.C.Genz and A.A.Malik, An adaptive algorithm for numerical
2041 C integration over an N-dimensional rectangular region,
2042 C J.Comp.Appl.Math., 6:295-302, 1980.
2043 C
2044 C T.O.Espelid, Integration Rules, Null Rules and Error
2045 C Estimation, Reports in Informatics 33, Dept. of Informatics,
2046 C Univ. of Bergen, 1988.
2047 C
2048 C***ROUTINES CALLED: DFSHRE, FUNSUB
2049 C
2050 C***END PROLOGUE DRLHRE
2051 C
2052 C Global variables.
2053 C
2054  EXTERNAL funsub
2055  INTEGER WTLENG,NUMFUN,NDIM
2056  DOUBLE PRECISION CENTER(NDIM),X(NDIM),HWIDTH(NDIM),BASVAL(NUMFUN),
2057  + RGNERR(NUMFUN),NULL(NUMFUN,8),W(5,WTLENG),
2058  + G(NDIM,WTLENG),ERRCOF(6),DIRECT,SCALES(3,WTLENG),
2059  + NORMS(3,WTLENG)
2060 C
2061 C Local variables.
2062 C
2063  DOUBLE PRECISION RGNVOL,DIFSUM,DIFMAX,FRTHDF
2064  INTEGER I,J,K,DIVAXN
2065  DOUBLE PRECISION SEARCH,RATIO
2066 C
2067 C***FIRST EXECUTABLE STATEMENT DRLHRE
2068 C
2069 C
2070 C Compute volume of subregion, initialize DIVAXN and rule sums;
2071 C compute fourth differences and new DIVAXN (RGNERR is used
2072 C for a work array here). The integrand values used for the
2073 C fourth divided differences are accumulated in rule arrays.
2074 C
2075  rgnvol = 1
2076  divaxn = 1
2077  DO 10 i = 1,ndim
2078  rgnvol = rgnvol*hwidth(i)
2079  x(i) = center(i)
2080  IF (hwidth(i).GT.hwidth(divaxn)) divaxn = i
2081 10 CONTINUE
2082  CALL funsub(ndim,x,numfun,rgnerr)
2083  DO 30 j = 1,numfun
2084  basval(j) = w(1,1)*rgnerr(j)
2085  DO 20 k = 1,4
2086  null(j,k) = w(k+1,1)*rgnerr(j)
2087 20 CONTINUE
2088 30 CONTINUE
2089  difmax = 0
2090  ratio = (g(1,3)/g(1,2))**2
2091  DO 60 i = 1,ndim
2092  x(i) = center(i) - hwidth(i)*g(1,2)
2093  CALL funsub(ndim,x,numfun,null(1,5))
2094  x(i) = center(i) + hwidth(i)*g(1,2)
2095  CALL funsub(ndim,x,numfun,null(1,6))
2096  x(i) = center(i) - hwidth(i)*g(1,3)
2097  CALL funsub(ndim,x,numfun,null(1,7))
2098  x(i) = center(i) + hwidth(i)*g(1,3)
2099  CALL funsub(ndim,x,numfun,null(1,8))
2100  x(i) = center(i)
2101  difsum = 0
2102  DO 50 j = 1,numfun
2103  frthdf = 2* (1-ratio)*rgnerr(j) - (null(j,7)+null(j,8)) +
2104  + ratio* (null(j,5)+null(j,6))
2105 C
2106 C Ignore differences below roundoff
2107 C
2108  IF (rgnerr(j)+frthdf/4.NE.rgnerr(j)) difsum = difsum +
2109  + abs(frthdf)
2110  DO 40 k = 1,4
2111  null(j,k) = null(j,k) + w(k+1,2)*
2112  + (null(j,5)+null(j,6)) +
2113  + w(k+1,3)* (null(j,7)+null(j,8))
2114 40 CONTINUE
2115  basval(j) = basval(j) + w(1,2)* (null(j,5)+null(j,6)) +
2116  + w(1,3)* (null(j,7)+null(j,8))
2117 50 CONTINUE
2118  IF (difsum.GT.difmax) THEN
2119  difmax = difsum
2120  divaxn = i
2121  END IF
2122 60 CONTINUE
2123  direct = divaxn
2124 C
2125 C Finish computing the rule values.
2126 C
2127  DO 90 i = 4,wtleng
2128  CALL dfshre(ndim,center,hwidth,x,g(1,i),numfun,funsub,rgnerr,
2129  + null(1,5))
2130  DO 80 j = 1,numfun
2131  basval(j) = basval(j) + w(1,i)*rgnerr(j)
2132  DO 70 k = 1,4
2133  null(j,k) = null(j,k) + w(k+1,i)*rgnerr(j)
2134 70 CONTINUE
2135 80 CONTINUE
2136 90 CONTINUE
2137 C
2138 C Compute errors.
2139 C
2140  DO 130 j = 1,numfun
2141 C
2142 C We search for the null rule, in the linear space spanned by two
2143 C successive null rules in our sequence, which gives the greatest
2144 C error estimate among all normalized (1-norm) null rules in this
2145 C space.
2146 C
2147  DO 110 i = 1,3
2148  search = 0
2149  DO 100 k = 1,wtleng
2150  search = max(search,abs(null(j,i+1)+scales(i,
2151  + k)*null(j,i))*norms(i,k))
2152 100 CONTINUE
2153  null(j,i) = search
2154 110 CONTINUE
2155  IF (errcof(1)*null(j,1).LE.null(j,2) .AND.
2156  + errcof(2)*null(j,2).LE.null(j,3)) THEN
2157  rgnerr(j) = errcof(3)*null(j,1)
2158  ELSE
2159  rgnerr(j) = errcof(4)*max(null(j,1),null(j,2),null(j,3))
2160  END IF
2161  rgnerr(j) = rgnvol*rgnerr(j)
2162  basval(j) = rgnvol*basval(j)
2163 130 CONTINUE
2164 C
2165 C***END DRLHRE
2166 C
2167  END
2168  SUBROUTINE dfshre(NDIM,CENTER,HWIDTH,X,G,NUMFUN,FUNSUB,FULSMS,
2169  + FUNVLS)
2170 C***BEGIN PROLOGUE DFSHRE
2171 C***KEYWORDS fully symmetric sum
2172 C***PURPOSE To compute fully symmetric basic rule sums
2173 C***AUTHOR Alan Genz, Computer Science Department, Washington
2174 C State University, Pullman, WA 99163-1210 USA
2175 C***LAST MODIFICATION 88-04-08
2176 C***DESCRIPTION DFSHRE computes a fully symmetric sum for a vector
2177 C of integrand values over a hyper-rectangular region.
2178 C The sum is fully symmetric with respect to the center of
2179 C the region and is taken over all sign changes and
2180 C permutations of the generators for the sum.
2181 C
2182 C ON ENTRY
2183 C
2184 C NDIM Integer.
2185 C Number of variables.
2186 C CENTER Real array of dimension NDIM.
2187 C The coordinates for the center of the region.
2188 C HWIDTH Real Array of dimension NDIM.
2189 C HWIDTH(I) is half of the width of dimension I of the region.
2190 C X Real Array of dimension NDIM.
2191 C A work array.
2192 C G Real Array of dimension NDIM.
2193 C The generators for the fully symmetric sum. These MUST BE
2194 C non-negative and non-increasing.
2195 C NUMFUN Integer.
2196 C Number of components for the vector integrand.
2197 C FUNSUB Externally declared subroutine.
2198 C For computing the components of the integrand at a point X.
2199 C It must have parameters (NDIM, X, NUMFUN, FUNVLS).
2200 C Input Parameters:
2201 C X Real array of dimension NDIM.
2202 C Defines the evaluation point.
2203 C NDIM Integer.
2204 C Number of variables for the integrand.
2205 C NUMFUN Integer.
2206 C Number of components for the vector integrand.
2207 C Output Parameters:
2208 C FUNVLS Real array of dimension NUMFUN.
2209 C The components of the integrand at the point X.
2210 C ON RETURN
2211 C
2212 C FULSMS Real array of dimension NUMFUN.
2213 C The values for the fully symmetric sums for each component
2214 C of the integrand.
2215 C FUNVLS Real array of dimension NUMFUN.
2216 C A work array.
2217 C
2218 C***ROUTINES CALLED: FUNSUB
2219 C
2220 C***END PROLOGUE DFSHRE
2221 C
2222 C Global variables.
2223 C
2224  EXTERNAL funsub
2225  INTEGER NDIM,NUMFUN
2226  DOUBLE PRECISION CENTER(NDIM),HWIDTH(NDIM),X(NDIM),G(NDIM),
2227  + FULSMS(NUMFUN),FUNVLS(NUMFUN)
2228 C
2229 C Local variables.
2230 C
2231  INTEGER IXCHNG,LXCHNG,I,J,L
2232  DOUBLE PRECISION GL,GI
2233 C
2234 C***FIRST EXECUTABLE STATEMENT DFSHRE
2235 C
2236  DO 10 j = 1,numfun
2237  fulsms(j) = 0
2238 10 CONTINUE
2239 C
2240 C Compute centrally symmetric sum for permutation of G
2241 C
2242 20 DO 30 i = 1,ndim
2243  x(i) = center(i) + g(i)*hwidth(i)
2244 30 CONTINUE
2245 40 CALL funsub(ndim,x,numfun,funvls)
2246  DO 50 j = 1,numfun
2247  fulsms(j) = fulsms(j) + funvls(j)
2248 50 CONTINUE
2249  DO 60 i = 1,ndim
2250  g(i) = - g(i)
2251  x(i) = center(i) + g(i)*hwidth(i)
2252  IF (g(i).LT.0) GO TO 40
2253 60 CONTINUE
2254 C
2255 C Find next distinct permutation of G and loop back for next sum.
2256 C Permutations are generated in reverse lexicographic order.
2257 C
2258  DO 80 i = 2,ndim
2259  IF (g(i-1).GT.g(i)) THEN
2260  gi = g(i)
2261  ixchng = i - 1
2262  DO 70 l = 1, (i-1)/2
2263  gl = g(l)
2264  g(l) = g(i-l)
2265  g(i-l) = gl
2266  IF (gl.LE.gi) ixchng = ixchng - 1
2267  IF (g(l).GT.gi) lxchng = l
2268 70 CONTINUE
2269  IF (g(ixchng).LE.gi) ixchng = lxchng
2270  g(i) = g(ixchng)
2271  g(ixchng) = gi
2272  GO TO 20
2273  END IF
2274 80 CONTINUE
2275 C
2276 C Restore original order to generators
2277 C
2278  DO 90 i = 1,ndim/2
2279  gi = g(i)
2280  g(i) = g(ndim-i+1)
2281  g(ndim-i+1) = gi
2282 90 CONTINUE
2283 C
2284 C***END DFSHRE
2285 C
2286  END
2287  SUBROUTINE dtrhre(DVFLAG,NDIM,NUMFUN,SBRGNS,VALUES,ERRORS,CENTRS,
2288  + HWIDTS,GREATE,ERROR,VALUE,CENTER,HWIDTH,DIR)
2289 C***BEGIN PROLOGUE DTRHRE
2290 C***PURPOSE DTRHRE maintains a heap of subregions.
2291 C***DESCRIPTION DTRHRE maintains a heap of subregions.
2292 C The subregions are ordered according to the size
2293 C of the greatest error estimates of each subregion(GREATE).
2294 C
2295 C PARAMETERS
2296 C
2297 C DVFLAG Integer.
2298 C If DVFLAG = 1, we remove the subregion with
2299 C greatest error from the heap.
2300 C If DVFLAG = 2, we insert a new subregion in the heap.
2301 C NDIM Integer.
2302 C Number of variables.
2303 C NUMFUN Integer.
2304 C Number of components of the integral.
2305 C SBRGNS Integer.
2306 C Number of subregions in the heap.
2307 C VALUES Real array of dimension (NUMFUN,SBRGNS).
2308 C Used to store estimated values of the integrals
2309 C over the subregions.
2310 C ERRORS Real array of dimension (NUMFUN,SBRGNS).
2311 C Used to store the corresponding estimated errors.
2312 C CENTRS Real array of dimension (NDIM,SBRGNS).
2313 C Used to store the center limits of the stored subregions.
2314 C HWIDTS Real array of dimension (NDIM,SBRGNS).
2315 C Used to store the hwidth limits of the stored subregions.
2316 C GREATE Real array of dimension SBRGNS.
2317 C Used to store the greatest estimated errors in
2318 C all subregions.
2319 C ERROR Real array of dimension NUMFUN.
2320 C Used as intermediate storage for the error of a subregion.
2321 C VALUE Real array of dimension NUMFUN.
2322 C Used as intermediate storage for the estimate
2323 C of the integral over a subregion.
2324 C CENTER Real array of dimension NDIM.
2325 C Used as intermediate storage for the center of
2326 C the subregion.
2327 C HWIDTH Real array of dimension NDIM.
2328 C Used as intermediate storage for the half width of
2329 C the subregion.
2330 C DIR Integer array of dimension SBRGNS.
2331 C DIR is used to store the directions for
2332 C further subdivision.
2333 C
2334 C***ROUTINES CALLED-NONE
2335 C***END PROLOGUE DTRHRE
2336 C
2337 C Global variables.
2338 C
2339  INTEGER DVFLAG,NDIM,NUMFUN,SBRGNS
2340  DOUBLE PRECISION VALUES(NUMFUN,*),ERRORS(NUMFUN,*)
2341  DOUBLE PRECISION CENTRS(NDIM,*)
2342  DOUBLE PRECISION HWIDTS(NDIM,*)
2343  DOUBLE PRECISION GREATE(*)
2344  DOUBLE PRECISION ERROR(NUMFUN),VALUE(NUMFUN)
2345  DOUBLE PRECISION CENTER(NDIM),HWIDTH(NDIM)
2346  DOUBLE PRECISION DIR(*)
2347 C
2348 C Local variables.
2349 C
2350 C GREAT is used as intermediate storage for the greatest error of a
2351 C subregion.
2352 C DIRECT is used as intermediate storage for the direction of further
2353 C subdivision.
2354 C SUBRGN Position of child/parent subregion in the heap.
2355 C SUBTMP Position of parent/child subregion in the heap.
2356 C
2357  INTEGER J,SUBRGN,SUBTMP
2358  DOUBLE PRECISION GREAT,DIRECT
2359 C
2360 C***FIRST EXECUTABLE STATEMENT DTRHRE
2361 C
2362 C Save values to be stored in their correct place in the heap.
2363 C
2364  GREAT = greate(sbrgns)
2365  direct = dir(sbrgns)
2366  DO 5 j = 1,numfun
2367  error(j) = errors(j,sbrgns)
2368  value(j) = values(j,sbrgns)
2369 5 CONTINUE
2370  DO 10 j = 1,ndim
2371  center(j) = centrs(j,sbrgns)
2372  hwidth(j) = hwidts(j,sbrgns)
2373 10 CONTINUE
2374 C
2375 C If DVFLAG = 1, we will remove the region
2376 C with greatest estimated error from the heap.
2377 C
2378  IF (dvflag.EQ.1) THEN
2379  sbrgns = sbrgns - 1
2380  subrgn = 1
2381 20 subtmp = 2*subrgn
2382  IF (subtmp.LE.sbrgns) THEN
2383  IF (subtmp.NE.sbrgns) THEN
2384 C
2385 C Find max. of left and right child.
2386 C
2387  IF (greate(subtmp).LT.greate(subtmp+1)) THEN
2388  subtmp = subtmp + 1
2389  END IF
2390  END IF
2391 C
2392 C Compare max.child with parent.
2393 C If parent is max., then done.
2394 C
2395  IF (great.LT.greate(subtmp)) THEN
2396 C
2397 C Move the values at position subtmp up the heap.
2398 C
2399  greate(subrgn) = greate(subtmp)
2400  DO 25 j = 1,numfun
2401  errors(j,subrgn) = errors(j,subtmp)
2402  values(j,subrgn) = values(j,subtmp)
2403 25 CONTINUE
2404  dir(subrgn) = dir(subtmp)
2405  DO 30 j = 1,ndim
2406  centrs(j,subrgn) = centrs(j,subtmp)
2407  hwidts(j,subrgn) = hwidts(j,subtmp)
2408 30 CONTINUE
2409  subrgn = subtmp
2410  GO TO 20
2411  END IF
2412  END IF
2413  ELSE IF (dvflag.EQ.2) THEN
2414 C
2415 C If DVFLAG = 2, then insert new region in the heap.
2416 C
2417  subrgn = sbrgns
2418 40 subtmp = subrgn/2
2419  IF (subtmp.GE.1) THEN
2420 C
2421 C Compare max.child with parent.
2422 C If parent is max, then done.
2423 C
2424  IF (great.GT.greate(subtmp)) THEN
2425 C
2426 C Move the values at position subtmp down the heap.
2427 C
2428  greate(subrgn) = greate(subtmp)
2429  DO 45 j = 1,numfun
2430  errors(j,subrgn) = errors(j,subtmp)
2431  values(j,subrgn) = values(j,subtmp)
2432 45 CONTINUE
2433  dir(subrgn) = dir(subtmp)
2434  DO 50 j = 1,ndim
2435  centrs(j,subrgn) = centrs(j,subtmp)
2436  hwidts(j,subrgn) = hwidts(j,subtmp)
2437 50 CONTINUE
2438  subrgn = subtmp
2439  GO TO 40
2440  END IF
2441  END IF
2442  END IF
2443 C
2444 C Insert the saved values in their correct places.
2445 C
2446  IF (sbrgns.GT.0) THEN
2447  greate(subrgn) = great
2448  DO 55 j = 1,numfun
2449  errors(j,subrgn) = error(j)
2450  values(j,subrgn) = value(j)
2451 55 CONTINUE
2452  dir(subrgn) = direct
2453  DO 60 j = 1,ndim
2454  centrs(j,subrgn) = center(j)
2455  hwidts(j,subrgn) = hwidth(j)
2456 60 CONTINUE
2457  END IF
2458 C
2459 C***END DTRHRE
2460 C
2461  RETURN
2462  END
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17