*DECK PCHIA
REAL FUNCTION PCHIA (N, X, F, D, INCFD, SKIP, A, B, IERR)
C***BEGIN PROLOGUE PCHIA
C***PURPOSE Evaluate the definite integral of a piecewise cubic
C Hermite function over an arbitrary interval.
C***LIBRARY SLATEC (PCHIP)
C***CATEGORY E3, H2A1B2
C***TYPE SINGLE PRECISION (PCHIA-S, DPCHIA-D)
C***KEYWORDS CUBIC HERMITE INTERPOLATION, NUMERICAL INTEGRATION, PCHIP,
C QUADRATURE
C***AUTHOR Fritsch, F. N., (LLNL)
C Lawrence Livermore National Laboratory
C P.O. Box 808 (L-316)
C Livermore, CA 94550
C FTS 532-4275, (510) 422-4275
C***DESCRIPTION
C
C PCHIA: Piecewise Cubic Hermite Integrator, Arbitrary limits
C
C Evaluates the definite integral of the cubic Hermite function
C defined by N, X, F, D over the interval [A, B].
C
C To provide compatibility with PCHIM and PCHIC, includes an
C increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C Calling sequence:
C
C PARAMETER (INCFD = ...)
C INTEGER N, IERR
C REAL X(N), F(INCFD,N), D(INCFD,N), A, B
C REAL VALUE, PCHIA
C LOGICAL SKIP
C
C VALUE = PCHIA (N, X, F, D, INCFD, SKIP, A, B, IERR)
C
C Parameters:
C
C VALUE -- (output) value of the requested integral.
C
C N -- (input) number of data points. (Error return if N.LT.2 .)
C
C X -- (input) real array of independent variable values. The
C elements of X must be strictly increasing:
C X(I-1) .LT. X(I), I = 2(1)N.
C (Error return if not.)
C
C F -- (input) real array of function values. F(1+(I-1)*INCFD) is
C the value corresponding to X(I).
C
C D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is
C the value corresponding to X(I).
C
C INCFD -- (input) increment between successive values in F and D.
C (Error return if INCFD.LT.1 .)
C
C SKIP -- (input/output) logical variable which should be set to
C .TRUE. if the user wishes to skip checks for validity of
C preceding parameters, or to .FALSE. otherwise.
C This will save time in case these checks have already
C been performed (say, in PCHIM or PCHIC).
C SKIP will be set to .TRUE. on return with IERR.GE.0 .
C
C A,B -- (input) the limits of integration.
C NOTE: There is no requirement that [A,B] be contained in
C [X(1),X(N)]. However, the resulting integral value
C will be highly suspect, if not.
C
C IERR -- (output) error flag.
C Normal return:
C IERR = 0 (no errors).
C Warning errors:
C IERR = 1 if A is outside the interval [X(1),X(N)].
C IERR = 2 if B is outside the interval [X(1),X(N)].
C IERR = 3 if both of the above are true. (Note that this
C means that either [A,B] contains data interval
C or the intervals do not intersect at all.)
C "Recoverable" errors:
C IERR = -1 if N.LT.2 .
C IERR = -2 if INCFD.LT.1 .
C IERR = -3 if the X-array is not strictly increasing.
C (VALUE will be zero in any of these cases.)
C NOTE: The above errors are checked in the order listed,
C and following arguments have **NOT** been validated.
C IERR = -4 in case of an error return from PCHID (which
C should never occur).
C
C***REFERENCES (NONE)
C***ROUTINES CALLED CHFIE, PCHID, XERMSG
C***REVISION HISTORY (YYMMDD)
C 820730 DATE WRITTEN
C 820804 Converted to SLATEC library version.
C 870707 Corrected double precision conversion instructions.
C 870813 Minor cosmetic changes.
C 890411 Added SAVE statements (Vers. 3.2).
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890703 Corrected category record. (WRB)
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 930503 Corrected to set VALUE=0 when IERR.lt.0. (FNF)
C 930504 Changed CHFIV to CHFIE. (FNF)
C***END PROLOGUE PCHIA
C
C Programming notes:
C 1. The error flag from PCHID is tested, because a logic flaw
C could conceivably result in IERD=-4, which should be reported.
C**End
C
C DECLARE ARGUMENTS.
C
INTEGER N, INCFD, IERR
REAL X(*), F(INCFD,*), D(INCFD,*), A, B
LOGICAL SKIP
C
C DECLARE LOCAL VARIABLES.
C
INTEGER I, IA, IB, IERD, IL, IR
REAL VALUE, XA, XB, ZERO
SAVE ZERO
REAL CHFIE, PCHID
C
C INITIALIZE.
C
DATA ZERO /0./
C***FIRST EXECUTABLE STATEMENT PCHIA
VALUE = ZERO
C
C VALIDITY-CHECK ARGUMENTS.
C
IF (SKIP) GO TO 5
C
IF ( N.LT.2 ) GO TO 5001
IF ( INCFD.LT.1 ) GO TO 5002
DO 1 I = 2, N
IF ( X(I).LE.X(I-1) ) GO TO 5003
1 CONTINUE
C
C FUNCTION DEFINITION IS OK, GO ON.
C
5 CONTINUE
SKIP = .TRUE.
IERR = 0
IF ( (A.LT.X(1)) .OR. (A.GT.X(N)) ) IERR = IERR + 1
IF ( (B.LT.X(1)) .OR. (B.GT.X(N)) ) IERR = IERR + 2
C
C COMPUTE INTEGRAL VALUE.
C
IF (A .NE. B) THEN
XA = MIN (A, B)
XB = MAX (A, B)
IF (XB .LE. X(2)) THEN
C INTERVAL IS TO LEFT OF X(2), SO USE FIRST CUBIC.
C --------------------------------------
VALUE = CHFIE (X(1),X(2), F(1,1),F(1,2),
+ D(1,1),D(1,2), A, B)
C --------------------------------------
ELSE IF (XA .GE. X(N-1)) THEN
C INTERVAL IS TO RIGHT OF X(N-1), SO USE LAST CUBIC.
C -----------------------------------------
VALUE = CHFIE(X(N-1),X(N), F(1,N-1),F(1,N),
+ D(1,N-1),D(1,N), A, B)
C -----------------------------------------
ELSE
C 'NORMAL' CASE -- XA.LT.XB, XA.LT.X(N-1), XB.GT.X(2).
C ......LOCATE IA AND IB SUCH THAT
C X(IA-1).LT.XA.LE.X(IA).LE.X(IB).LE.XB.LE.X(IB+1)
IA = 1
DO 10 I = 1, N-1
IF (XA .GT. X(I)) IA = I + 1
10 CONTINUE
C IA = 1 IMPLIES XA.LT.X(1) . OTHERWISE,
C IA IS LARGEST INDEX SUCH THAT X(IA-1).LT.XA,.
C
IB = N
DO 20 I = N, IA, -1
IF (XB .LT. X(I)) IB = I - 1
20 CONTINUE
C IB = N IMPLIES XB.GT.X(N) . OTHERWISE,
C IB IS SMALLEST INDEX SUCH THAT XB.LT.X(IB+1) .
C
C ......COMPUTE THE INTEGRAL.
IF (IB .LT. IA) THEN
C THIS MEANS IB = IA-1 AND
C (A,B) IS A SUBSET OF (X(IB),X(IA)).
C ------------------------------------------
VALUE = CHFIE (X(IB),X(IA), F(1,IB),F(1,IA),
+ D(1,IB),D(1,IA), A, B)
C ------------------------------------------
ELSE
C
C FIRST COMPUTE INTEGRAL OVER (X(IA),X(IB)).
C (Case (IB .EQ. IA) is taken care of by initialization
C of VALUE to ZERO.)
IF (IB .GT. IA) THEN
C ---------------------------------------------
VALUE = PCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERD)
C ---------------------------------------------
IF (IERD .LT. 0) GO TO 5004
ENDIF
C
C THEN ADD ON INTEGRAL OVER (XA,X(IA)).
IF (XA .LT. X(IA)) THEN
IL = MAX(1, IA-1)
IR = IL + 1
C -------------------------------------
VALUE = VALUE + CHFIE (X(IL),X(IR), F(1,IL),F(1,IR),
+ D(1,IL),D(1,IR), XA, X(IA))
C -------------------------------------
ENDIF
C
C THEN ADD ON INTEGRAL OVER (X(IB),XB).
IF (XB .GT. X(IB)) THEN
IR = MIN (IB+1, N)
IL = IR - 1
C -------------------------------------
VALUE = VALUE + CHFIE (X(IL),X(IR), F(1,IL),F(1,IR),
+ D(1,IL),D(1,IR), X(IB), XB)
C -------------------------------------
ENDIF
C
C FINALLY, ADJUST SIGN IF NECESSARY.
IF (A .GT. B) VALUE = -VALUE
ENDIF
ENDIF
ENDIF
C
C NORMAL RETURN.
C
5000 CONTINUE
PCHIA = VALUE
RETURN
C
C ERROR RETURNS.
C
5001 CONTINUE
C N.LT.2 RETURN.
IERR = -1
CALL XERMSG ('SLATEC', 'PCHIA',
+ 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
GO TO 5000
C
5002 CONTINUE
C INCFD.LT.1 RETURN.
IERR = -2
CALL XERMSG ('SLATEC', 'PCHIA', 'INCREMENT LESS THAN ONE', IERR,
+ 1)
GO TO 5000
C
5003 CONTINUE
C X-ARRAY NOT STRICTLY INCREASING.
IERR = -3
CALL XERMSG ('SLATEC', 'PCHIA',
+ 'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
GO TO 5000
C
5004 CONTINUE
C TROUBLE IN PCHID. (SHOULD NEVER OCCUR.)
IERR = -4
CALL XERMSG ('SLATEC', 'PCHIA', 'TROUBLE IN PCHID', IERR, 1)
GO TO 5000
C------------- LAST LINE OF PCHIA FOLLOWS ------------------------------
END