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
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
|
*DECK PCHCE
SUBROUTINE PCHCE (IC, VC, N, X, H, SLOPE, D, INCFD, IERR)
C***BEGIN PROLOGUE PCHCE
C***SUBSIDIARY
C***PURPOSE Set boundary conditions for PCHIC
C***LIBRARY SLATEC (PCHIP)
C***TYPE SINGLE PRECISION (PCHCE-S, DPCHCE-D)
C***AUTHOR Fritsch, F. N., (LLNL)
C***DESCRIPTION
C
C PCHCE: PCHIC End Derivative Setter.
C
C Called by PCHIC to set end derivatives as requested by the user.
C It must be called after interior derivative values have been set.
C -----
C
C To facilitate two-dimensional applications, includes an increment
C between successive values of the D-array.
C
C ----------------------------------------------------------------------
C
C Calling sequence:
C
C PARAMETER (INCFD = ...)
C INTEGER IC(2), N, IERR
C REAL VC(2), X(N), H(N), SLOPE(N), D(INCFD,N)
C
C CALL PCHCE (IC, VC, N, X, H, SLOPE, D, INCFD, IERR)
C
C Parameters:
C
C IC -- (input) integer array of length 2 specifying desired
C boundary conditions:
C IC(1) = IBEG, desired condition at beginning of data.
C IC(2) = IEND, desired condition at end of data.
C ( see prologue to PCHIC for details. )
C
C VC -- (input) real array of length 2 specifying desired boundary
C values. VC(1) need be set only if IC(1) = 2 or 3 .
C VC(2) need be set only if IC(2) = 2 or 3 .
C
C N -- (input) number of data points. (assumes N.GE.2)
C
C X -- (input) real array of independent variable values. (the
C elements of X are assumed to be strictly increasing.)
C
C H -- (input) real array of interval lengths.
C SLOPE -- (input) real array of data slopes.
C If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
C H(I) = X(I+1)-X(I),
C SLOPE(I) = (Y(I+1)-Y(I))/H(I), I=1(1)N-1.
C
C D -- (input) real array of derivative values at the data points.
C The value corresponding to X(I) must be stored in
C D(1+(I-1)*INCFD), I=1(1)N.
C (output) the value of D at X(1) and/or X(N) is changed, if
C necessary, to produce the requested boundary conditions.
C no other entries in D are changed.
C
C INCFD -- (input) increment between successive values in D.
C This argument is provided primarily for 2-D applications.
C
C IERR -- (output) error flag.
C Normal return:
C IERR = 0 (no errors).
C Warning errors:
C IERR = 1 if IBEG.LT.0 and D(1) had to be adjusted for
C monotonicity.
C IERR = 2 if IEND.LT.0 and D(1+(N-1)*INCFD) had to be
C adjusted for monotonicity.
C IERR = 3 if both of the above are true.
C
C -------
C WARNING: This routine does no validity-checking of arguments.
C -------
C
C Fortran intrinsics used: ABS.
C
C***SEE ALSO PCHIC
C***ROUTINES CALLED PCHDF, PCHST, XERMSG
C***REVISION HISTORY (YYMMDD)
C 820218 DATE WRITTEN
C 820805 Converted to SLATEC library version.
C 870707 Minor corrections made to prologue..
C 890411 Added SAVE statements (Vers. 3.2).
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 900328 Added TYPE section. (WRB)
C 910408 Updated AUTHOR section in prologue. (WRB)
C 930503 Improved purpose. (FNF)
C***END PROLOGUE PCHCE
C
C Programming notes:
C 1. The function PCHST(ARG1,ARG2) is assumed to return zero if
C either argument is zero, +1 if they are of the same sign, and
C -1 if they are of opposite sign.
C 2. One could reduce the number of arguments and amount of local
C storage, at the expense of reduced code clarity, by passing in
C the array WK (rather than splitting it into H and SLOPE) and
C increasing its length enough to incorporate STEMP and XTEMP.
C 3. The two monotonicity checks only use the sufficient conditions.
C Thus, it is possible (but unlikely) for a boundary condition to
C be changed, even though the original interpolant was monotonic.
C (At least the result is a continuous function of the data.)
C**End
C
C DECLARE ARGUMENTS.
C
INTEGER IC(2), N, INCFD, IERR
REAL VC(2), X(*), H(*), SLOPE(*), D(INCFD,*)
C
C DECLARE LOCAL VARIABLES.
C
INTEGER IBEG, IEND, IERF, INDEX, J, K
REAL HALF, STEMP(3), THREE, TWO, XTEMP(4), ZERO
SAVE ZERO, HALF, TWO, THREE
REAL PCHDF, PCHST
C
C INITIALIZE.
C
DATA ZERO /0./, HALF /0.5/, TWO /2./, THREE /3./
C
C***FIRST EXECUTABLE STATEMENT PCHCE
IBEG = IC(1)
IEND = IC(2)
IERR = 0
C
C SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
C
IF ( ABS(IBEG).GT.N ) IBEG = 0
IF ( ABS(IEND).GT.N ) IEND = 0
C
C TREAT BEGINNING BOUNDARY CONDITION.
C
IF (IBEG .EQ. 0) GO TO 2000
K = ABS(IBEG)
IF (K .EQ. 1) THEN
C BOUNDARY VALUE PROVIDED.
D(1,1) = VC(1)
ELSE IF (K .EQ. 2) THEN
C BOUNDARY SECOND DERIVATIVE PROVIDED.
D(1,1) = HALF*( (THREE*SLOPE(1) - D(1,2)) - HALF*VC(1)*H(1) )
ELSE IF (K .LT. 5) THEN
C USE K-POINT DERIVATIVE FORMULA.
C PICK UP FIRST K POINTS, IN REVERSE ORDER.
DO 10 J = 1, K
INDEX = K-J+1
C INDEX RUNS FROM K DOWN TO 1.
XTEMP(J) = X(INDEX)
IF (J .LT. K) STEMP(J) = SLOPE(INDEX-1)
10 CONTINUE
C -----------------------------
D(1,1) = PCHDF (K, XTEMP, STEMP, IERF)
C -----------------------------
IF (IERF .NE. 0) GO TO 5001
ELSE
C USE 'NOT A KNOT' CONDITION.
D(1,1) = ( THREE*(H(1)*SLOPE(2) + H(2)*SLOPE(1))
* - TWO*(H(1)+H(2))*D(1,2) - H(1)*D(1,3) ) / H(2)
ENDIF
C
IF (IBEG .GT. 0) GO TO 2000
C
C CHECK D(1,1) FOR COMPATIBILITY WITH MONOTONICITY.
C
IF (SLOPE(1) .EQ. ZERO) THEN
IF (D(1,1) .NE. ZERO) THEN
D(1,1) = ZERO
IERR = IERR + 1
ENDIF
ELSE IF ( PCHST(D(1,1),SLOPE(1)) .LT. ZERO) THEN
D(1,1) = ZERO
IERR = IERR + 1
ELSE IF ( ABS(D(1,1)) .GT. THREE*ABS(SLOPE(1)) ) THEN
D(1,1) = THREE*SLOPE(1)
IERR = IERR + 1
ENDIF
C
C TREAT END BOUNDARY CONDITION.
C
2000 CONTINUE
IF (IEND .EQ. 0) GO TO 5000
K = ABS(IEND)
IF (K .EQ. 1) THEN
C BOUNDARY VALUE PROVIDED.
D(1,N) = VC(2)
ELSE IF (K .EQ. 2) THEN
C BOUNDARY SECOND DERIVATIVE PROVIDED.
D(1,N) = HALF*( (THREE*SLOPE(N-1) - D(1,N-1)) +
* HALF*VC(2)*H(N-1) )
ELSE IF (K .LT. 5) THEN
C USE K-POINT DERIVATIVE FORMULA.
C PICK UP LAST K POINTS.
DO 2010 J = 1, K
INDEX = N-K+J
C INDEX RUNS FROM N+1-K UP TO N.
XTEMP(J) = X(INDEX)
IF (J .LT. K) STEMP(J) = SLOPE(INDEX)
2010 CONTINUE
C -----------------------------
D(1,N) = PCHDF (K, XTEMP, STEMP, IERF)
C -----------------------------
IF (IERF .NE. 0) GO TO 5001
ELSE
C USE 'NOT A KNOT' CONDITION.
D(1,N) = ( THREE*(H(N-1)*SLOPE(N-2) + H(N-2)*SLOPE(N-1))
* - TWO*(H(N-1)+H(N-2))*D(1,N-1) - H(N-1)*D(1,N-2) )
* / H(N-2)
ENDIF
C
IF (IEND .GT. 0) GO TO 5000
C
C CHECK D(1,N) FOR COMPATIBILITY WITH MONOTONICITY.
C
IF (SLOPE(N-1) .EQ. ZERO) THEN
IF (D(1,N) .NE. ZERO) THEN
D(1,N) = ZERO
IERR = IERR + 2
ENDIF
ELSE IF ( PCHST(D(1,N),SLOPE(N-1)) .LT. ZERO) THEN
D(1,N) = ZERO
IERR = IERR + 2
ELSE IF ( ABS(D(1,N)) .GT. THREE*ABS(SLOPE(N-1)) ) THEN
D(1,N) = THREE*SLOPE(N-1)
IERR = IERR + 2
ENDIF
C
C NORMAL RETURN.
C
5000 CONTINUE
RETURN
C
C ERROR RETURN.
C
5001 CONTINUE
C ERROR RETURN FROM PCHDF.
C *** THIS CASE SHOULD NEVER OCCUR ***
IERR = -1
CALL XERMSG ('SLATEC', 'PCHCE', 'ERROR RETURN FROM PCHDF', IERR,
+ 1)
RETURN
C------------- LAST LINE OF PCHCE FOLLOWS ------------------------------
END
|
