qawoe.f 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547
  1. *DECK QAWOE
  2. SUBROUTINE QAWOE (F, A, B, OMEGA, INTEGR, EPSABS, EPSREL, LIMIT,
  3. + ICALL, MAXP1, RESULT, ABSERR, NEVAL, IER, LAST, ALIST, BLIST,
  4. + RLIST, ELIST, IORD, NNLOG, MOMCOM, CHEBMO)
  5. C***BEGIN PROLOGUE QAWOE
  6. C***PURPOSE Calculate an approximation to a given definite integral
  7. C I = Integral of F(X)*W(X) over (A,B), where
  8. C W(X) = COS(OMEGA*X)
  9. C or W(X) = SIN(OMEGA*X),
  10. C hopefully satisfying the following claim for accuracy
  11. C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  12. C***LIBRARY SLATEC (QUADPACK)
  13. C***CATEGORY H2A2A1
  14. C***TYPE SINGLE PRECISION (QAWOE-S, DQAWOE-D)
  15. C***KEYWORDS AUTOMATIC INTEGRATOR, CLENSHAW-CURTIS METHOD,
  16. C EXTRAPOLATION, GLOBALLY ADAPTIVE,
  17. C INTEGRAND WITH OSCILLATORY COS OR SIN FACTOR, QUADPACK,
  18. C QUADRATURE, SPECIAL-PURPOSE
  19. C***AUTHOR Piessens, Robert
  20. C Applied Mathematics and Programming Division
  21. C K. U. Leuven
  22. C de Doncker, Elise
  23. C Applied Mathematics and Programming Division
  24. C K. U. Leuven
  25. C***DESCRIPTION
  26. C
  27. C Computation of Oscillatory integrals
  28. C Standard fortran subroutine
  29. C Real version
  30. C
  31. C PARAMETERS
  32. C ON ENTRY
  33. C F - Real
  34. C Function subprogram defining the integrand
  35. C function F(X). The actual name for F needs to be
  36. C declared E X T E R N A L in the driver program.
  37. C
  38. C A - Real
  39. C Lower limit of integration
  40. C
  41. C B - Real
  42. C Upper limit of integration
  43. C
  44. C OMEGA - Real
  45. C Parameter in the integrand weight function
  46. C
  47. C INTEGR - Integer
  48. C Indicates which of the WEIGHT functions is to be
  49. C used
  50. C INTEGR = 1 W(X) = COS(OMEGA*X)
  51. C INTEGR = 2 W(X) = SIN(OMEGA*X)
  52. C If INTEGR.NE.1 and INTEGR.NE.2, the routine
  53. C will end with IER = 6.
  54. C
  55. C EPSABS - Real
  56. C Absolute accuracy requested
  57. C EPSREL - Real
  58. C Relative accuracy requested
  59. C If EPSABS.LE.0
  60. C and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  61. C the routine will end with IER = 6.
  62. C
  63. C LIMIT - Integer
  64. C Gives an upper bound on the number of subdivisions
  65. C in the partition of (A,B), LIMIT.GE.1.
  66. C
  67. C ICALL - Integer
  68. C If QAWOE is to be used only once, ICALL must
  69. C be set to 1. Assume that during this call, the
  70. C Chebyshev moments (for CLENSHAW-CURTIS integration
  71. C of degree 24) have been computed for intervals of
  72. C lengths (ABS(B-A))*2**(-L), L=0,1,2,...MOMCOM-1.
  73. C If ICALL.GT.1 this means that QAWOE has been
  74. C called twice or more on intervals of the same
  75. C length ABS(B-A). The Chebyshev moments already
  76. C computed are then re-used in subsequent calls.
  77. C If ICALL.LT.1, the routine will end with IER = 6.
  78. C
  79. C MAXP1 - Integer
  80. C Gives an upper bound on the number of Chebyshev
  81. C moments which can be stored, i.e. for the
  82. C intervals of lengths ABS(B-A)*2**(-L),
  83. C L=0,1, ..., MAXP1-2, MAXP1.GE.1.
  84. C If MAXP1.LT.1, the routine will end with IER = 6.
  85. C
  86. C ON RETURN
  87. C RESULT - Real
  88. C Approximation to the integral
  89. C
  90. C ABSERR - Real
  91. C Estimate of the modulus of the absolute error,
  92. C which should equal or exceed ABS(I-RESULT)
  93. C
  94. C NEVAL - Integer
  95. C Number of integrand evaluations
  96. C
  97. C IER - Integer
  98. C IER = 0 Normal and reliable termination of the
  99. C routine. It is assumed that the
  100. C requested accuracy has been achieved.
  101. C - IER.GT.0 Abnormal termination of the routine.
  102. C The estimates for integral and error are
  103. C less reliable. It is assumed that the
  104. C requested accuracy has not been achieved.
  105. C ERROR MESSAGES
  106. C IER = 1 Maximum number of subdivisions allowed
  107. C has been achieved. One can allow more
  108. C subdivisions by increasing the value of
  109. C LIMIT (and taking according dimension
  110. C adjustments into account). However, if
  111. C this yields no improvement it is advised
  112. C to analyze the integrand, in order to
  113. C determine the integration difficulties.
  114. C If the position of a local difficulty can
  115. C be determined (e.g. SINGULARITY,
  116. C DISCONTINUITY within the interval) one
  117. C will probably gain from splitting up the
  118. C interval at this point and calling the
  119. C integrator on the subranges. If possible,
  120. C an appropriate special-purpose integrator
  121. C should be used which is designed for
  122. C handling the type of difficulty involved.
  123. C = 2 The occurrence of roundoff error is
  124. C detected, which prevents the requested
  125. C tolerance from being achieved.
  126. C The error may be under-estimated.
  127. C = 3 Extremely bad integrand behaviour occurs
  128. C at some points of the integration
  129. C interval.
  130. C = 4 The algorithm does not converge.
  131. C Roundoff error is detected in the
  132. C extrapolation table.
  133. C It is presumed that the requested
  134. C tolerance cannot be achieved due to
  135. C roundoff in the extrapolation table,
  136. C and that the returned result is the
  137. C best which can be obtained.
  138. C = 5 The integral is probably divergent, or
  139. C slowly convergent. It must be noted that
  140. C divergence can occur with any other value
  141. C of IER.GT.0.
  142. C = 6 The input is invalid, because
  143. C (EPSABS.LE.0 and
  144. C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
  145. C or (INTEGR.NE.1 and INTEGR.NE.2) or
  146. C ICALL.LT.1 or MAXP1.LT.1.
  147. C RESULT, ABSERR, NEVAL, LAST, RLIST(1),
  148. C ELIST(1), IORD(1) and NNLOG(1) are set
  149. C to ZERO. ALIST(1) and BLIST(1) are set
  150. C to A and B respectively.
  151. C
  152. C LAST - Integer
  153. C On return, LAST equals the number of
  154. C subintervals produces in the subdivision
  155. C process, which determines the number of
  156. C significant elements actually in the
  157. C WORK ARRAYS.
  158. C ALIST - Real
  159. C Vector of dimension at least LIMIT, the first
  160. C LAST elements of which are the left
  161. C end points of the subintervals in the partition
  162. C of the given integration range (A,B)
  163. C
  164. C BLIST - Real
  165. C Vector of dimension at least LIMIT, the first
  166. C LAST elements of which are the right
  167. C end points of the subintervals in the partition
  168. C of the given integration range (A,B)
  169. C
  170. C RLIST - Real
  171. C Vector of dimension at least LIMIT, the first
  172. C LAST elements of which are the integral
  173. C approximations on the subintervals
  174. C
  175. C ELIST - Real
  176. C Vector of dimension at least LIMIT, the first
  177. C LAST elements of which are the moduli of the
  178. C absolute error estimates on the subintervals
  179. C
  180. C IORD - Integer
  181. C Vector of dimension at least LIMIT, the first K
  182. C elements of which are pointers to the error
  183. C estimates over the subintervals,
  184. C such that ELIST(IORD(1)), ...,
  185. C ELIST(IORD(K)) form a decreasing sequence, with
  186. C K = LAST if LAST.LE.(LIMIT/2+2), and
  187. C K = LIMIT+1-LAST otherwise.
  188. C
  189. C NNLOG - Integer
  190. C Vector of dimension at least LIMIT, containing the
  191. C subdivision levels of the subintervals, i.e.
  192. C IWORK(I) = L means that the subinterval
  193. C numbered I is of length ABS(B-A)*2**(1-L)
  194. C
  195. C ON ENTRY AND RETURN
  196. C MOMCOM - Integer
  197. C Indicating that the Chebyshev moments
  198. C have been computed for intervals of lengths
  199. C (ABS(B-A))*2**(-L), L=0,1,2, ..., MOMCOM-1,
  200. C MOMCOM.LT.MAXP1
  201. C
  202. C CHEBMO - Real
  203. C Array of dimension (MAXP1,25) containing the
  204. C Chebyshev moments
  205. C
  206. C***REFERENCES (NONE)
  207. C***ROUTINES CALLED QC25F, QELG, QPSRT, R1MACH
  208. C***REVISION HISTORY (YYMMDD)
  209. C 800101 DATE WRITTEN
  210. C 890531 Changed all specific intrinsics to generic. (WRB)
  211. C 890831 Modified array declarations. (WRB)
  212. C 890831 REVISION DATE from Version 3.2
  213. C 891214 Prologue converted to Version 4.0 format. (BAB)
  214. C***END PROLOGUE QAWOE
  215. C
  216. REAL A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
  217. 1 A2,B,BLIST,B1,B2,CHEBMO,CORREC,DEFAB1,DEFAB2,DEFABS,
  218. 2 DOMEGA,R1MACH,DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,
  219. 3 ERLAST,ERRBND,ERRMAX,ERROR1,ERRO12,ERROR2,ERRSUM,ERTEST,F,OFLOW,
  220. 4 OMEGA,RESABS,RESEPS,RESULT,RES3LA,RLIST,RLIST2,SMALL,UFLOW,WIDTH
  221. INTEGER ICALL,ID,IER,IERRO,INTEGR,IORD,IROFF1,IROFF2,IROFF3,
  222. 1 JUPBND,K,KSGN,KTMIN,LAST,LIMIT,MAXERR,MAXP1,MOMCOM,NEV,
  223. 2 NEVAL,NNLOG,NRES,NRMAX,NRMOM,NUMRL2
  224. LOGICAL EXTRAP,NOEXT,EXTALL
  225. C
  226. DIMENSION ALIST(*),BLIST(*),RLIST(*),ELIST(*),
  227. 1 IORD(*),RLIST2(52),RES3LA(3),CHEBMO(MAXP1,25),NNLOG(*)
  228. C
  229. EXTERNAL F
  230. C
  231. C THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
  232. C LIMEXP IN SUBROUTINE QELG (RLIST2 SHOULD BE OF
  233. C DIMENSION (LIMEXP+2) AT LEAST).
  234. C
  235. C LIST OF MAJOR VARIABLES
  236. C -----------------------
  237. C
  238. C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
  239. C CONSIDERED UP TO NOW
  240. C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
  241. C CONSIDERED UP TO NOW
  242. C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER
  243. C (ALIST(I),BLIST(I))
  244. C RLIST2 - ARRAY OF DIMENSION AT LEAST LIMEXP+2
  245. C CONTAINING THE PART OF THE EPSILON TABLE
  246. C WHICH IS STILL NEEDED FOR FURTHER COMPUTATIONS
  247. C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I)
  248. C MAXERR - POINTER TO THE INTERVAL WITH LARGEST
  249. C ERROR ESTIMATE
  250. C ERRMAX - ELIST(MAXERR)
  251. C ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
  252. C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
  253. C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS
  254. C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
  255. C ABS(RESULT))
  256. C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL
  257. C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL
  258. C LAST - INDEX FOR SUBDIVISION
  259. C NRES - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
  260. C NUMRL2 - NUMBER OF ELEMENTS IN RLIST2. IF AN APPROPRIATE
  261. C APPROXIMATION TO THE COMPOUNDED INTEGRAL HAS
  262. C BEEN OBTAINED IT IS PUT IN RLIST2(NUMRL2) AFTER
  263. C NUMRL2 HAS BEEN INCREASED BY ONE
  264. C SMALL - LENGTH OF THE SMALLEST INTERVAL CONSIDERED
  265. C UP TO NOW, MULTIPLIED BY 1.5
  266. C ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER
  267. C THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
  268. C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE IS
  269. C ATTEMPTING TO PERFORM EXTRAPOLATION, I.E. BEFORE
  270. C SUBDIVIDING THE SMALLEST INTERVAL WE TRY TO
  271. C DECREASE THE VALUE OF ERLARG
  272. C NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION
  273. C IS NO LONGER ALLOWED (TRUE VALUE)
  274. C
  275. C MACHINE DEPENDENT CONSTANTS
  276. C ---------------------------
  277. C
  278. C EPMACH IS THE LARGEST RELATIVE SPACING.
  279. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  280. C OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  281. C
  282. C***FIRST EXECUTABLE STATEMENT QAWOE
  283. EPMACH = R1MACH(4)
  284. C
  285. C TEST ON VALIDITY OF PARAMETERS
  286. C ------------------------------
  287. C
  288. IER = 0
  289. NEVAL = 0
  290. LAST = 0
  291. RESULT = 0.0E+00
  292. ABSERR = 0.0E+00
  293. ALIST(1) = A
  294. BLIST(1) = B
  295. RLIST(1) = 0.0E+00
  296. ELIST(1) = 0.0E+00
  297. IORD(1) = 0
  298. NNLOG(1) = 0
  299. IF((INTEGR.NE.1.AND.INTEGR.NE.2).OR.(EPSABS.LE.0.0E+00.AND.
  300. 1 EPSREL.LT.MAX(0.5E+02*EPMACH,0.5E-14)).OR.ICALL.LT.1.OR.
  301. 2 MAXP1.LT.1) IER = 6
  302. IF(IER.EQ.6) GO TO 999
  303. C
  304. C FIRST APPROXIMATION TO THE INTEGRAL
  305. C -----------------------------------
  306. C
  307. DOMEGA = ABS(OMEGA)
  308. NRMOM = 0
  309. IF (ICALL.GT.1) GO TO 5
  310. MOMCOM = 0
  311. 5 CALL QC25F(F,A,B,DOMEGA,INTEGR,NRMOM,MAXP1,0,RESULT,ABSERR,
  312. 1 NEVAL,DEFABS,RESABS,MOMCOM,CHEBMO)
  313. C
  314. C TEST ON ACCURACY.
  315. C
  316. DRES = ABS(RESULT)
  317. ERRBND = MAX(EPSABS,EPSREL*DRES)
  318. RLIST(1) = RESULT
  319. ELIST(1) = ABSERR
  320. IORD(1) = 1
  321. IF(ABSERR.LE.0.1E+03*EPMACH*DEFABS.AND.ABSERR.GT.
  322. 1 ERRBND) IER = 2
  323. IF(LIMIT.EQ.1) IER = 1
  324. IF(IER.NE.0.OR.ABSERR.LE.ERRBND) GO TO 200
  325. C
  326. C INITIALIZATIONS
  327. C ---------------
  328. C
  329. UFLOW = R1MACH(1)
  330. OFLOW = R1MACH(2)
  331. ERRMAX = ABSERR
  332. MAXERR = 1
  333. AREA = RESULT
  334. ERRSUM = ABSERR
  335. ABSERR = OFLOW
  336. NRMAX = 1
  337. EXTRAP = .FALSE.
  338. NOEXT = .FALSE.
  339. IERRO = 0
  340. IROFF1 = 0
  341. IROFF2 = 0
  342. IROFF3 = 0
  343. KTMIN = 0
  344. SMALL = ABS(B-A)*0.75E+00
  345. NRES = 0
  346. NUMRL2 = 0
  347. EXTALL = .FALSE.
  348. IF(0.5E+00*ABS(B-A)*DOMEGA.GT.0.2E+01) GO TO 10
  349. NUMRL2 = 1
  350. EXTALL = .TRUE.
  351. RLIST2(1) = RESULT
  352. 10 IF(0.25E+00*ABS(B-A)*DOMEGA.LE.0.2E+01) EXTALL = .TRUE.
  353. KSGN = -1
  354. IF(DRES.GE.(0.1E+01-0.5E+02*EPMACH)*DEFABS) KSGN = 1
  355. C
  356. C MAIN DO-LOOP
  357. C ------------
  358. C
  359. DO 140 LAST = 2,LIMIT
  360. C
  361. C BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST
  362. C ERROR ESTIMATE.
  363. C
  364. NRMOM = NNLOG(MAXERR)+1
  365. A1 = ALIST(MAXERR)
  366. B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR))
  367. A2 = B1
  368. B2 = BLIST(MAXERR)
  369. ERLAST = ERRMAX
  370. CALL QC25F(F,A1,B1,DOMEGA,INTEGR,NRMOM,MAXP1,0,
  371. 1 AREA1,ERROR1,NEV,RESABS,DEFAB1,MOMCOM,CHEBMO)
  372. NEVAL = NEVAL+NEV
  373. CALL QC25F(F,A2,B2,DOMEGA,INTEGR,NRMOM,MAXP1,1,
  374. 1 AREA2,ERROR2,NEV,RESABS,DEFAB2,MOMCOM,CHEBMO)
  375. NEVAL = NEVAL+NEV
  376. C
  377. C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
  378. C AND ERROR AND TEST FOR ACCURACY.
  379. C
  380. AREA12 = AREA1+AREA2
  381. ERRO12 = ERROR1+ERROR2
  382. ERRSUM = ERRSUM+ERRO12-ERRMAX
  383. AREA = AREA+AREA12-RLIST(MAXERR)
  384. IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 25
  385. IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1E-04*ABS(AREA12)
  386. 1 .OR.ERRO12.LT.0.99E+00*ERRMAX) GO TO 20
  387. IF(EXTRAP) IROFF2 = IROFF2+1
  388. IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
  389. 20 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
  390. 25 RLIST(MAXERR) = AREA1
  391. RLIST(LAST) = AREA2
  392. NNLOG(MAXERR) = NRMOM
  393. NNLOG(LAST) = NRMOM
  394. ERRBND = MAX(EPSABS,EPSREL*ABS(AREA))
  395. C
  396. C TEST FOR ROUNDOFF ERROR AND EVENTUALLY
  397. C SET ERROR FLAG
  398. C
  399. IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
  400. IF(IROFF2.GE.5) IERRO = 3
  401. C
  402. C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
  403. C SUBINTERVALS EQUALS LIMIT.
  404. C
  405. IF(LAST.EQ.LIMIT) IER = 1
  406. C
  407. C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
  408. C AT A POINT OF THE INTEGRATION RANGE.
  409. C
  410. IF(MAX(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03*EPMACH)
  411. 1 *(ABS(A2)+0.1E+04*UFLOW)) IER = 4
  412. C
  413. C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  414. C
  415. IF(ERROR2.GT.ERROR1) GO TO 30
  416. ALIST(LAST) = A2
  417. BLIST(MAXERR) = B1
  418. BLIST(LAST) = B2
  419. ELIST(MAXERR) = ERROR1
  420. ELIST(LAST) = ERROR2
  421. GO TO 40
  422. 30 ALIST(MAXERR) = A2
  423. ALIST(LAST) = A1
  424. BLIST(LAST) = B1
  425. RLIST(MAXERR) = AREA2
  426. RLIST(LAST) = AREA1
  427. ELIST(MAXERR) = ERROR2
  428. ELIST(LAST) = ERROR1
  429. C
  430. C CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING
  431. C IN THE LIST OF ERROR ESTIMATES AND SELECT THE
  432. C SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE
  433. C BISECTED NEXT).
  434. C
  435. 40 CALL QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
  436. C ***JUMP OUT OF DO-LOOP
  437. IF(ERRSUM.LE.ERRBND) GO TO 170
  438. IF(IER.NE.0) GO TO 150
  439. IF(LAST.EQ.2.AND.EXTALL) GO TO 120
  440. IF(NOEXT) GO TO 140
  441. IF(.NOT.EXTALL) GO TO 50
  442. ERLARG = ERLARG-ERLAST
  443. IF(ABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
  444. IF(EXTRAP) GO TO 70
  445. C
  446. C TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
  447. C SMALLEST INTERVAL.
  448. C
  449. 50 WIDTH = ABS(BLIST(MAXERR)-ALIST(MAXERR))
  450. IF(WIDTH.GT.SMALL) GO TO 140
  451. IF(EXTALL) GO TO 60
  452. C
  453. C TEST WHETHER WE CAN START WITH THE EXTRAPOLATION
  454. C PROCEDURE (WE DO THIS IF WE INTEGRATE OVER THE
  455. C NEXT INTERVAL WITH USE OF A GAUSS-KRONROD RULE - SEE
  456. C SUBROUTINE QC25F).
  457. C
  458. SMALL = SMALL*0.5E+00
  459. IF(0.25E+00*WIDTH*DOMEGA.GT.0.2E+01) GO TO 140
  460. EXTALL = .TRUE.
  461. GO TO 130
  462. 60 EXTRAP = .TRUE.
  463. NRMAX = 2
  464. 70 IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 90
  465. C
  466. C THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
  467. C BEFORE BISECTING DECREASE THE SUM OF THE ERRORS
  468. C OVER THE LARGER INTERVALS (ERLARG) AND PERFORM
  469. C EXTRAPOLATION.
  470. C
  471. JUPBND = LAST
  472. IF (LAST.GT.(LIMIT/2+2)) JUPBND = LIMIT+3-LAST
  473. ID = NRMAX
  474. DO 80 K = ID,JUPBND
  475. MAXERR = IORD(NRMAX)
  476. ERRMAX = ELIST(MAXERR)
  477. IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 140
  478. NRMAX = NRMAX+1
  479. 80 CONTINUE
  480. C
  481. C PERFORM EXTRAPOLATION.
  482. C
  483. 90 NUMRL2 = NUMRL2+1
  484. RLIST2(NUMRL2) = AREA
  485. IF(NUMRL2.LT.3) GO TO 110
  486. CALL QELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
  487. KTMIN = KTMIN+1
  488. IF(KTMIN.GT.5.AND.ABSERR.LT.0.1E-02*ERRSUM) IER = 5
  489. IF(ABSEPS.GE.ABSERR) GO TO 100
  490. KTMIN = 0
  491. ABSERR = ABSEPS
  492. RESULT = RESEPS
  493. CORREC = ERLARG
  494. ERTEST = MAX(EPSABS,EPSREL*ABS(RESEPS))
  495. C ***JUMP OUT OF DO-LOOP
  496. IF(ABSERR.LE.ERTEST) GO TO 150
  497. C
  498. C PREPARE BISECTION OF THE SMALLEST INTERVAL.
  499. C
  500. 100 IF(NUMRL2.EQ.1) NOEXT = .TRUE.
  501. IF(IER.EQ.5) GO TO 150
  502. 110 MAXERR = IORD(1)
  503. ERRMAX = ELIST(MAXERR)
  504. NRMAX = 1
  505. EXTRAP = .FALSE.
  506. SMALL = SMALL*0.5E+00
  507. ERLARG = ERRSUM
  508. GO TO 140
  509. 120 SMALL = SMALL*0.5E+00
  510. NUMRL2 = NUMRL2+1
  511. RLIST2(NUMRL2) = AREA
  512. 130 ERTEST = ERRBND
  513. ERLARG = ERRSUM
  514. 140 CONTINUE
  515. C
  516. C SET THE FINAL RESULT.
  517. C ---------------------
  518. C
  519. 150 IF(ABSERR.EQ.OFLOW.OR.NRES.EQ.0) GO TO 170
  520. IF(IER+IERRO.EQ.0) GO TO 165
  521. IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
  522. IF(IER.EQ.0) IER = 3
  523. IF(RESULT.NE.0.0E+00.AND.AREA.NE.0.0E+00) GO TO 160
  524. IF(ABSERR.GT.ERRSUM) GO TO 170
  525. IF(AREA.EQ.0.0E+00) GO TO 190
  526. GO TO 165
  527. 160 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA)) GO TO 170
  528. C
  529. C TEST ON DIVERGENCE.
  530. C
  531. 165 IF(KSGN.EQ.(-1).AND.MAX(ABS(RESULT),ABS(AREA)).LE.
  532. 1 DEFABS*0.1E-01) GO TO 190
  533. IF(0.1E-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1E+03
  534. 1 .OR.ERRSUM.GE.ABS(AREA)) IER = 6
  535. GO TO 190
  536. C
  537. C COMPUTE GLOBAL INTEGRAL SUM.
  538. C
  539. 170 RESULT = 0.0E+00
  540. DO 180 K=1,LAST
  541. RESULT = RESULT+RLIST(K)
  542. 180 CONTINUE
  543. ABSERR = ERRSUM
  544. 190 IF (IER.GT.2) IER=IER-1
  545. 200 IF (INTEGR.EQ.2.AND.OMEGA.LT.0.0E+00) RESULT=-RESULT
  546. 999 RETURN
  547. END