Coverage Report: pchipMod.f90

Generated from GCOV analysis of Fortran source code

0.0%
Lines Executed
73 total lines
0.0%
Branches Executed
0 total branches
0.0%
Calls Executed
0 total calls
0
-
Source:pchipMod.f90
0
-
Graph:pchipMod.gcno
0
-
Data:pchipMod.gcda
0
-
Runs:29
1
-
!! Piecewise Cubic Hermite Interpolating Polynomial calculations
2
-
3
-
MODULE pchipmod
4
-
!! Subroutines and functions related to the calculation of the Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)
5
-
6
-
USE splinec
7
-
8
-
IMPLICIT NONE
9
-
PUBLIC
10
-
SAVE
11
-
12
-
REAL (KIND=8), PRIVATE :: h
13
-
REAL (KIND=8), PRIVATE :: fprime_r, fprime_i
14
-
15
-
CONTAINS
16
-
17
#####
SUBROUTINE PCHIP( x, y, N, PolyCoef, csWork )
18
-
19
-
! This implements the monotone piecewise cubic Hermite interpolating
20
-
! polynomial (PCHIP) algorithm. This is a new variant of monotone PCHIP
21
-
! (paper submitted to JASA). Also see;
22
-
!
23
-
! F. N. Fritsch and J. Butland, "A Method for Constructing Local Monotone
24
-
! Piecewise Cubic Interpolants", SIAM Journal on Scientific and Statistical
25
-
! Computing, 5(2):300-304, (1984) https://doi.org/10.1137/0905021
26
-
!
27
-
! F. N. Fritsch and R. E. Carlson. "Monotone Piecewise Cubic Interpolation",
28
-
! SIAM Journal on Numerical Analysis, 17(2):238-246, (1980) https://doi.org/10.1137/0717021
29
-
!
30
-
! N is the number of nodes
31
-
! x is a vector of the abscissa values
32
-
! y is a vector of the associated ordinate values
33
-
! PolyCoef are the coefficients of the standard polynomial
34
-
! csWork is a temporary work space for the cubic spline
35
-
36
-
INTEGER, INTENT( IN ) :: N
37
-
REAL (KIND=8), INTENT( IN ) :: x(:)
38
-
COMPLEX (KIND=8), INTENT( IN ) :: y(:)
39
-
COMPLEX (KIND=8), INTENT( INOUT ) :: PolyCoef( 4, N ), csWork( 4, N )
40
-
41
-
INTEGER :: ix, iBCBeg, iBCEnd
42
-
REAL (KIND=8) :: h1, h2
43
-
COMPLEX (KIND=8) :: del1, del2, f1, f2, f1prime, f2prime, fprimeT
44
-
45
-
! Precompute estimates of the derivatives at the nodes
46
-
! The vector PolyCoef(1,*) holds the ordinate values at the nodes
47
-
! The vector PolyCoef(2,*) holds the ordinate derivatives at the nodes
48
-
49
#####
IF ( N == 2 ) THEN
50
-
51
-
! handle special case of two data points seperately (linear interpolation)
52
-
53
#####
PolyCoef( 1, 1 ) = y( 1 )
54
#####
PolyCoef( 2, 1 ) = ( y( 2 ) - y( 1 ) ) / ( x( 2 ) - x( 1 ) )
55
#####
PolyCoef( 3, 1 ) = 0.0D0
56
#####
PolyCoef( 4, 1 ) = 0.0D0
57
-
58
-
ELSE
59
-
60
-
! general case of more than two data points
61
-
62
#####
PolyCoef( 1, 1 : N ) = y( 1 : N )
63
-
64
-
! left endpoint (non-centered 3-point difference formula)
65
-
66
#####
CALL h_del( x, y, 2, h1, h2, del1, del2 )
67
#####
fprimeT = ( ( 2.0D0 * h1 + h2 ) * del1 - h1 * del2 ) / ( h1 + h2 )
68
#####
PolyCoef( 2, 1 ) = fprime_left_end_Cmplx( del1, del2, fprimeT )
69
-
70
-
! right endpoint (non-centered 3-point difference formula)
71
-
72
#####
CALL h_del( x, y, N - 1, h1, h2, del1, del2 )
73
#####
fprimeT = ( -h2 * del1 + ( h1 + 2.0D0 * h2 ) * del2 ) / ( h1 + h2 )
74
#####
PolyCoef( 2, N ) = fprime_right_end_Cmplx( del1, del2, fprimeT )
75
-
76
-
! compute coefficients of the cubic spline interpolating polynomial
77
-
78
#####
iBCBeg = 1 ! specified derivatives at the end points
79
#####
iBCEnd = 1
80
#####
csWork( 1, 1 : N ) = PolyCoef( 1, 1 : N )
81
#####
csWork( 2, 1 ) = PolyCoef( 2, 1 )
82
#####
csWork( 2, N ) = PolyCoef( 2, N )
83
#####
CALL CSpline( x, csWork( 1, 1 ), N, iBCBeg, iBCEnd, N )
84
-
85
-
! interior nodes (use derivatives from the cubic spline as initial estimate)
86
-
87
#####
DO ix = 2, N - 1
88
#####
CALL h_del( x, y, ix, h1, h2, del1, del2 )
89
-
! check if the derivative from the cubic spline satisfies monotonicity
90
#####
PolyCoef( 2, ix ) = fprime_interior_Cmplx( del1, del2, csWork( 2, ix ) )
91
-
END DO
92
-
93
-
! 2 3
94
-
! compute coefficients of std cubic polynomial: c0 + c1*x + c2*x + c3*x
95
-
!
96
-
97
#####
DO ix = 1, N - 1
98
#####
h = x( ix + 1 ) - x( ix )
99
-
100
#####
f1 = PolyCoef( 1, ix )
101
#####
f2 = PolyCoef( 1, ix + 1 )
102
-
103
#####
f1prime = PolyCoef( 2, ix )
104
#####
f2prime = PolyCoef( 2, ix + 1 )
105
-
106
#####
PolyCoef( 3, ix ) = ( 3.0D0 * ( f2 - f1 ) - h * ( 2.0D0 * f1prime + f2prime ) ) / h**2
107
#####
PolyCoef( 4, ix ) = ( h * ( f1prime + f2prime ) - 2.0D0 * ( f2 - f1 ) ) / h**3
108
-
END DO
109
-
110
-
END IF
111
-
112
#####
RETURN
113
-
END SUBROUTINE PCHIP
114
-
115
-
! **********************************************************************!
116
-
117
#####
SUBROUTINE h_del( x, y, ix, h1, h2, del1, del2 )
118
-
119
-
INTEGER, INTENT( IN ) :: ix ! index of the center point
120
-
REAL (KIND=8), INTENT( IN ) :: x(:)
121
-
COMPLEX (KIND=8), INTENT( IN ) :: y(:)
122
-
REAL (KIND=8), INTENT( OUT ) :: h1, h2
123
-
COMPLEX (KIND=8), INTENT( OUT ) :: del1, del2
124
-
125
#####
h1 = x( ix ) - x( ix - 1 )
126
#####
h2 = x( ix + 1 ) - x( ix )
127
-
128
#####
del1 = ( y( ix ) - y( ix - 1 ) ) / h1
129
#####
del2 = ( y( ix + 1 ) - y( ix ) ) / h2
130
-
131
#####
RETURN
132
-
END SUBROUTINE h_del
133
-
134
-
! **********************************************************************!
135
-
136
#####
FUNCTION fprime_interior_Cmplx( del1, del2, fprime )
137
-
138
-
COMPLEX (KIND=8), INTENT( IN ) :: del1, del2, fprime
139
-
COMPLEX (KIND=8) :: fprime_interior_Cmplx
140
-
141
#####
fprime_r = fprime_interior( REAL( del1 ), REAL( del2 ), REAL( fprime ) )
142
#####
fprime_i = fprime_interior( AIMAG( del1 ), AIMAG( del2 ), AIMAG( fprime ) )
143
-
144
#####
fprime_interior_Cmplx = CMPLX( fprime_r, fprime_i, KIND=8 )
145
-
146
#####
END FUNCTION fprime_interior_Cmplx
147
-
148
-
! **********************************************************************!
149
-
150
#####
FUNCTION fprime_left_end_Cmplx( del1, del2, fprime )
151
-
152
-
COMPLEX (KIND=8), INTENT( IN ) :: del1, del2, fprime
153
-
COMPLEX (KIND=8) :: fprime_left_end_Cmplx
154
-
155
#####
fprime_r = fprime_left_end( REAL( del1 ), REAL( del2 ), REAL( fprime ) )
156
#####
fprime_i = fprime_left_end( AIMAG( del1 ), AIMAG( del2 ), AIMAG( fprime ) )
157
-
158
#####
fprime_left_end_Cmplx = CMPLX( fprime_r, fprime_i, KIND=8 )
159
-
160
#####
END FUNCTION fprime_left_end_Cmplx
161
-
162
-
! **********************************************************************!
163
-
164
#####
FUNCTION fprime_right_end_Cmplx( del1, del2, fprime )
165
-
166
-
COMPLEX (KIND=8), INTENT( IN ) :: del1, del2, fprime
167
-
COMPLEX (KIND=8) :: fprime_right_end_Cmplx
168
-
169
#####
fprime_r = fprime_right_end( REAL( del1 ), REAL( del2 ), REAL( fprime ) )
170
#####
fprime_i = fprime_right_end( AIMAG( del1 ), AIMAG( del2 ), AIMAG( fprime ) )
171
-
172
#####
fprime_right_end_Cmplx = CMPLX( fprime_r, fprime_i, KIND=8 )
173
-
174
#####
END FUNCTION fprime_right_end_Cmplx
175
-
176
-
! **********************************************************************!
177
-
178
#####
FUNCTION fprime_interior( del1, del2, fprime )
179
-
180
-
REAL (KIND=8), INTENT( IN ) :: del1, del2, fprime
181
-
REAL (KIND=8) :: fprime_interior
182
-
183
-
! check if derivative is within the trust region, project into it if not
184
-
185
#####
IF ( del1 * del2 > 0.0 ) THEN
186
-
! adjacent secant slopes have the same sign, enforce monotonicity
187
#####
IF ( del1 > 0.0 ) THEN
188
#####
fprime_interior = MIN( MAX(fprime, 0.0D0), 3.0D0 * MIN(del1, del2) )
189
-
ELSE
190
#####
fprime_interior = MAX( MIN(fprime, 0.0D0), 3.0D0 * MAX(del1, del2) )
191
-
END IF
192
-
ELSE
193
-
! force the interpolant to have an extrema here
194
#####
fprime_interior = 0.0D0;
195
-
END IF
196
-
197
#####
END FUNCTION fprime_interior
198
-
199
-
! **********************************************************************!
200
-
201
#####
FUNCTION fprime_left_end( del1, del2, fprime )
202
-
203
-
REAL (KIND=8), INTENT( IN ) :: del1, del2, fprime
204
-
REAL (KIND=8) :: fprime_left_end
205
-
206
#####
fprime_left_end = fprime
207
-
208
#####
IF ( del1 * fprime <= 0.0D0 ) THEN
209
-
! set derivative to zero if the sign differs from sign of secant slope
210
#####
fprime_left_end = 0.0;
211
#####
ELSE IF ( ( del1 * del2 <= 0.0D0 ) .AND. ( ABS( fprime ) > ABS( 3.0D0 * del1 ) ) ) THEN
212
-
! adjust derivative value to enforce monotonicity
213
#####
fprime_left_end = 3.0D0 * del1;
214
-
END IF
215
-
216
#####
END FUNCTION fprime_left_end
217
-
218
-
! **********************************************************************!
219
-
220
#####
FUNCTION fprime_right_end( del1, del2, fprime )
221
-
222
-
! This is essentially the same as fprime_left_end( del2, del1, fprime )
223
-
! Written separately for clarity
224
-
225
-
REAL (KIND=8), INTENT( IN ) :: del1, del2, fprime
226
-
REAL (KIND=8) :: fprime_right_end
227
-
228
#####
fprime_right_end = fprime
229
-
230
#####
IF ( del2 * fprime <= 0.0D0 ) THEN
231
-
! set derivative to zero if the sign differs from sign of secant slope
232
#####
fprime_right_end = 0.0;
233
#####
ELSE IF ( ( del1 * del2 <= 0.0D0 ) .AND. ( ABS( fprime ) > ABS( 3.0D0 * del2 ) ) ) THEN
234
-
! adjust derivative value to enforce monotonicity
235
#####
fprime_right_end = 3.0D0 * del2;
236
-
END IF
237
-
238
#####
END FUNCTION fprime_right_end
239
-
240
-
END MODULE pchipmod

Legend:

Execution Count: Number of times each line was executed. ##### indicates never executed, - indicates non-executable.