Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module interp_1d_pm ! piecewise monotonic algorithms
21 :
22 : use const_def, only: dp
23 :
24 : implicit none
25 :
26 : contains
27 :
28 : ! the following produce piecewise monotonic interpolants
29 : ! rather than monotonicity preserving
30 : ! this stricter limit never introduces interpolated values exceeding the given values,
31 : ! even in places where the given values are not monotonic.
32 : ! the downside is reduced accuracy on smooth data compared to the mp routines.
33 :
34 :
35 76990 : subroutine mk_pmcub(x, nx, f1, slope_only, nwork, work1, str, ierr)
36 : ! make piecewise monotonic cubic interpolant
37 : use interp_1d_def
38 : integer, intent(in) :: nx ! length of x vector (nx >= 2)
39 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
40 : real(dp), intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
41 : logical, intent(in) :: slope_only
42 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
43 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
44 : character (len=*) :: str
45 : integer, intent(out) :: ierr
46 :
47 76990 : real(dp), dimension(:), pointer :: h, s, p
48 : integer :: i
49 : logical, parameter :: dbg = .true.
50 76990 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
51 76990 : f(1:4,1:nx) => f1(1:4*nx)
52 :
53 : include 'formats'
54 :
55 76990 : ierr = 0
56 :
57 76990 : if (nx < 2) then
58 27 : return
59 : end if
60 :
61 76990 : if (nx == 2) then
62 2 : call mk_pmlinear(x, f1, slope_only, nwork, work1, str, ierr)
63 2 : return
64 : end if
65 :
66 76988 : if (nx == 3) then
67 13 : call mk_pmquad(x, f1, slope_only, nwork, work1, str, ierr)
68 13 : return
69 : end if
70 :
71 76975 : if (nx == 4) then
72 12 : call mk_pmcub4(x, f1, slope_only, nwork, work1, str, ierr)
73 12 : return
74 : end if
75 :
76 76963 : if (nwork < pm_work_size) then
77 0 : ierr = -1
78 0 : return
79 : end if
80 :
81 76963 : h(1:nx) => work1(1:nx)
82 76963 : s(1:nx) => work1(1+nx:2*nx)
83 76963 : p(1:nx) => work1(1+2*nx:3*nx)
84 :
85 : if (dbg) then
86 3075808 : h(:) = 0; s(:) = 0; p(:) = 0
87 : end if
88 :
89 999615 : do i=1,nx-1
90 999615 : h(i) = x(i+1) - x(i) ! width of interval
91 : end do
92 999615 : do i = 1, nx-1
93 999615 : if (h(i) == 0) then
94 : write(*, '(a,1x,2i5,1x,a)') &
95 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
96 0 : ierr = -1
97 0 : return
98 : end if
99 : end do
100 :
101 999615 : do i=1,nx-1
102 999615 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
103 : end do
104 :
105 922652 : do i=2,nx-1
106 922652 : p(i) = (s(i-1)*h(i) + s(i)*h(i-1))/(h(i-1)+h(i))
107 : ! slope at i of parabola through i-1, i, and i+1
108 : end do
109 922652 : do i=2,nx-1
110 : f(2,i) = (dsign(1d0, s(i-1))+dsign(1d0, s(i)))* &
111 922652 : min(abs(s(i-1)), abs(s(i)), 0.5d0*abs(p(i)))
112 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
113 : end do
114 :
115 76963 : p(1) = s(1)*(1 + h(1) / (h(1) + h(2))) - s(2) * h(1) / (h(1) + h(2))
116 : ! slope at 1 of parabola through 1st 3 points
117 76963 : if (p(1)*s(1) <= 0) then
118 61026 : f(2, 1) = 0
119 15937 : else if (abs(p(1)) > 2*abs(s(1))) then
120 240 : f(2, 1) = 2*s(1)
121 : else
122 15697 : f(2, 1) = p(1)
123 : end if
124 :
125 : p(nx) = s(nx-1)*(1 + h(nx-1) / (h(nx-1) + h(nx-2))) &
126 76963 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
127 : ! slope at nx of parabola through last 3 points
128 76963 : if (p(nx)*s(nx-1) <= 0) then
129 53231 : f(2,nx) = 0
130 23732 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
131 1100 : f(2,nx) = 2*s(nx-1)
132 : else
133 22632 : f(2,nx) = p(nx)
134 : end if
135 :
136 76963 : if (slope_only) return
137 :
138 999615 : do i=1,nx-1
139 922652 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
140 999615 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
141 : end do
142 76963 : f(3,nx) = 0
143 76963 : f(4,nx) = 0
144 :
145 76990 : end subroutine mk_pmcub
146 :
147 :
148 : ! optimize special case for nx = 4
149 12 : subroutine mk_pmcub4(x, f1, slope_only, nwork, work1, str, ierr)
150 : use interp_1d_def
151 : integer, parameter :: nx = 4 ! length of x vector (nx >= 2)
152 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
153 : real(dp), intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
154 : logical, intent(in) :: slope_only
155 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
156 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
157 : character (len=*) :: str
158 : integer, intent(out) :: ierr
159 :
160 12 : real(dp), dimension(:), pointer :: h, s, p
161 : integer :: i
162 12 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
163 12 : f(1:4,1:nx) => f1(1:4*nx)
164 :
165 12 : ierr = 0
166 :
167 12 : if (nwork < pm_work_size) then
168 0 : ierr = -1
169 0 : return
170 : end if
171 :
172 12 : h(1:nx) => work1(1:nx)
173 12 : s(1:nx) => work1(1+nx:2*nx)
174 12 : p(1:nx) => work1(1+2*nx:3*nx)
175 :
176 48 : do i=1,nx-1
177 48 : h(i) = x(i+1) - x(i) ! width of interval
178 : end do
179 48 : do i = 1, nx-1
180 48 : if (h(i) == 0) then
181 : write(*, '(a,1x,2i5,1x,a)') &
182 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
183 0 : ierr = -1
184 0 : return
185 : end if
186 : end do
187 :
188 48 : do i=1,nx-1
189 48 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
190 : end do
191 36 : do i=2,nx-1
192 36 : p(i) = (s(i-1)*h(i) + s(i)*h(i-1))/(h(i-1)+h(i))
193 : ! slope at i of parabola through i-1, i, and i+1
194 : end do
195 36 : do i=2,nx-1
196 : f(2,i) = (dsign(1d0, s(i-1))+dsign(1d0, s(i)))* &
197 36 : min(abs(s(i-1)), abs(s(i)), 0.5d0*abs(p(i)))
198 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
199 : end do
200 :
201 12 : p(1) = s(1)*(1 + h(1) / (h(1) + h(2))) - s(2) * h(1) / (h(1) + h(2))
202 : ! slope at 1 of parabola through 1st 3 points
203 12 : if (p(1)*s(1) <= 0) then
204 0 : f(2, 1) = 0
205 12 : else if (abs(p(1)) > 2*abs(s(1))) then
206 0 : f(2, 1) = 2*s(1)
207 : else
208 12 : f(2, 1) = p(1)
209 : end if
210 :
211 : p(nx) = s(nx-1)*(1 + h(nx-1) / (h(nx-1) + h(nx-2))) &
212 12 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
213 : ! slope at nx of parabola through last 3 points
214 12 : if (p(nx)*s(nx-1) <= 0) then
215 2 : f(2, nx) = 0
216 10 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
217 0 : f(2, nx) = 2*s(nx-1)
218 : else
219 10 : f(2, nx) = p(nx)
220 : end if
221 :
222 12 : if (slope_only) return
223 :
224 48 : do i=1,nx-1
225 36 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
226 48 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
227 : end do
228 12 : f(3,nx) = 0
229 12 : f(4,nx) = 0
230 :
231 12 : end subroutine mk_pmcub4
232 :
233 :
234 : ! optimize special case for nx = 3
235 13 : subroutine mk_pmquad(x, f1, slope_only, nwork, work1, str, ierr)
236 : ! make piecewise monotonic quadratic interpolant
237 : use interp_1d_def
238 : integer, parameter :: nx = 3
239 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
240 : real(dp), intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
241 : logical, intent(in) :: slope_only
242 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
243 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
244 : character (len=*) :: str
245 : integer, intent(out) :: ierr
246 :
247 13 : real(dp), dimension(:), pointer :: h, s, p
248 : integer :: i
249 13 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
250 13 : f(1:4,1:nx) => f1(1:4*nx)
251 :
252 13 : if (nwork < pm_work_size) then
253 0 : ierr = -1
254 0 : return
255 : end if
256 13 : ierr = 0
257 :
258 13 : h(1:nx) => work1(1:nx)
259 13 : s(1:nx) => work1(1+nx:2*nx)
260 13 : p(1:nx) => work1(1+2*nx:3*nx)
261 :
262 39 : do i=1,nx-1
263 39 : h(i) = x(i+1) - x(i) ! width of interval
264 : end do
265 39 : do i = 1, nx-1
266 39 : if (h(i) == 0) then
267 : write(*, '(a,1x,2i5,1x,a)') &
268 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
269 0 : ierr = -1
270 0 : return
271 : end if
272 : end do
273 :
274 39 : do i=1,nx-1
275 39 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
276 : end do
277 26 : do i=2,nx-1
278 26 : p(i) = (s(i-1)*h(i) + s(i)*h(i-1))/(h(i-1)+h(i))
279 : ! slope at i of parabola through i-1, i, and i+1
280 : end do
281 26 : do i=2,nx-1
282 : f(2,i) = (dsign(1d0, s(i-1))+dsign(1d0, s(i)))* &
283 26 : min(abs(s(i-1)), abs(s(i)), 0.5d0*abs(p(i)))
284 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
285 : end do
286 :
287 13 : p(1) = s(1)*(1 + h(1) / (h(1) + h(2))) - s(2) * h(1) / (h(1) + h(2))
288 : ! slope at 1 of parabola through 1st 3 points
289 13 : if (p(1)*s(1) <= 0) then
290 2 : f(2, 1) = 0
291 11 : else if (abs(p(1)) > 2*abs(s(1))) then
292 2 : f(2, 1) = 2*s(1)
293 : else
294 9 : f(2, 1) = p(1)
295 : end if
296 :
297 : p(nx) = s(nx-1)*(1 + h(nx-1) / (h(nx-1) + h(nx-2))) &
298 13 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
299 : ! slope at nx of parabola through last 3 points
300 13 : if (p(nx)*s(nx-1) <= 0) then
301 2 : f(2, nx) = 0
302 11 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
303 0 : f(2, nx) = 2*s(nx-1)
304 : else
305 11 : f(2, nx) = p(nx)
306 : end if
307 :
308 13 : if (slope_only) return
309 :
310 39 : do i=1,nx-1
311 26 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
312 39 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
313 : end do
314 13 : f(3,nx) = 0
315 13 : f(4,nx) = 0
316 :
317 13 : end subroutine mk_pmquad
318 :
319 :
320 : ! optimize special case for nx = 2
321 2 : subroutine mk_pmlinear(x, f1, slope_only, nwork, work1, str, ierr)
322 : use interp_1d_def
323 : integer, parameter :: nx = 2
324 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
325 : real(dp), intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
326 : logical, intent(in) :: slope_only
327 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
328 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
329 : character (len=*) :: str ! for debugging
330 : integer, intent(out) :: ierr
331 :
332 2 : real(dp) :: h, s
333 2 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
334 2 : f(1:4,1:nx) => f1(1:4*nx)
335 :
336 2 : ierr = 0
337 :
338 2 : if (nwork < pm_work_size) then
339 0 : ierr = -1
340 0 : return
341 : end if
342 :
343 2 : h = x(2) - x(1) ! width of interval
344 2 : if (h == 0) then
345 0 : ierr = -1
346 : write(*, '(a,1x,2i5,1x,a)') &
347 0 : 'same interpolation x values at', 1, 2, 'for ' // trim(str)
348 0 : return
349 : end if
350 :
351 2 : s = (f(1, 2) - f(1, 1)) / h ! slope across interval
352 2 : f(2, 1) = s
353 2 : f(2, 2) = 0
354 :
355 2 : if (slope_only) return
356 :
357 6 : f(3, 1:2) = 0
358 6 : f(4, 1:2) = 0
359 :
360 2 : end subroutine mk_pmlinear
361 :
362 :
363 0 : subroutine mk_pmcub_uniform(dx, n, f1, slope_only, nwork, work1, str, ierr)
364 : ! make piecewise monotonic cubic interpolant on unit spaced mesh
365 : use interp_1d_def
366 : integer, intent(in) :: n ! length of vector
367 : real(dp), intent(in) :: dx
368 : real(dp), intent(inout), pointer :: f1(:) ! =(4,n) ! data & interpolation coefficients
369 : logical, intent(in) :: slope_only
370 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
371 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
372 : character (len=*) :: str ! for debugging
373 : integer, intent(out) :: ierr
374 :
375 0 : real(dp), dimension(:), pointer :: s, p
376 0 : real(dp) :: x(2)
377 : integer :: i
378 0 : real(dp), pointer :: f(:,:) ! (4, n) ! data & interpolation coefficients
379 0 : f(1:4,1:n) => f1(1:4*n)
380 :
381 0 : ierr = 0
382 :
383 0 : if (n < 2) then
384 0 : return
385 : end if
386 :
387 0 : if (n == 2) then
388 0 : x(1) = 0
389 0 : x(2) = dx
390 0 : call mk_pmlinear(x, f1, slope_only, nwork, work1, str, ierr)
391 0 : return
392 : end if
393 :
394 0 : if (dx == 0) then
395 0 : ierr = -1
396 0 : return
397 : end if
398 :
399 0 : if (nwork < pm_work_size) then
400 0 : ierr = -1
401 0 : return
402 : end if
403 : ierr = 0
404 :
405 0 : s(1:n) => work1(1:n)
406 0 : p(1:n) => work1(1+n:2*n)
407 :
408 : ierr = 0
409 :
410 0 : do i=1,n-1
411 0 : s(i) = f(1,i+1) - f(1,i) ! slope across interval
412 : end do
413 0 : do i=2,n-1
414 0 : p(i) = 0.5d0*(s(i-1) + s(i))
415 : ! slope at i of parabola through i-1, i, and i+1
416 : end do
417 0 : do i=2,n-1
418 : f(2,i) = (sign(1d0, s(i-1))+sign(1d0, s(i)))* &
419 0 : min(abs(s(i-1)), abs(s(i)), 0.5d0*abs(p(i)))
420 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
421 : end do
422 :
423 0 : p(1) = 1.5d0 * s(1) - 0.5d0 * s(2)
424 : ! slope at 1 of parabola through 1st 3 points
425 0 : if (p(1)*s(1) <= 0) then
426 0 : f(2, 1) = 0
427 0 : else if (abs(p(1)) > 2*abs(s(1))) then
428 0 : f(2, 1) = 2*s(1)
429 : else
430 0 : f(2, 1) = p(1)
431 : end if
432 :
433 0 : p(n) = 1.5d0 * s(n-1) - 0.5d0 * s(n-2)
434 : ! slope at n of parabola through last 3 points
435 0 : if (p(n)*s(n-1) <= 0) then
436 0 : f(2, n) = 0
437 0 : else if (abs(p(n)) > 2*abs(s(n-1))) then
438 0 : f(2, n) = 2*s(n-1)
439 : else
440 0 : f(2, n) = p(n)
441 : end if
442 :
443 0 : f(2, 1:n) = f(2, 1:n) / dx
444 :
445 0 : if (slope_only) return
446 :
447 : ! 2nd and 3rd derivatives
448 0 : do i=1,n-1
449 0 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / dx
450 0 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (dx*dx)
451 : end do
452 0 : f(3,n) = (3*f(1, n-1) - 3*f(1, n) + (f(2, n-1) + 2*f(2, n)) * dx) / (dx*dx*dx)
453 0 : f(4,n) = (-2*f(1, n-1) + 2*f(1, n) - (f(2, n-1) + f(2, n))*dx) / (dx*dx*dx)
454 :
455 0 : end subroutine mk_pmcub_uniform
456 :
457 :
458 : end module interp_1d_pm
|
