GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dpchim.f
Go to the documentation of this file.
1 *DECK DPCHIM
2  SUBROUTINE dpchim (N, X, F, D, INCFD, IERR)
3 C***BEGIN PROLOGUE DPCHIM
4 C***PURPOSE Set derivatives needed to determine a monotone piecewise
5 C cubic Hermite interpolant to given data. Boundary values
6 C are provided which are compatible with monotonicity. The
7 C interpolant will have an extremum at each point where mono-
8 C tonicity switches direction. (See DPCHIC if user control
9 C is desired over boundary or switch conditions.)
10 C***LIBRARY SLATEC (PCHIP)
11 C***CATEGORY E1A
12 C***TYPE DOUBLE PRECISION (PCHIM-S, DPCHIM-D)
13 C***KEYWORDS CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
14 C PCHIP, PIECEWISE CUBIC INTERPOLATION
15 C***AUTHOR Fritsch, F. N., (LLNL)
16 C Lawrence Livermore National Laboratory
17 C P.O. Box 808 (L-316)
18 C Livermore, CA 94550
19 C FTS 532-4275, (510) 422-4275
20 C***DESCRIPTION
21 C
22 C DPCHIM: Piecewise Cubic Hermite Interpolation to
23 C Monotone data.
24 C
25 C Sets derivatives needed to determine a monotone piecewise cubic
26 C Hermite interpolant to the data given in X and F.
27 C
28 C Default boundary conditions are provided which are compatible
29 C with monotonicity. (See DPCHIC if user control of boundary con-
30 C ditions is desired.)
31 C
32 C If the data are only piecewise monotonic, the interpolant will
33 C have an extremum at each point where monotonicity switches direc-
34 C tion. (See DPCHIC if user control is desired in such cases.)
35 C
36 C To facilitate two-dimensional applications, includes an increment
37 C between successive values of the F- and D-arrays.
38 C
39 C The resulting piecewise cubic Hermite function may be evaluated
40 C by DPCHFE or DPCHFD.
41 C
42 C ----------------------------------------------------------------------
43 C
44 C Calling sequence:
45 C
46 C PARAMETER (INCFD = ...)
47 C INTEGER N, IERR
48 C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N)
49 C
50 C CALL DPCHIM (N, X, F, D, INCFD, IERR)
51 C
52 C Parameters:
53 C
54 C N -- (input) number of data points. (Error return if N.LT.2 .)
55 C If N=2, simply does linear interpolation.
56 C
57 C X -- (input) real*8 array of independent variable values. The
58 C elements of X must be strictly increasing:
59 C X(I-1) .LT. X(I), I = 2(1)N.
60 C (Error return if not.)
61 C
62 C F -- (input) real*8 array of dependent variable values to be
63 C interpolated. F(1+(I-1)*INCFD) is value corresponding to
64 C X(I). DPCHIM is designed for monotonic data, but it will
65 C work for any F-array. It will force extrema at points where
66 C monotonicity switches direction. If some other treatment of
67 C switch points is desired, DPCHIC should be used instead.
68 C -----
69 C D -- (output) real*8 array of derivative values at the data
70 C points. If the data are monotonic, these values will
71 C determine a monotone cubic Hermite function.
72 C The value corresponding to X(I) is stored in
73 C D(1+(I-1)*INCFD), I=1(1)N.
74 C No other entries in D are changed.
75 C
76 C INCFD -- (input) increment between successive values in F and D.
77 C This argument is provided primarily for 2-D applications.
78 C (Error return if INCFD.LT.1 .)
79 C
80 C IERR -- (output) error flag.
81 C Normal return:
82 C IERR = 0 (no errors).
83 C Warning error:
84 C IERR.GT.0 means that IERR switches in the direction
85 C of monotonicity were detected.
86 C "Recoverable" errors:
87 C IERR = -1 if N.LT.2 .
88 C IERR = -2 if INCFD.LT.1 .
89 C IERR = -3 if the X-array is not strictly increasing.
90 C (The D-array has not been changed in any of these cases.)
91 C NOTE: The above errors are checked in the order listed,
92 C and following arguments have **NOT** been validated.
93 C
94 C***REFERENCES 1. F. N. Fritsch and J. Butland, A method for construc-
95 C ting local monotone piecewise cubic interpolants, SIAM
96 C Journal on Scientific and Statistical Computing 5, 2
97 C (June 1984), pp. 300-304.
98 C 2. F. N. Fritsch and R. E. Carlson, Monotone piecewise
99 C cubic interpolation, SIAM Journal on Numerical Ana-
100 C lysis 17, 2 (April 1980), pp. 238-246.
101 C***ROUTINES CALLED DPCHST, XERMSG
102 C***REVISION HISTORY (YYMMDD)
103 C 811103 DATE WRITTEN
104 C 820201 1. Introduced DPCHST to reduce possible over/under-
105 C flow problems.
106 C 2. Rearranged derivative formula for same reason.
107 C 820602 1. Modified end conditions to be continuous functions
108 C of data when monotonicity switches in next interval.
109 C 2. Modified formulas so end conditions are less prone
110 C of over/underflow problems.
111 C 820803 Minor cosmetic changes for release 1.
112 C 870707 Corrected XERROR calls for d.p. name(s).
113 C 870813 Updated Reference 1.
114 C 890206 Corrected XERROR calls.
115 C 890411 Added SAVE statements (Vers. 3.2).
116 C 890531 Changed all specific intrinsics to generic. (WRB)
117 C 890703 Corrected category record. (WRB)
118 C 890831 Modified array declarations. (WRB)
119 C 891006 Cosmetic changes to prologue. (WRB)
120 C 891006 REVISION DATE from Version 3.2
121 C 891214 Prologue converted to Version 4.0 format. (BAB)
122 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
123 C 920429 Revised format and order of references. (WRB,FNF)
124 C***END PROLOGUE DPCHIM
125 C Programming notes:
126 C
127 C 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if
128 C either argument is zero, +1 if they are of the same sign, and
129 C -1 if they are of opposite sign.
130 C 2. To produce a single precision version, simply:
131 C a. Change DPCHIM to PCHIM wherever it occurs,
132 C b. Change DPCHST to PCHST wherever it occurs,
133 C c. Change all references to the Fortran intrinsics to their
134 C single precision equivalents,
135 C d. Change the double precision declarations to real, and
136 C e. Change the constants ZERO and THREE to single precision.
137 C
138 C DECLARE ARGUMENTS.
139 C
140  INTEGER N, INCFD, IERR
141  DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*)
142 C
143 C DECLARE LOCAL VARIABLES.
144 C
145  INTEGER I, NLESS1
146  DOUBLE PRECISION DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE,
147  * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO
148  SAVE zero, three
149  DOUBLE PRECISION DPCHST
150  DATA zero /0.d0/, three/3.d0/
151 C
152 C VALIDITY-CHECK ARGUMENTS.
153 C
154 C***FIRST EXECUTABLE STATEMENT DPCHIM
155  IF ( n.LT.2 ) GO TO 5001
156  IF ( incfd.LT.1 ) GO TO 5002
157  DO 1 i = 2, n
158  IF ( x(i).LE.x(i-1) ) GO TO 5003
159  1 CONTINUE
160 C
161 C FUNCTION DEFINITION IS OK, GO ON.
162 C
163  ierr = 0
164  nless1 = n - 1
165  h1 = x(2) - x(1)
166  del1 = (f(1,2) - f(1,1))/h1
167  dsave = del1
168 C
169 C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
170 C
171  IF (nless1 .GT. 1) GO TO 10
172  d(1,1) = del1
173  d(1,n) = del1
174  GO TO 5000
175 C
176 C NORMAL CASE (N .GE. 3).
177 C
178  10 CONTINUE
179  h2 = x(3) - x(2)
180  del2 = (f(1,3) - f(1,2))/h2
181 C
182 C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
183 C SHAPE-PRESERVING.
184 C
185  hsum = h1 + h2
186  w1 = (h1 + hsum)/hsum
187  w2 = -h1/hsum
188  d(1,1) = w1*del1 + w2*del2
189  IF ( dpchst(d(1,1),del1) .LE. zero) THEN
190  d(1,1) = zero
191  ELSE IF ( dpchst(del1,del2) .LT. zero) THEN
192 C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
193  dmax = three*del1
194  IF (abs(d(1,1)) .GT. abs(dmax)) d(1,1) = dmax
195  ENDIF
196 C
197 C LOOP THROUGH INTERIOR POINTS.
198 C
199  DO 50 i = 2, nless1
200  IF (i .EQ. 2) GO TO 40
201 C
202  h1 = h2
203  h2 = x(i+1) - x(i)
204  hsum = h1 + h2
205  del1 = del2
206  del2 = (f(1,i+1) - f(1,i))/h2
207  40 CONTINUE
208 C
209 C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
210 C
211  d(1,i) = zero
212  IF ( dpchst(del1,del2) .LT. 0.) GO TO 42
213  IF ( dpchst(del1,del2) .EQ. 0.) GO TO 41
214  GO TO 45
215 C
216 C COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY.
217 C
218  41 CONTINUE
219  IF (del2 .EQ. zero) GO TO 50
220  IF ( dpchst(dsave,del2) .LT. zero) ierr = ierr + 1
221  dsave = del2
222  GO TO 50
223 C
224  42 CONTINUE
225  ierr = ierr + 1
226  dsave = del2
227  GO TO 50
228 C
229 C USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
230 C
231  45 CONTINUE
232  hsumt3 = hsum+hsum+hsum
233  w1 = (hsum + h1)/hsumt3
234  w2 = (hsum + h2)/hsumt3
235  dmax = max( abs(del1), abs(del2) )
236  dmin = min( abs(del1), abs(del2) )
237  drat1 = del1/dmax
238  drat2 = del2/dmax
239  d(1,i) = dmin/(w1*drat1 + w2*drat2)
240 C
241  50 CONTINUE
242 C
243 C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
244 C SHAPE-PRESERVING.
245 C
246  w1 = -h2/hsum
247  w2 = (h2 + hsum)/hsum
248  d(1,n) = w1*del1 + w2*del2
249  IF ( dpchst(d(1,n),del2) .LE. zero) THEN
250  d(1,n) = zero
251  ELSE IF ( dpchst(del1,del2) .LT. zero) THEN
252 C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
253  dmax = three*del2
254  IF (abs(d(1,n)) .GT. abs(dmax)) d(1,n) = dmax
255  ENDIF
256 C
257 C NORMAL RETURN.
258 C
259  5000 CONTINUE
260  RETURN
261 C
262 C ERROR RETURNS.
263 C
264  5001 CONTINUE
265 C N.LT.2 RETURN.
266  ierr = -1
267  CALL xermsg ('SLATEC', 'DPCHIM',
268  + 'NUMBER OF DATA POINTS LESS THAN TWO', ierr, 1)
269  RETURN
270 C
271  5002 CONTINUE
272 C INCFD.LT.1 RETURN.
273  ierr = -2
274  CALL xermsg ('SLATEC', 'DPCHIM', 'INCREMENT LESS THAN ONE', ierr,
275  + 1)
276  RETURN
277 C
278  5003 CONTINUE
279 C X-ARRAY NOT STRICTLY INCREASING.
280  ierr = -3
281  CALL xermsg ('SLATEC', 'DPCHIM',
282  + 'X-ARRAY NOT STRICTLY INCREASING', ierr, 1)
283  RETURN
284 C------------- LAST LINE OF DPCHIM FOLLOWS -----------------------------
285  END
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
subroutine dpchim(N, X, F, D, INCFD, IERR)
Definition: dpchim.f:3
static T abs(T x)
Definition: pr-output.cc:1678
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:3