1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
|
*DECK DPCHBS
SUBROUTINE DPCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF,
+ NDIM, KORD, IERR)
C***BEGIN PROLOGUE DPCHBS
C***PURPOSE Piecewise Cubic Hermite to B-Spline converter.
C***LIBRARY SLATEC (PCHIP)
C***CATEGORY E3
C***TYPE DOUBLE PRECISION (PCHBS-S, DPCHBS-D)
C***KEYWORDS B-SPLINES, CONVERSION, CUBIC HERMITE INTERPOLATION,
C PIECEWISE CUBIC INTERPOLATION
C***AUTHOR Fritsch, F. N., (LLNL)
C Computing and 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 INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR
C PARAMETER (INCFD = ...)
C DOUBLE PRECISION X(nmax), F(INCFD,nmax), D(INCFD,nmax),
C * T(2*nmax+4), BCOEF(2*nmax)
C
C CALL DPCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF,
C * NDIM, KORD, IERR)
C
C *Arguments:
C
C N:IN is the number of data points, N.ge.2 . (not checked)
C
C X:IN is the 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. (not checked)
C nmax, the dimension of X, must be .ge.N.
C
C F:IN is the real array of dependent variable values.
C F(1+(I-1)*INCFD) is the value corresponding to X(I).
C nmax, the second dimension of F, must be .ge.N.
C
C D:IN is the real array of derivative values at the data points.
C D(1+(I-1)*INCFD) is the value corresponding to X(I).
C nmax, the second dimension of D, must be .ge.N.
C
C INCFD:IN is the increment between successive values in F and D.
C This argument is provided primarily for 2-D applications.
C It may have the value 1 for one-dimensional applications,
C in which case F and D may be singly-subscripted arrays.
C
C KNOTYP:IN is a flag to control the knot sequence.
C The knot sequence T is normally computed from X by putting
C a double knot at each X and setting the end knot pairs
C according to the value of KNOTYP:
C KNOTYP = 0: Quadruple knots at X(1) and X(N). (default)
C KNOTYP = 1: Replicate lengths of extreme subintervals:
C T( 1 ) = T( 2 ) = X(1) - (X(2)-X(1)) ;
C T(M+4) = T(M+3) = X(N) + (X(N)-X(N-1)).
C KNOTYP = 2: Periodic placement of boundary knots:
C T( 1 ) = T( 2 ) = X(1) - (X(N)-X(N-1));
C T(M+4) = T(M+3) = X(N) + (X(2)-X(1)) .
C Here M=NDIM=2*N.
C If the input value of KNOTYP is negative, however, it is
C assumed that NKNOTS and T were set in a previous call.
C This option is provided for improved efficiency when used
C in a parametric setting.
C
C NKNOTS:INOUT is the number of knots.
C If KNOTYP.GE.0, then NKNOTS will be set to NDIM+4.
C If KNOTYP.LT.0, then NKNOTS is an input variable, and an
C error return will be taken if it is not equal to NDIM+4.
C
C T:INOUT is the array of 2*N+4 knots for the B-representation.
C If KNOTYP.GE.0, T will be returned by DPCHBS with the
C interior double knots equal to the X-values and the
C boundary knots set as indicated above.
C If KNOTYP.LT.0, it is assumed that T was set by a
C previous call to DPCHBS. (This routine does **not**
C verify that T forms a legitimate knot sequence.)
C
C BCOEF:OUT is the array of 2*N B-spline coefficients.
C
C NDIM:OUT is the dimension of the B-spline space. (Set to 2*N.)
C
C KORD:OUT is the order of the B-spline. (Set to 4.)
C
C IERR:OUT is an error flag.
C Normal return:
C IERR = 0 (no errors).
C "Recoverable" errors:
C IERR = -4 if KNOTYP.GT.2 .
C IERR = -5 if KNOTYP.LT.0 and NKNOTS.NE.(2*N+4).
C
C *Description:
C DPCHBS computes the B-spline representation of the PCH function
C determined by N,X,F,D. To be compatible with the rest of PCHIP,
C DPCHBS includes INCFD, the increment between successive values of
C the F- and D-arrays.
C
C The output is the B-representation for the function: NKNOTS, T,
C BCOEF, NDIM, KORD.
C
C *Caution:
C Since it is assumed that the input PCH function has been
C computed by one of the other routines in the package PCHIP,
C input arguments N, X, INCFD are **not** checked for validity.
C
C *Restrictions/assumptions:
C 1. N.GE.2 . (not checked)
C 2. X(i).LT.X(i+1), i=1,...,N . (not checked)
C 3. INCFD.GT.0 . (not checked)
C 4. KNOTYP.LE.2 . (error return if not)
C *5. NKNOTS = NDIM+4 = 2*N+4 . (error return if not)
C *6. T(2*k+1) = T(2*k) = X(k), k=1,...,N . (not checked)
C
C * Indicates this applies only if KNOTYP.LT.0 .
C
C *Portability:
C Argument INCFD is used only to cause the compiler to generate
C efficient code for the subscript expressions (1+(I-1)*INCFD) .
C The normal usage, in which DPCHBS is called with one-dimensional
C arrays F and D, is probably non-Fortran 77, in the strict sense,
C but it works on all systems on which DPCHBS has been tested.
C
C *See Also:
C PCHIC, PCHIM, or PCHSP can be used to determine an interpolating
C PCH function from a set of data.
C The B-spline routine DBVALU can be used to evaluate the
C B-representation that is output by DPCHBS.
C (See BSPDOC for more information.)
C
C***REFERENCES F. N. Fritsch, "Representations for parametric cubic
C splines," Computer Aided Geometric Design 6 (1989),
C pp.79-82.
C***ROUTINES CALLED DPCHKT, XERMSG
C***REVISION HISTORY (YYMMDD)
C 870701 DATE WRITTEN
C 900405 Converted Fortran to upper case.
C 900405 Removed requirement that X be dimensioned N+1.
C 900406 Modified to make PCHKT a subsidiary routine to simplify
C usage. In the process, added argument INCFD to be com-
C patible with the rest of PCHIP.
C 900410 Converted prologue to SLATEC 4.0 format.
C 900410 Added calls to XERMSG and changed constant 3. to 3 to
C reduce single/double differences.
C 900411 Added reference.
C 900430 Produced double precision version.
C 900501 Corrected declarations.
C 930317 Minor cosmetic changes. (FNF)
C 930514 Corrected problems with dimensioning of arguments and
C clarified DESCRIPTION. (FNF)
C 930604 Removed NKNOTS from DPCHKT call list. (FNF)
C***END PROLOGUE DPCHBS
C
C*Internal Notes:
C
C**End
C
C Declare arguments.
C
INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR
DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*), T(*), BCOEF(*)
C
C Declare local variables.
C
INTEGER K, KK
DOUBLE PRECISION DOV3, HNEW, HOLD
CHARACTER*8 LIBNAM, SUBNAM
C***FIRST EXECUTABLE STATEMENT DPCHBS
C
C Initialize.
C
NDIM = 2*N
KORD = 4
IERR = 0
LIBNAM = 'SLATEC'
SUBNAM = 'DPCHBS'
C
C Check argument validity. Set up knot sequence if OK.
C
IF ( KNOTYP.GT.2 ) THEN
IERR = -1
CALL XERMSG (LIBNAM, SUBNAM, 'KNOTYP GREATER THAN 2', IERR, 1)
RETURN
ENDIF
IF ( KNOTYP.LT.0 ) THEN
IF ( NKNOTS.NE.NDIM+4 ) THEN
IERR = -2
CALL XERMSG (LIBNAM, SUBNAM,
* 'KNOTYP.LT.0 AND NKNOTS.NE.(2*N+4)', IERR, 1)
RETURN
ENDIF
ELSE
C Set up knot sequence.
NKNOTS = NDIM + 4
CALL DPCHKT (N, X, KNOTYP, T)
ENDIF
C
C Compute B-spline coefficients.
C
HNEW = T(3) - T(1)
DO 40 K = 1, N
KK = 2*K
HOLD = HNEW
C The following requires mixed mode arithmetic.
DOV3 = D(1,K)/3
BCOEF(KK-1) = F(1,K) - HOLD*DOV3
C The following assumes T(2*K+1) = X(K).
HNEW = T(KK+3) - T(KK+1)
BCOEF(KK) = F(1,K) + HNEW*DOV3
40 CONTINUE
C
C Terminate.
C
RETURN
C------------- LAST LINE OF DPCHBS FOLLOWS -----------------------------
END
|
