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