Check a cubic Hermite function for monotonicity
Checks the piecewise cubic Hermite function defined by N,X,F,D for monotonicity.
To provide compatibility with dpchim and dpchic, includes an increment between successive values of the F- and D-arrays.
This provides the same capability as old DPCHMC, except that a
new output value, -3, was added February 1989. (Formerly, -3
and +3 were lumped together in the single value 3.) Codes that
flag nonmonotonicity by "IF (ISMON==2)" need not be changed.
Codes that check via "IF (ISMON>=3)" should change the test to
"IF (IABS(ISMON)>=3)". Codes that declare monotonicity via
"IF (ISMON<=1)" should change to "IF (IABS(ISMON)<=1)".
An alternate organization would have separate loops for computing ISMON(i), i=1,...,NSEG, and for the computation of ISMON(N). The first loop can be readily parallelized, since the NSEG calls to CHFCM are independent. The second loop can be cut short if ISMON(N) is ever equal to 2, for it cannot be changed further.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n |
the number of data points. (Error return if N<2). |
||
| real(kind=wp), | intent(in) | :: | x(n) |
array of independent variable values. The elements of X must be strictly increasing: X(I-1) < X(I), I = 2(1)N. (Error return if not.) |
||
| real(kind=wp), | intent(in) | :: | f(incfd,n) |
array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). |
||
| real(kind=wp), | intent(in) | :: | d(incfd,n) |
array of derivative values. D(1+(I-1)*INCFD) is is the value corresponding to X(I). |
||
| integer, | intent(in) | :: | incfd |
the increment between successive values in F and D. (Error return if INCFD<1). |
||
| logical, | intent(inout) | :: | skip |
logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed. SKIP will be set to .TRUE. on normal return. |
||
| integer, | intent(out) | :: | ismon(n) |
array indicating on which intervals the PCH function defined by N, X, F, D is monotonic. For data interval [X(I),X(I+1)]:
If ABS(ISMON)=3, this means that the D-values are near the boundary of the monotonicity region. A small increase produces non-monotonicity; decrease, strict monotonicity. The above applies to I=1(1)N-1. ISMON(N) indicates whether the entire function is monotonic on [X(1),X(N)]. |
||
| integer, | intent(out) | :: | ierr |
error flag.
(The ISMON-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have NOT been validated. |
subroutine dpchcm (n, x, f, d, incfd, skip, ismon, ierr) integer,intent(in) :: n !! the number of data points. (Error return if N<2). real(wp),intent(in) :: x(n) !! array of independent variable values. The !! elements of X must be strictly increasing: !! X(I-1) < X(I), I = 2(1)N. !! (Error return if not.) real(wp),intent(in) :: f(incfd,n) !! array of function values. F(1+(I-1)*INCFD) is !! the value corresponding to X(I). real(wp),intent(in) :: d(incfd,n) !! array of derivative values. D(1+(I-1)*INCFD) is !! is the value corresponding to X(I). integer,intent(in) :: incfd !! the increment between successive values in F and D. !! (Error return if INCFD<1). logical,intent(inout) :: skip !! logical variable which should be set to !! .TRUE. if the user wishes to skip checks for validity of !! preceding parameters, or to .FALSE. otherwise. !! This will save time in case these checks have already !! been performed. !! SKIP will be set to .TRUE. on normal return. integer,intent(out) :: ismon(n) !! array indicating on which intervals the !! PCH function defined by N, X, F, D is monotonic. !! For data interval [X(I),X(I+1)]: !! !! * ISMON(I) = -3 if function is probably decreasing; !! * ISMON(I) = -1 if function is strictly decreasing; !! * ISMON(I) = 0 if function is constant; !! * ISMON(I) = 1 if function is strictly increasing; !! * ISMON(I) = 2 if function is non-monotonic; !! * ISMON(I) = 3 if function is probably increasing. !! !! If ABS(ISMON)=3, this means that the D-values are near !! the boundary of the monotonicity region. A small !! increase produces non-monotonicity; decrease, strict !! monotonicity. !! !! The above applies to I=1(1)N-1. ISMON(N) indicates whether !! the entire function is monotonic on [X(1),X(N)]. integer,intent(out) :: ierr !! error flag. !! !! * Normal return: !! * IERR = 0 (no errors). !! * "Recoverable" errors: !! * IERR = -1 if N<2 . !! * IERR = -2 if INCFD<1 . !! * IERR = -3 if the X-array is not strictly increasing. !! !! (The ISMON-array has not been changed in any of these cases.) !! NOTE: The above errors are checked in the order listed, !! and following arguments have **NOT** been validated. integer :: i, nseg real(wp) :: delta ! validity-check arguments. if (.not. skip) then if ( n<2 ) then ierr = -1 call xermsg ('PCHIP', 'dpchcm', 'number of data points less than two', ierr, 1) return end if if ( incfd<1 ) then ierr = -2 call xermsg ('PCHIP', 'dpchcm', 'increment less than one', ierr, 1) return end if do i = 2, n if ( x(i)<=x(i-1) ) then ! x-array not strictly increasing. ierr = -3 call xermsg ('PCHIP', 'dpchcm', 'x-array not strictly increasing', ierr, 1) return end if end do skip = .true. end if ! function definition is ok -- go on. nseg = n - 1 do i = 1, nseg delta = (f(1,i+1)-f(1,i))/(x(i+1)-x(i)) ! ------------------------------- ismon(i) = dchfcm (d(1,i), d(1,i+1), delta) ! ------------------------------- if (i == 1) then ismon(n) = ismon(1) else ! Need to figure out cumulative monotonicity from following ! "multiplication table": ! ! + I S M O N (I) ! + -3 -1 0 1 3 2 ! +------------------------+ ! I -3 I -3 -3 -3 2 2 2 I ! S -1 I -3 -1 -1 2 2 2 I ! M 0 I -3 -1 0 1 3 2 I ! O 1 I 2 2 1 1 3 2 I ! N 3 I 2 2 3 3 3 2 I ! (N) 2 I 2 2 2 2 2 2 I ! +------------------------+ ! Note that the 2 row and column are out of order so as not ! to obscure the symmetry in the rest of the table. ! ! No change needed if equal or constant on this interval or ! already declared nonmonotonic. if ( (ismon(i)/=ismon(n)) .and. (ismon(i)/=0) .and. (ismon(n)/=2) ) then if ( (ismon(i)==2) .or. (ismon(n)==0) ) then ismon(n) = ismon(i) else if (ismon(i)*ismon(n) < 0) then ! This interval has opposite sense from curve so far. ismon(n) = 2 else ! At this point, both are nonzero with same sign, and ! we have already eliminated case both +-1. ismon(n) = sign (3, ismon(n)) endif endif endif end do ierr = 0 end subroutine dpchcm