2
SUBROUTINE DPCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF,
4
C***BEGIN PROLOGUE DPCHBS
5
C***PURPOSE Piecewise Cubic Hermite to B-Spline converter.
6
C***LIBRARY SLATEC (PCHIP)
8
C***TYPE DOUBLE PRECISION (PCHBS-S, DPCHBS-D)
9
C***KEYWORDS B-SPLINES, CONVERSION, CUBIC HERMITE INTERPOLATION,
10
C PIECEWISE CUBIC INTERPOLATION
11
C***AUTHOR Fritsch, F. N., (LLNL)
12
C Computing and Mathematics Research Division
13
C Lawrence Livermore National Laboratory
14
C P.O. Box 808 (L-316)
16
C FTS 532-4275, (510) 422-4275
21
C INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR
22
C PARAMETER (INCFD = ...)
23
C DOUBLE PRECISION X(nmax), F(INCFD,nmax), D(INCFD,nmax),
24
C * T(2*nmax+4), BCOEF(2*nmax)
26
C CALL DPCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF,
31
C N:IN is the number of data points, N.ge.2 . (not checked)
33
C X:IN is the real array of independent variable values. The
34
C elements of X must be strictly increasing:
35
C X(I-1) .LT. X(I), I = 2(1)N. (not checked)
36
C nmax, the dimension of X, must be .ge.N.
38
C F:IN is the real array of dependent variable values.
39
C F(1+(I-1)*INCFD) is the value corresponding to X(I).
40
C nmax, the second dimension of F, must be .ge.N.
42
C D:IN is the real array of derivative values at the data points.
43
C D(1+(I-1)*INCFD) is the value corresponding to X(I).
44
C nmax, the second dimension of D, must be .ge.N.
46
C INCFD:IN is the increment between successive values in F and D.
47
C This argument is provided primarily for 2-D applications.
48
C It may have the value 1 for one-dimensional applications,
49
C in which case F and D may be singly-subscripted arrays.
51
C KNOTYP:IN is a flag to control the knot sequence.
52
C The knot sequence T is normally computed from X by putting
53
C a double knot at each X and setting the end knot pairs
54
C according to the value of KNOTYP:
55
C KNOTYP = 0: Quadruple knots at X(1) and X(N). (default)
56
C KNOTYP = 1: Replicate lengths of extreme subintervals:
57
C T( 1 ) = T( 2 ) = X(1) - (X(2)-X(1)) ;
58
C T(M+4) = T(M+3) = X(N) + (X(N)-X(N-1)).
59
C KNOTYP = 2: Periodic placement of boundary knots:
60
C T( 1 ) = T( 2 ) = X(1) - (X(N)-X(N-1));
61
C T(M+4) = T(M+3) = X(N) + (X(2)-X(1)) .
63
C If the input value of KNOTYP is negative, however, it is
64
C assumed that NKNOTS and T were set in a previous call.
65
C This option is provided for improved efficiency when used
66
C in a parametric setting.
68
C NKNOTS:INOUT is the number of knots.
69
C If KNOTYP.GE.0, then NKNOTS will be set to NDIM+4.
70
C If KNOTYP.LT.0, then NKNOTS is an input variable, and an
71
C error return will be taken if it is not equal to NDIM+4.
73
C T:INOUT is the array of 2*N+4 knots for the B-representation.
74
C If KNOTYP.GE.0, T will be returned by DPCHBS with the
75
C interior double knots equal to the X-values and the
76
C boundary knots set as indicated above.
77
C If KNOTYP.LT.0, it is assumed that T was set by a
78
C previous call to DPCHBS. (This routine does **not**
79
C verify that T forms a legitimate knot sequence.)
81
C BCOEF:OUT is the array of 2*N B-spline coefficients.
83
C NDIM:OUT is the dimension of the B-spline space. (Set to 2*N.)
85
C KORD:OUT is the order of the B-spline. (Set to 4.)
87
C IERR:OUT is an error flag.
89
C IERR = 0 (no errors).
90
C "Recoverable" errors:
91
C IERR = -4 if KNOTYP.GT.2 .
92
C IERR = -5 if KNOTYP.LT.0 and NKNOTS.NE.(2*N+4).
95
C DPCHBS computes the B-spline representation of the PCH function
96
C determined by N,X,F,D. To be compatible with the rest of PCHIP,
97
C DPCHBS includes INCFD, the increment between successive values of
98
C the F- and D-arrays.
100
C The output is the B-representation for the function: NKNOTS, T,
104
C Since it is assumed that the input PCH function has been
105
C computed by one of the other routines in the package PCHIP,
106
C input arguments N, X, INCFD are **not** checked for validity.
108
C *Restrictions/assumptions:
109
C 1. N.GE.2 . (not checked)
110
C 2. X(i).LT.X(i+1), i=1,...,N . (not checked)
111
C 3. INCFD.GT.0 . (not checked)
112
C 4. KNOTYP.LE.2 . (error return if not)
113
C *5. NKNOTS = NDIM+4 = 2*N+4 . (error return if not)
114
C *6. T(2*k+1) = T(2*k) = X(k), k=1,...,N . (not checked)
116
C * Indicates this applies only if KNOTYP.LT.0 .
119
C Argument INCFD is used only to cause the compiler to generate
120
C efficient code for the subscript expressions (1+(I-1)*INCFD) .
121
C The normal usage, in which DPCHBS is called with one-dimensional
122
C arrays F and D, is probably non-Fortran 77, in the strict sense,
123
C but it works on all systems on which DPCHBS has been tested.
126
C PCHIC, PCHIM, or PCHSP can be used to determine an interpolating
127
C PCH function from a set of data.
128
C The B-spline routine DBVALU can be used to evaluate the
129
C B-representation that is output by DPCHBS.
130
C (See BSPDOC for more information.)
132
C***REFERENCES F. N. Fritsch, "Representations for parametric cubic
133
C splines," Computer Aided Geometric Design 6 (1989),
135
C***ROUTINES CALLED DPCHKT, XERMSG
136
C***REVISION HISTORY (YYMMDD)
137
C 870701 DATE WRITTEN
138
C 900405 Converted Fortran to upper case.
139
C 900405 Removed requirement that X be dimensioned N+1.
140
C 900406 Modified to make PCHKT a subsidiary routine to simplify
141
C usage. In the process, added argument INCFD to be com-
142
C patible with the rest of PCHIP.
143
C 900410 Converted prologue to SLATEC 4.0 format.
144
C 900410 Added calls to XERMSG and changed constant 3. to 3 to
145
C reduce single/double differences.
146
C 900411 Added reference.
147
C 900430 Produced double precision version.
148
C 900501 Corrected declarations.
149
C 930317 Minor cosmetic changes. (FNF)
150
C 930514 Corrected problems with dimensioning of arguments and
151
C clarified DESCRIPTION. (FNF)
152
C 930604 Removed NKNOTS from DPCHKT call list. (FNF)
153
C***END PROLOGUE DPCHBS
161
INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR
162
DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*), T(*), BCOEF(*)
164
C Declare local variables.
167
DOUBLE PRECISION DOV3, HNEW, HOLD
168
CHARACTER*8 LIBNAM, SUBNAM
169
C***FIRST EXECUTABLE STATEMENT DPCHBS
179
C Check argument validity. Set up knot sequence if OK.
181
IF ( KNOTYP.GT.2 ) THEN
183
CALL XERMSG (LIBNAM, SUBNAM, 'KNOTYP GREATER THAN 2', IERR, 1)
186
IF ( KNOTYP.LT.0 ) THEN
187
IF ( NKNOTS.NE.NDIM+4 ) THEN
189
CALL XERMSG (LIBNAM, SUBNAM,
190
* 'KNOTYP.LT.0 AND NKNOTS.NE.(2*N+4)', IERR, 1)
194
C Set up knot sequence.
196
CALL DPCHKT (N, X, KNOTYP, T)
199
C Compute B-spline coefficients.
205
C The following requires mixed mode arithmetic.
207
BCOEF(KK-1) = F(1,K) - HOLD*DOV3
208
C The following assumes T(2*K+1) = X(K).
209
HNEW = T(KK+3) - T(KK+1)
210
BCOEF(KK) = F(1,K) + HNEW*DOV3
216
C------------- LAST LINE OF DPCHBS FOLLOWS -----------------------------
added
removed
deps/openlibm/slatec/dpchbs.f
