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
|
*DECK DPCHCM
SUBROUTINE DPCHCM (N, X, F, D, INCFD, SKIP, ISMON, IERR)
C***BEGIN PROLOGUE DPCHCM
C***PURPOSE Check a cubic Hermite function for monotonicity.
C***LIBRARY SLATEC (PCHIP)
C***CATEGORY E3
C***TYPE DOUBLE PRECISION (PCHCM-S, DPCHCM-D)
C***KEYWORDS CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
C PCHIP, PIECEWISE CUBIC INTERPOLATION, UTILITY ROUTINE
C***AUTHOR Fritsch, F. N., (LLNL)
C Computing & 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 PARAMETER (INCFD = ...)
C INTEGER N, ISMON(N), IERR
C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N)
C LOGICAL SKIP
C
C CALL DPCHCM (N, X, F, D, INCFD, SKIP, ISMON, IERR)
C
C *Arguments:
C
C N:IN is the number of data points. (Error return if N.LT.2 .)
C
C X:IN is a real*8 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.
C (Error return if not.)
C
C F:IN is a real*8 array of function values. F(1+(I-1)*INCFD) is
C the value corresponding to X(I).
C
C D:IN is a real*8 array of derivative values. D(1+(I-1)*INCFD) is
C is the value corresponding to X(I).
C
C INCFD:IN is the increment between successive values in F and D.
C (Error return if INCFD.LT.1 .)
C
C SKIP:INOUT is a logical variable which should be set to
C .TRUE. if the user wishes to skip checks for validity of
C preceding parameters, or to .FALSE. otherwise.
C This will save time in case these checks have already
C been performed.
C SKIP will be set to .TRUE. on normal return.
C
C ISMON:OUT is an integer array indicating on which intervals the
C PCH function defined by N, X, F, D is monotonic.
C For data interval [X(I),X(I+1)],
C ISMON(I) = -3 if function is probably decreasing;
C ISMON(I) = -1 if function is strictly decreasing;
C ISMON(I) = 0 if function is constant;
C ISMON(I) = 1 if function is strictly increasing;
C ISMON(I) = 2 if function is non-monotonic;
C ISMON(I) = 3 if function is probably increasing.
C If ABS(ISMON)=3, this means that the D-values are near
C the boundary of the monotonicity region. A small
C increase produces non-monotonicity; decrease, strict
C monotonicity.
C The above applies to I=1(1)N-1. ISMON(N) indicates whether
C the entire function is monotonic on [X(1),X(N)].
C
C IERR:OUT is an error flag.
C Normal return:
C IERR = 0 (no errors).
C "Recoverable" errors:
C IERR = -1 if N.LT.2 .
C IERR = -2 if INCFD.LT.1 .
C IERR = -3 if the X-array is not strictly increasing.
C (The ISMON-array has not been changed in any of these cases.)
C NOTE: The above errors are checked in the order listed,
C and following arguments have **NOT** been validated.
C
C *Description:
C
C DPCHCM: Piecewise Cubic Hermite -- Check Monotonicity.
C
C Checks the piecewise cubic Hermite function defined by N,X,F,D
C for monotonicity.
C
C To provide compatibility with DPCHIM and DPCHIC, includes an
C increment between successive values of the F- and D-arrays.
C
C *Cautions:
C This provides the same capability as old DPCHMC, except that a
C new output value, -3, was added February 1989. (Formerly, -3
C and +3 were lumped together in the single value 3.) Codes that
C flag nonmonotonicity by "IF (ISMON.EQ.2)" need not be changed.
C Codes that check via "IF (ISMON.GE.3)" should change the test to
C "IF (IABS(ISMON).GE.3)". Codes that declare monotonicity via
C "IF (ISMON.LE.1)" should change to "IF (IABS(ISMON).LE.1)".
C
C***REFERENCES F. N. Fritsch and R. E. Carlson, Monotone piecewise
C cubic interpolation, SIAM Journal on Numerical Ana-
C lysis 17, 2 (April 1980), pp. 238-246.
C***ROUTINES CALLED DCHFCM, XERMSG
C***REVISION HISTORY (YYMMDD)
C 820518 DATE WRITTEN
C 820804 Converted to SLATEC library version.
C 831201 Reversed order of subscripts of F and D, so that the
C routine will work properly when INCFD.GT.1 . (Bug!)
C 870707 Corrected XERROR calls for d.p. name(s).
C 890206 Corrected XERROR calls.
C 890209 Added possible ISMON value of -3 and modified code so
C that 1,3,-1 produces ISMON(N)=2, rather than 3.
C 890306 Added caution about changed output.
C 890407 Changed name from DPCHMC to DPCHCM, as requested at the
C March 1989 SLATEC CML meeting, and made a few other
C minor modifications necessitated by this change.
C 890407 Converted to new SLATEC format.
C 890407 Modified DESCRIPTION to LDOC format.
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 920429 Revised format and order of references. (WRB,FNF)
C***END PROLOGUE DPCHCM
C
C Fortran intrinsics used: ISIGN.
C Other routines used: CHFCM, XERMSG.
C
C ----------------------------------------------------------------------
C
C Programming notes:
C
C An alternate organization would have separate loops for computing
C ISMON(i), i=1,...,NSEG, and for the computation of ISMON(N). The
C first loop can be readily parallelized, since the NSEG calls to
C CHFCM are independent. The second loop can be cut short if
C ISMON(N) is ever equal to 2, for it cannot be changed further.
C
C To produce a single precision version, simply:
C a. Change DPCHCM to PCHCM wherever it occurs,
C b. Change DCHFCM to CHFCM wherever it occurs, and
C c. Change the double precision declarations to real.
C
C DECLARE ARGUMENTS.
C
INTEGER N, INCFD, ISMON(N), IERR
DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N)
LOGICAL SKIP
C
C DECLARE LOCAL VARIABLES.
C
INTEGER I, NSEG
DOUBLE PRECISION DELTA
INTEGER DCHFCM
C
C VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT DPCHCM
IF (SKIP) GO TO 5
C
IF ( N.LT.2 ) GO TO 5001
IF ( INCFD.LT.1 ) GO TO 5002
DO 1 I = 2, N
IF ( X(I).LE.X(I-1) ) GO TO 5003
1 CONTINUE
SKIP = .TRUE.
C
C FUNCTION DEFINITION IS OK -- GO ON.
C
5 CONTINUE
NSEG = N - 1
DO 90 I = 1, NSEG
DELTA = (F(1,I+1)-F(1,I))/(X(I+1)-X(I))
C -------------------------------
ISMON(I) = DCHFCM (D(1,I), D(1,I+1), DELTA)
C -------------------------------
IF (I .EQ. 1) THEN
ISMON(N) = ISMON(1)
ELSE
C Need to figure out cumulative monotonicity from following
C "multiplication table":
C
C + I S M O N (I)
C + -3 -1 0 1 3 2
C +------------------------+
C I -3 I -3 -3 -3 2 2 2 I
C S -1 I -3 -1 -1 2 2 2 I
C M 0 I -3 -1 0 1 3 2 I
C O 1 I 2 2 1 1 3 2 I
C N 3 I 2 2 3 3 3 2 I
C (N) 2 I 2 2 2 2 2 2 I
C +------------------------+
C Note that the 2 row and column are out of order so as not
C to obscure the symmetry in the rest of the table.
C
C No change needed if equal or constant on this interval or
C already declared nonmonotonic.
IF ( (ISMON(I).NE.ISMON(N)) .AND. (ISMON(I).NE.0)
. .AND. (ISMON(N).NE.2) ) THEN
IF ( (ISMON(I).EQ.2) .OR. (ISMON(N).EQ.0) ) THEN
ISMON(N) = ISMON(I)
ELSE IF (ISMON(I)*ISMON(N) .LT. 0) THEN
C This interval has opposite sense from curve so far.
ISMON(N) = 2
ELSE
C At this point, both are nonzero with same sign, and
C we have already eliminated case both +-1.
ISMON(N) = ISIGN (3, ISMON(N))
ENDIF
ENDIF
ENDIF
90 CONTINUE
C
C NORMAL RETURN.
C
IERR = 0
RETURN
C
C ERROR RETURNS.
C
5001 CONTINUE
C N.LT.2 RETURN.
IERR = -1
CALL XERMSG ('SLATEC', 'DPCHCM',
+ 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
RETURN
C
5002 CONTINUE
C INCFD.LT.1 RETURN.
IERR = -2
CALL XERMSG ('SLATEC', 'DPCHCM', 'INCREMENT LESS THAN ONE', IERR,
+ 1)
RETURN
C
5003 CONTINUE
C X-ARRAY NOT STRICTLY INCREASING.
IERR = -3
CALL XERMSG ('SLATEC', 'DPCHCM',
+ 'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
RETURN
C------------- LAST LINE OF DPCHCM FOLLOWS -----------------------------
END
|
