*DECK PCHCM
SUBROUTINE PCHCM (N, X, F, D, INCFD, SKIP, ISMON, IERR)
C***BEGIN PROLOGUE PCHCM
C***PURPOSE Check a cubic Hermite function for monotonicity.
C***LIBRARY SLATEC (PCHIP)
C***CATEGORY E3
C***TYPE SINGLE PRECISION (PCHCM-S, DPCHCM-D)
C***KEYWORDS CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
C PCHIP, PIECEWISE CUBIC INTERPOLATION, UTILITY ROUTINE
C***AUTHOR Fritsch, F. N., (LLNL)
C Computing & Mathematics Research Division
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 *Usage:
C
C PARAMETER (INCFD = ...)
C INTEGER N, ISMON(N), IERR
C REAL X(N), F(INCFD,N), D(INCFD,N)
C LOGICAL SKIP
C
C CALL PCHCM (N, X, F, D, INCFD, SKIP, ISMON, IERR)
C
C *Arguments:
C
C N:IN is the number of data points. (Error return if N.LT.2 .)
C
C X:IN is a 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:IN is a real array of function values. F(1+(I-1)*INCFD) is
C the value corresponding to X(I).
C
C D:IN is a real array of derivative values. D(1+(I-1)*INCFD) is
C the value corresponding to X(I).
C
C INCFD:IN is the increment between successive values in F and D.
C (Error return if INCFD.LT.1 .)
C
C SKIP:INOUT is a 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.
C SKIP will be set to .TRUE. on normal return.
C
C ISMON:OUT is an integer array indicating on which intervals the
C PCH function defined by N, X, F, D is monotonic.
C For data interval [X(I),X(I+1)],
C ISMON(I) = -3 if function is probably decreasing;
C ISMON(I) = -1 if function is strictly decreasing;
C ISMON(I) = 0 if function is constant;
C ISMON(I) = 1 if function is strictly increasing;
C ISMON(I) = 2 if function is non-monotonic;
C ISMON(I) = 3 if function is probably increasing.
C If ABS(ISMON)=3, this means that the D-values are near
C the boundary of the monotonicity region. A small
C increase produces non-monotonicity; decrease, strict
C monotonicity.
C The above applies to I=1(1)N-1. ISMON(N) indicates whether
C the entire function is monotonic on [X(1),X(N)].
C
C IERR:OUT is an error flag.
C Normal return:
C IERR = 0 (no errors).
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 (The ISMON-array has not been changed 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
C *Description:
C
C PCHCM: Piecewise Cubic Hermite -- Check Monotonicity.
C
C Checks the piecewise cubic Hermite function defined by N,X,F,D
C for monotonicity.
C
C To provide compatibility with PCHIM and PCHIC, includes an
C increment between successive values of the F- and D-arrays.
C
C *Cautions:
C This provides the same capability as old PCHMC, except that a
C new output value, -3, was added February 1989. (Formerly, -3
C and +3 were lumped together in the single value 3.) Codes that
C flag nonmonotonicity by "IF (ISMON.EQ.2)" need not be changed.
C Codes that check via "IF (ISMON.GE.3)" should change the test to
C "IF (IABS(ISMON).GE.3)". Codes that declare monotonicity via
C "IF (ISMON.LE.1)" should change to "IF (IABS(ISMON).LE.1)".
C
C***REFERENCES F. N. Fritsch and R. E. Carlson, Monotone piecewise
C cubic interpolation, SIAM Journal on Numerical Ana-
C lysis 17, 2 (April 1980), pp. 238-246.
C***ROUTINES CALLED CHFCM, XERMSG
C***REVISION HISTORY (YYMMDD)
C 820518 DATE WRITTEN
C 820804 Converted to SLATEC library version.
C 831201 Reversed order of subscripts of F and D, so that the
C routine will work properly when INCFD.GT.1 . (Bug!!)
C 870707 Minor cosmetic changes to prologue.
C 890208 Added possible ISMON value of -3 and modified code so
C that 1,3,-1 produces ISMON(N)=2, rather than 3.
C 890306 Added caution about changed output.
C 890407 Changed name from PCHMC to PCHCM, as requested at the
C March 1989 SLATEC CML meeting, and made a few other
C minor modifications necessitated by this change.
C 890407 Converted to new SLATEC format.
C 890407 Modified DESCRIPTION to LDOC format.
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 920429 Revised format and order of references. (WRB,FNF)
C***END PROLOGUE PCHCM
C
C Fortran intrinsics used: ISIGN.
C Other routines used: CHFCM, XERMSG.
C
C ----------------------------------------------------------------------
C
C Programming notes:
C
C An alternate organization would have separate loops for computing
C ISMON(i), i=1,...,NSEG, and for the computation of ISMON(N). The
C first loop can be readily parallelized, since the NSEG calls to
C CHFCM are independent. The second loop can be cut short if
C ISMON(N) is ever equal to 2, for it cannot be changed further.
C
C To produce a double precision version, simply:
C a. Change PCHCM to DPCHCM wherever it occurs,
C b. Change CHFCM to DCHFCM wherever it occurs, and
C c. Change the real declarations to double precision.
C
C DECLARE ARGUMENTS.
C
INTEGER N, INCFD, ISMON(N), IERR
REAL X(N), F(INCFD,N), D(INCFD,N)
LOGICAL SKIP
C
C DECLARE LOCAL VARIABLES.
C
INTEGER I, NSEG
REAL DELTA
INTEGER CHFCM
C
C VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT PCHCM
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
SKIP = .TRUE.
C
C FUNCTION DEFINITION IS OK -- GO ON.
C
5 CONTINUE
NSEG = N - 1
DO 90 I = 1, NSEG
DELTA = (F(1,I+1)-F(1,I))/(X(I+1)-X(I))
C -------------------------------
ISMON(I) = CHFCM (D(1,I), D(1,I+1), DELTA)
C -------------------------------
IF (I .EQ. 1) THEN
ISMON(N) = ISMON(1)
ELSE
C Need to figure out cumulative monotonicity from following
C "multiplication table":
C
C + I S M O N (I)
C + -3 -1 0 1 3 2
C +------------------------+
C I -3 I -3 -3 -3 2 2 2 I
C S -1 I -3 -1 -1 2 2 2 I
C M 0 I -3 -1 0 1 3 2 I
C O 1 I 2 2 1 1 3 2 I
C N 3 I 2 2 3 3 3 2 I
C (N) 2 I 2 2 2 2 2 2 I
C +------------------------+
C Note that the 2 row and column are out of order so as not
C to obscure the symmetry in the rest of the table.
C
C No change needed if equal or constant on this interval or
C already declared nonmonotonic.
IF ( (ISMON(I).NE.ISMON(N)) .AND. (ISMON(I).NE.0)
. .AND. (ISMON(N).NE.2) ) THEN
IF ( (ISMON(I).EQ.2) .OR. (ISMON(N).EQ.0) ) THEN
ISMON(N) = ISMON(I)
ELSE IF (ISMON(I)*ISMON(N) .LT. 0) THEN
C This interval has opposite sense from curve so far.
ISMON(N) = 2
ELSE
C At this point, both are nonzero with same sign, and
C we have already eliminated case both +-1.
ISMON(N) = ISIGN (3, ISMON(N))
ENDIF
ENDIF
ENDIF
90 CONTINUE
C
C NORMAL RETURN.
C
IERR = 0
RETURN
C
C ERROR RETURNS.
C
5001 CONTINUE
C N.LT.2 RETURN.
IERR = -1
CALL XERMSG ('SLATEC', 'PCHCM',
+ 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
RETURN
C
5002 CONTINUE
C INCFD.LT.1 RETURN.
IERR = -2
CALL XERMSG ('SLATEC', 'PCHCM', 'INCREMENT LESS THAN ONE', IERR,
+ 1)
RETURN
C
5003 CONTINUE
C X-ARRAY NOT STRICTLY INCREASING.
IERR = -3
CALL XERMSG ('SLATEC', 'PCHCM', 'X-ARRAY NOT STRICTLY INCREASING'
+ , IERR, 1)
RETURN
C------------- LAST LINE OF PCHCM FOLLOWS ------------------------------
END