Piecewise Cubic Hermite Spline
Computes the Hermite representation of the cubic spline interpolant to the data given in X and F satisfying the boundary conditions specified by IC and VC.
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.
Note
This is a modified version of C. de Boor's cubic spline routine CUBSPL.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | ic(2) |
integer array of length 2 specifying desired boundary conditions:
Valid values for IBEG:
NOTES:
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:
|
||
| real(kind=wp), | intent(in) | :: | vc(2) |
array of length 2 specifying desired boundary values, as indicated above.
|
||
| 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:
|
||
| 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 the cubic spline
interpolant with the requested boundary conditions.
The value corresponding to X(I) is stored in
|
||
| 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(2,*) |
array of working storage. |
||
| integer, | intent(in) | :: | nwk |
length of work array. (Error return if NWK<2*N). |
||
| integer, | intent(out) | :: | ierr |
error flag.
|
subroutine dpchsp (ic, vc, n, x, f, d, incfd, wk, nwk, ierr) integer,intent(in) :: ic(2) !! integer 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 to set D(1) so that the third derivative is !! continuous at X(2). This is the "not a knot" condition !! provided by de Boor's cubic spline routine CUBSPL. !! **This is the default boundary condition.** !! * 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). !! !! NOTES: !! !! 1. An error return is taken if IBEG is out of range. !! 2. For the "natural" boundary condition, use IBEG=2 and !! VC(1)=0. !! !! 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: !! !! 1. An error return is taken if IEND is out of range. !! 2. For the "natural" boundary condition, use IEND=2 and !! VC(2)=0. 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). 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. !! * IERR = -4 if IBEG<0 or IBEG>4 . !! * IERR = -5 if IEND<0 of IEND>4 . !! * IERR = -6 if both of the above are true. !! * IERR = -7 if NWK is too small. !! !! NOTE: The above errors are checked in the order listed, !! and following arguments have **NOT** been validated. !! (The D-array has not been changed in any of these cases.) !! !! * IERR = -8 in case of trouble solving the linear system !! for the interior derivative values. !! (The D-array may have been changed in this case.) !! ( Do **NOT** use it! ) 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) :: 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 the cubic spline !! interpolant with the requested boundary conditions. !! 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(2,*) !! array of working storage. integer :: ibeg, iend, index, j, nm1 real(wp) :: g, stemp(3), xtemp(4) logical :: forward ! validity-check arguments. if ( n<2 ) then ierr = -1 call xermsg ('PCHIP', 'dpchsp', 'number of data points less than two', ierr, 1) return end if if ( incfd<1 ) then ierr = -2 call xermsg ('PCHIP', 'dpchsp', 'increment less than one', ierr, 1) return end if do j = 2, n if ( x(j)<=x(j-1) ) then ! x-array not strictly increasing. ierr = -3 call xermsg ('PCHIP', 'dpchsp', 'x-array not strictly increasing', ierr, 1) return end if end do ibeg = ic(1) iend = ic(2) ierr = 0 if ( (ibeg<0).or.(ibeg>4) ) ierr = ierr - 1 if ( (iend<0).or.(iend>4) ) ierr = ierr - 2 if ( ierr<0 ) then ! ic out of range return. ierr = ierr - 3 call xermsg ('PCHIP', 'dpchsp', 'ic out of range', ierr, 1) return end if ! function definition is ok -- go on. if ( nwk < 2*n ) then ! nwk too small return. ierr = -7 call xermsg ('PCHIP', 'dpchsp', 'work array too small', ierr, 1) return end if ! compute first differences of x sequence and store in wk(1,.). also, ! compute first divided difference of data and store in wk(2,.). do j=2,n wk(1,j) = x(j) - x(j-1) wk(2,j) = (f(1,j) - f(1,j-1))/wk(1,j) end do ! set to default boundary conditions if n is too small. if ( ibeg>n ) ibeg = 0 if ( iend>n ) iend = 0 ! set up for boundary conditions. if ( (ibeg==1).or.(ibeg==2) ) then d(1,1) = vc(1) else if (ibeg > 2) then ! pick up first ibeg points, in reverse order. do j = 1, ibeg index = ibeg-j+1 ! index runs from ibeg down to 1. xtemp(j) = x(index) if (j < ibeg) stemp(j) = wk(2,index) end do d(1,1) = dpchdf (ibeg, xtemp, stemp, ierr) if (ierr /= 0) then ! error return from dpchdf. ! *** this case should never occur *** ierr = -9 call xermsg ('PCHIP', 'dpchsp', 'error return from dpchdf', ierr, 1) end if ibeg = 1 endif if ( (iend==1).or.(iend==2) ) then d(1,n) = vc(2) else if (iend > 2) then ! pick up last iend points. do j = 1, iend index = n-iend+j ! index runs from n+1-iend up to n. xtemp(j) = x(index) if (j < iend) stemp(j) = wk(2,index+1) end do d(1,n) = dpchdf (iend, xtemp, stemp, ierr) if (ierr /= 0) then ! error return from dpchdf. ! *** this case should never occur *** ierr = -9 call xermsg ('PCHIP', 'dpchsp', 'error return from dpchdf', ierr, 1) end if iend = 1 endif ! --------------------( begin coding from cubspl )-------------------- ! ! **** a tridiagonal linear system for the unknown slopes s(j) of ! f at x(j), j=1,...,n, is generated and then solved by gauss elim- ! ination, with s(j) ending up in d(1,j), all j. ! wk(1,.) and wk(2,.) are used for temporary storage. ! ! construct first equation from first boundary condition, of the form ! wk(2,1)*s(1) + wk(1,1)*s(2) = d(1,1) if (ibeg == 0) then if (n == 2) then ! no condition at left end and n = 2. wk(2,1) = one wk(1,1) = one d(1,1) = two*wk(2,2) else ! not-a-knot condition at left end and n > 2. wk(2,1) = wk(1,3) wk(1,1) = wk(1,2) + wk(1,3) d(1,1) =((wk(1,2) + two*wk(1,1))*wk(2,2)*wk(1,3) + wk(1,2)**2*wk(2,3)) / wk(1,1) endif else if (ibeg == 1) then ! slope prescribed at left end. wk(2,1) = one wk(1,1) = zero else ! second derivative prescribed at left end. wk(2,1) = two wk(1,1) = one d(1,1) = three*wk(2,2) - half*wk(1,2)*d(1,1) endif ! if there are interior knots, generate the corresponding equations and ! carry out the forward pass of gauss elimination, after which the j-th ! equation reads `wk(2,j)*s(j) + wk(1,j)*s(j+1) = d(1,j)`. nm1 = n-1 if (nm1 > 1) then do j=2,nm1 if (wk(2,j-1) == zero) then ! singular system. ! *** theoretically, this can only occur if successive x-values *** ! *** are equal, which should already have been caught (ierr=-3). *** ierr = -8 call xermsg ('PCHIP', 'dpchsp', 'singular linear system', ierr, 1) return end if g = -wk(1,j+1)/wk(2,j-1) d(1,j) = g*d(1,j-1) + three*(wk(1,j)*wk(2,j+1) + wk(1,j+1)*wk(2,j)) wk(2,j) = g*wk(1,j-1) + two*(wk(1,j) + wk(1,j+1)) end do endif ! construct last equation from second boundary condition, of the form ! (-g*wk(2,n-1))*s(n-1) + wk(2,n)*s(n) = d(1,n) ! ! if slope is prescribed at right end, one can go directly to back- ! substitution, since arrays happen to be set up just right for it ! at this point. if (iend /= 1) then forward = .true. if (iend == 0) then if (n==2 .and. ibeg==0) then ! not-a-knot at right endpoint and at left endpoint and n = 2. d(1,2) = wk(2,2) forward = .false. else if ((n==2) .or. (n==3 .and. ibeg==0)) then ! either (n=3 and not-a-knot also at left) or (n=2 and *not* ! not-a-knot at left end point). d(1,n) = two*wk(2,n) wk(2,n) = one if (wk(2,n-1) == zero) then ! singular system. ! *** theoretically, this can only occur if successive x-values *** ! *** are equal, which should already have been caught (ierr=-3). *** ierr = -8 call xermsg ('PCHIP', 'dpchsp', 'singular linear system', ierr, 1) return end if g = -one/wk(2,n-1) else ! not-a-knot and n >= 3, and either n>3 or also not-a- ! knot at left end point. g = wk(1,n-1) + wk(1,n) ! do not need to check following denominators (x-differences). d(1,n) = ((wk(1,n)+two*g)*wk(2,n)*wk(1,n-1) & + wk(1,n)**2*(f(1,n-1)-f(1,n-2))/wk(1,n-1))/g if (wk(2,n-1) == zero) then ! singular system. ! *** theoretically, this can only occur if successive x-values *** ! *** are equal, which should already have been caught (ierr=-3). *** ierr = -8 call xermsg ('PCHIP', 'dpchsp', 'singular linear system', ierr, 1) return end if g = -g/wk(2,n-1) wk(2,n) = wk(1,n-1) endif else ! second derivative prescribed at right endpoint. d(1,n) = three*wk(2,n) + half*wk(1,n)*d(1,n) wk(2,n) = two if (wk(2,n-1) == zero) then ! singular system. ! *** theoretically, this can only occur if successive x-values *** ! *** are equal, which should already have been caught (ierr=-3). *** ierr = -8 call xermsg ('PCHIP', 'dpchsp', 'singular linear system', ierr, 1) return end if g = -one/wk(2,n-1) endif ! complete forward pass of gauss elimination. if (forward) then wk(2,n) = g*wk(1,n-1) + wk(2,n) if (wk(2,n) == zero) then ! singular system. ! *** theoretically, this can only occur if successive x-values *** ! *** are equal, which should already have been caught (ierr=-3). *** ierr = -8 call xermsg ('PCHIP', 'dpchsp', 'singular linear system', ierr, 1) return end if d(1,n) = (g*d(1,n-1) + d(1,n))/wk(2,n) end if end if ! carry out back substitution do j=nm1,1,-1 if (wk(2,j) == zero) then ! singular system. ! *** theoretically, this can only occur if successive x-values *** ! *** are equal, which should already have been caught (ierr=-3). *** ierr = -8 call xermsg ('PCHIP', 'dpchsp', 'singular linear system', ierr, 1) return end if d(1,j) = (d(1,j) - wk(1,j)*d(1,j+1))/wk(2,j) end do end subroutine dpchsp