Piecewise Cubic Hermite Interpolation Coefficients.
Sets derivatives needed to determine a piecewise monotone piecewise cubic interpolant to the data given in X and F satisfying the boundary conditions specified by IC and VC.
The treatment of points where monotonicity switches direction is controlled by argument SWITCH.
To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays.
The resulting piecewise cubic Hermite function may be evaluated by dpchfe or dpchfd.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | ic(2) |
array of length 2 specifying desired boundary conditions:
Valid values for IBEG:
Allowable values for the magnitude of IBEG are:
NOTES (IBEG):
IEND may take on the same values as IBEG, but applied to derivative at X(N). In case IEND = 1 or 2, the value is given in VC(2). NOTES (IEND):
|
||
| real(kind=wp), | intent(in) | :: | vc(2) |
array of length 2 specifying desired boundary values, as indicated above.
|
||
| real(kind=wp), | intent(in) | :: | switch |
indicates desired treatment of points where direction of monotonicity switches: Set SWITCH to zero if interpolant is required to be mono- tonic in each interval, regardless of monotonicity of data. NOTES:
Set SWITCH nonzero to use a formula based on the 3-point difference formula in the vicinity of switch points. If SWITCH is positive, the interpolant on each interval containing an extremum is controlled to not deviate from the data by more than SWITCH*DFLOC, where DFLOC is the maximum of the change of F on this interval and its two immediate neighbors. If SWITCH is negative, no such control is to be imposed. |
||
| integer, | intent(in) | :: | n |
number of data points. (Error return if N<2). |
||
| real(kind=wp), | intent(in) | :: | x(*) |
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,*) |
array of dependent variable values to be interpolated. F(1+(I-1)*INCFD) is value corresponding to X(I). |
||
| real(kind=wp), | intent(out) | :: | d(incfd,*) |
array of derivative values at the data points. These values will determine a monotone cubic Hermite function on each subinterval on which the data are monotonic, except possibly adjacent to switches in monotonicity. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed. |
||
| integer, | intent(in) | :: | incfd |
increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD<1). |
||
| real(kind=wp), | intent(inout) | :: | wk(nwk) |
array of working storage. The user may wish to know that the returned values are:
for |
||
| integer, | intent(in) | :: | nwk |
length of work array. (Error return if NWK<2*(N-1)). |
||
| integer, | intent(out) | :: | ierr |
error flag.
(The D-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 dpchic (ic, vc, switch, n, x, f, d, incfd, wk, nwk, ierr) integer,intent(in) :: ic(2) !! array of length 2 specifying desired boundary conditions: !! !! * IC(1) = IBEG, desired condition at beginning of data. !! * IC(2) = IEND, desired condition at end of data. !! !! Valid values for IBEG: !! !! * IBEG = 0 for the default boundary condition (the same as used by [[dpchim]]). !! * If IBEG/=0, then its sign indicates whether the boundary !! derivative is to be adjusted, if necessary, to be !! compatible with monotonicity: !! !! * IBEG>0 if no adjustment is to be performed. !! * IBEG<0 if the derivative is to be adjusted for monotonicity. !! !! Allowable values for the magnitude of IBEG are: !! !! * IBEG = 1 if first derivative at X(1) is given in VC(1). !! * IBEG = 2 if second derivative at X(1) is given in VC(1). !! * IBEG = 3 to use the 3-point difference formula for D(1). !! (Reverts to the default b.c. if N<3). !! * IBEG = 4 to use the 4-point difference formula for D(1). !! (Reverts to the default b.c. if N<4). !! * IBEG = 5 to set D(1) so that the second derivative is continuous !! at X(2). (Reverts to the default b.c. if N<4.) !! This option is somewhat analogous to the "not a knot" !! boundary condition provided by DPCHSP. !! !! NOTES (IBEG): !! !! 1. An error return is taken if ABS(IBEG)>5 . !! 2. Only in case IBEG<=0 is it guaranteed that the !! interpolant will be monotonic in the first interval. !! If the returned value of D(1) lies between zero and !! 3*SLOPE(1), the interpolant will be monotonic. This !! is **NOT** checked if IBEG>0 . !! 3. If IBEG<0 and D(1) had to be changed to achieve mono- !! tonicity, a warning error is returned. !! !! IEND may take on the same values as IBEG, but applied to !! derivative at X(N). In case IEND = 1 or 2, the value is !! given in VC(2). !! !! NOTES (IEND): !! !! 1. An error return is taken if ABS(IEND)>5 . !! 2. Only in case IEND<=0 is it guaranteed that the !! interpolant will be monotonic in the last interval. !! If the returned value of D(1+(N-1)*INCFD) lies between !! zero and 3*SLOPE(N-1), the interpolant will be monotonic. !! This is **NOT** checked if IEND>0 . !! 3. If IEND<0 and D(1+(N-1)*INCFD) had to be changed to !! achieve monotonicity, a warning error is returned. integer,intent(in) :: n !! number of data points. (Error return if N<2). integer,intent(in) :: incfd !! increment between successive values in F and D. !! This argument is provided primarily for 2-D applications. !! (Error return if INCFD<1). integer,intent(in) :: nwk !! length of work array. (Error return if NWK<2*(N-1)). integer,intent(out) :: ierr !! error flag. !! !! * Normal return: !! !! * IERR = 0 (no errors). !! !! * Warning errors: !! !! * IERR = 1 if IBEG<0 and D(1) had to be adjusted for !! monotonicity. !! * IERR = 2 if IEND<0 and D(1+(N-1)*INCFD) had to be !! adjusted for monotonicity. !! * IERR = 3 if both of the above are true. !! !! * "Recoverable" errors: !! !! * IERR = -1 if N<2 . !! * IERR = -2 if INCFD<1 . !! * IERR = -3 if the X-array is not strictly increasing. !! * IERR = -4 if ABS(IBEG)>5 . !! * IERR = -5 if ABS(IEND)>5 . !! * IERR = -6 if both of the above are true. !! * IERR = -7 if NWK<2*(N-1) . !! !! (The D-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. real(wp),intent(in) :: vc(2) !! array of length 2 specifying desired boundary !! values, as indicated above. !! !! * VC(1) need be set only if IC(1) = 1 or 2 . !! * VC(2) need be set only if IC(2) = 1 or 2 . real(wp),intent(in) :: switch !! indicates desired treatment of points where !! direction of monotonicity switches: !! Set SWITCH to zero if interpolant is required to be mono- !! tonic in each interval, regardless of monotonicity of data. !! !! NOTES: !! !! 1. This will cause D to be set to zero at all switch !! points, thus forcing extrema there. !! 2. The result of using this option with the default boun- !! dary conditions will be identical to using DPCHIM, but !! will generally cost more compute time. !! This option is provided only to facilitate comparison !! of different switch and/or boundary conditions. !! !! Set SWITCH nonzero to use a formula based on the 3-point !! difference formula in the vicinity of switch points. !! !! If SWITCH is positive, the interpolant on each interval !! containing an extremum is controlled to not deviate from !! the data by more than SWITCH*DFLOC, where DFLOC is the !! maximum of the change of F on this interval and its two !! immediate neighbors. !! !! If SWITCH is negative, no such control is to be imposed. real(wp),intent(in) :: x(*) !! 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,*) !! array of dependent variable values to be !! interpolated. F(1+(I-1)*INCFD) is value corresponding to !! X(I). real(wp),intent(out) :: d(incfd,*) !! array of derivative values at the data !! points. These values will determine a monotone cubic !! Hermite function on each subinterval on which the data !! are monotonic, except possibly adjacent to switches in !! monotonicity. The value corresponding to X(I) is stored in !! D(1+(I-1)*INCFD), I=1(1)N. !! No other entries in D are changed. real(wp),intent(inout) :: wk(nwk) !! array of working storage. The user may !! wish to know that the returned values are: !! !! * `WK(I) = H(I) = X(I+1) - X(I)` ; !! * `WK(N-1+I) = SLOPE(I) = (F(1,I+1) - F(1,I)) / H(I)` !! !! for `I = 1(1)N-1`. integer :: i, ibeg, iend, nless1 ! validity-check arguments. if ( n<2 ) then ierr = -1 call xermsg ('PCHIP', 'dpchic', 'number of data points less than two', ierr, 1) return end if if ( incfd<1 ) then ierr = -2 call xermsg ('PCHIP', 'dpchic', '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', 'dpchic', 'x-array not strictly increasing', ierr, 1) return end if end do ibeg = ic(1) iend = ic(2) ierr = 0 if (abs(ibeg) > 5) ierr = ierr - 1 if (abs(iend) > 5) ierr = ierr - 2 if (ierr < 0) then ! ic out of range return. ierr = ierr - 3 call xermsg ('PCHIP', 'dpchic', 'ic out of range', ierr, 1) return end if ! function definition is ok -- go on. nless1 = n - 1 if ( nwk < 2*nless1 ) then ierr = -7 call xermsg ('PCHIP', 'dpchic', 'work array too small', ierr, 1) return end if ! set up h and slope arrays. do i = 1, nless1 wk(i) = x(i+1) - x(i) wk(nless1+i) = (f(1,i+1) - f(1,i)) / wk(i) end do ! special case n=2 -- use linear interpolation. if (nless1 > 1) then ! normal case (n >= 3) . ! set interior derivatives and default end conditions. call dpchci (n, wk(1), wk(n), d, incfd) ! set derivatives at points where monotonicity switches direction. if (switch /= zero) then call dpchcs (switch, n, wk(1), wk(n), d, incfd, ierr) if (ierr /= 0) then ! error return from dpchcs. ierr = -8 call xermsg ('PCHIP', 'dpchic', 'error return from dpchcs', ierr, 1) return end if end if else d(1,1) = wk(2) d(1,n) = wk(2) end if ! set end conditions. if ( (ibeg/=0) .or. (iend/=0) ) then call dpchce (ic, vc, n, x, wk(1), wk(n), d, incfd, ierr) if (ierr < 0) then ! error return from dpchce. ! *** this case should never occur *** ierr = -9 call xermsg ('PCHIP', 'dpchic', 'error return from dpchce', ierr, 1) end if end if end subroutine dpchic