relibc_openlibm/slatec/davint.f

*DECK DAVINT
      SUBROUTINE DAVINT (X, Y, N, XLO, XUP, ANS, IERR)
C***BEGIN PROLOGUE  DAVINT
C***PURPOSE  Integrate a function tabulated at arbitrarily spaced
C            abscissas using overlapping parabolas.
C***LIBRARY   SLATEC
C***CATEGORY  H2A1B2
C***TYPE      DOUBLE PRECISION (AVINT-S, DAVINT-D)
C***KEYWORDS  INTEGRATION, QUADRATURE, TABULATED DATA
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C         DAVINT integrates a function tabulated at arbitrarily spaced
C         abscissas.  The limits of integration need not coincide
C         with the tabulated abscissas.
C
C         A method of overlapping parabolas fitted to the data is used
C         provided that there are at least 3 abscissas between the
C         limits of integration.  DAVINT also handles two special cases.
C         If the limits of integration are equal, DAVINT returns a
C         result of zero regardless of the number of tabulated values.
C         If there are only two function values, DAVINT uses the
C         trapezoid rule.
C
C     Description of Parameters
C         The user must dimension all arrays appearing in the call list
C              X(N), Y(N)
C
C         Input--
C      X    - DOUBLE PRECISION array of abscissas, which must be in
C             increasing order.
C      Y    - DOUBLE PRECISION array of function values. i.e.,
C                Y(I)=FUNC(X(I))
C      N    - The integer number of function values supplied.
C                N .GE. 2 unless XLO = XUP.
C      XLO  - DOUBLE PRECISION lower limit of integration
C      XUP  - DOUBLE PRECISION upper limit of integration.  Must have
C              XLO.LE.XUP
C
C         Output--
C      ANS  - Double Precision computed approximate value of integral
C      IERR - A status code
C           --Normal Code
C                =1 Means the requested integration was performed.
C           --Abnormal Codes
C                =2 Means XUP was less than XLO.
C                =3 Means the number of X(I) between XLO and XUP
C                   (inclusive) was less than 3 and neither of the two
C                   special cases described in the abstract occurred.
C                   No integration was performed.
C                =4 Means the restriction X(I+1).GT.X(I) was violated.
C                =5 Means the number N of function values was .lt. 2.
C                   ANS is set to zero if IERR=2,3,4,or 5.
C
C    DAVINT is documented completely in SC-M-69-335
C    Original program from *Numerical Integration* by Davis & Rabinowitz
C    Adaptation and modifications by Rondall E Jones.
C
C***REFERENCES  R. E. Jones, Approximate integrator of functions
C                 tabulated at arbitrarily spaced abscissas,
C                 Report SC-M-69-335, Sandia Laboratories, 1969.
C***ROUTINES CALLED  XERMSG
C***REVISION HISTORY  (YYMMDD)
C   690901  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DAVINT
C
      INTEGER I, IERR, INLFT, INRT, ISTART, ISTOP, N
      DOUBLE PRECISION A, ANS, B, C, CA, CB, CC, FL, FR, R3, RP5,
     1     SLOPE, SUM, SYL, SYL2, SYL3, SYU, SYU2, SYU3, TERM1, TERM2,
     2     TERM3, X, X1, X12, X13, X2, X23, X3, XLO, XUP, Y
      DIMENSION X(*),Y(*)
C     BEGIN BLOCK PERMITTING ...EXITS TO 190
C        BEGIN BLOCK PERMITTING ...EXITS TO 180
C***FIRST EXECUTABLE STATEMENT  DAVINT
            IERR = 1
            ANS = 0.0D0
            IF (XLO .GT. XUP) GO TO 160
               IF (XLO .EQ. XUP) GO TO 150
                  IF (N .GE. 2) GO TO 10
                     IERR = 5
                     CALL XERMSG ('SLATEC', 'DAVINT',
     +                  'LESS THAN TWO FUNCTION VALUES WERE SUPPLIED.',
     +                  4, 1)
C     ...............EXIT
                     GO TO 190
   10             CONTINUE
                  DO 20 I = 2, N
C        ............EXIT
                     IF (X(I) .LE. X(I-1)) GO TO 180
C                 ...EXIT
                     IF (X(I) .GT. XUP) GO TO 30
   20             CONTINUE
   30             CONTINUE
                  IF (N .GE. 3) GO TO 40
C
C                    SPECIAL N=2 CASE
                     SLOPE = (Y(2) - Y(1))/(X(2) - X(1))
                     FL = Y(1) + SLOPE*(XLO - X(1))
                     FR = Y(2) + SLOPE*(XUP - X(2))
                     ANS = 0.5D0*(FL + FR)*(XUP - XLO)
C     ...............EXIT
                     GO TO 190
   40             CONTINUE
                  IF (X(N-2) .GE. XLO) GO TO 50
                     IERR = 3
                     CALL XERMSG ('SLATEC', 'DAVINT',
     +                  'THERE WERE LESS THAN THREE FUNCTION VALUES ' //
     +                  'BETWEEN THE LIMITS OF INTEGRATION.', 4, 1)
C     ...............EXIT
                     GO TO 190
   50             CONTINUE
                  IF (X(3) .LE. XUP) GO TO 60
                     IERR = 3
                     CALL XERMSG ('SLATEC', 'DAVINT',
     +                  'THERE WERE LESS THAN THREE FUNCTION VALUES ' //
     +                  'BETWEEN THE LIMITS OF INTEGRATION.', 4, 1)
C     ...............EXIT
                     GO TO 190
   60             CONTINUE
                  I = 1
   70             IF (X(I) .GE. XLO) GO TO 80
                     I = I + 1
                  GO TO 70
   80             CONTINUE
                  INLFT = I
                  I = N
   90             IF (X(I) .LE. XUP) GO TO 100
                     I = I - 1
                  GO TO 90
  100             CONTINUE
                  INRT = I
                  IF ((INRT - INLFT) .GE. 2) GO TO 110
                     IERR = 3
                     CALL XERMSG ('SLATEC', 'DAVINT',
     +                  'THERE WERE LESS THAN THREE FUNCTION VALUES ' //
     +                  'BETWEEN THE LIMITS OF INTEGRATION.', 4, 1)
C     ...............EXIT
                     GO TO 190
  110             CONTINUE
                  ISTART = INLFT
                  IF (INLFT .EQ. 1) ISTART = 2
                  ISTOP = INRT
                  IF (INRT .EQ. N) ISTOP = N - 1
C
                  R3 = 3.0D0
                  RP5 = 0.5D0
                  SUM = 0.0D0
                  SYL = XLO
                  SYL2 = SYL*SYL
                  SYL3 = SYL2*SYL
C
                  DO 140 I = ISTART, ISTOP
                     X1 = X(I-1)
                     X2 = X(I)
                     X3 = X(I+1)
                     X12 = X1 - X2
                     X13 = X1 - X3
                     X23 = X2 - X3
                     TERM1 = Y(I-1)/(X12*X13)
                     TERM2 = -Y(I)/(X12*X23)
                     TERM3 = Y(I+1)/(X13*X23)
                     A = TERM1 + TERM2 + TERM3
                     B = -(X2 + X3)*TERM1 - (X1 + X3)*TERM2
     1                   - (X1 + X2)*TERM3
                     C = X2*X3*TERM1 + X1*X3*TERM2 + X1*X2*TERM3
                     IF (I .GT. ISTART) GO TO 120
                        CA = A
                        CB = B
                        CC = C
                     GO TO 130
  120                CONTINUE
                        CA = 0.5D0*(A + CA)
                        CB = 0.5D0*(B + CB)
                        CC = 0.5D0*(C + CC)
  130                CONTINUE
                     SYU = X2
                     SYU2 = SYU*SYU
                     SYU3 = SYU2*SYU
                     SUM = SUM + CA*(SYU3 - SYL3)/R3
     1                     + CB*RP5*(SYU2 - SYL2) + CC*(SYU - SYL)
                     CA = A
                     CB = B
                     CC = C
                     SYL = SYU
                     SYL2 = SYU2
                     SYL3 = SYU3
  140             CONTINUE
                  SYU = XUP
                  ANS = SUM + CA*(SYU**3 - SYL3)/R3
     1                  + CB*RP5*(SYU**2 - SYL2) + CC*(SYU - SYL)
  150          CONTINUE
            GO TO 170
  160       CONTINUE
               IERR = 2
               CALL XERMSG ('SLATEC', 'DAVINT',
     +            'THE UPPER LIMIT OF INTEGRATION WAS NOT GREATER ' //
     +            'THAN THE LOWER LIMIT.', 4, 1)
  170       CONTINUE
C     ......EXIT
            GO TO 190
  180    CONTINUE
         IERR = 4
         CALL XERMSG ('SLATEC', 'DAVINT',
     +      'THE ABSCISSAS WERE NOT STRICTLY INCREASING.  MUST HAVE ' //
     +      'X(I-1) .LT. X(I) FOR ALL I.', 4, 1)
  190 CONTINUE
      RETURN
      END