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 46269 : 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 46269 : real(dp), dimension(:), pointer :: h, s, p
48 : integer :: i
49 : logical, parameter :: dbg = .true.
50 46269 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
51 46269 : f(1:4,1:nx) => f1(1:4*nx)
52 :
53 : include 'formats'
54 :
55 46269 : ierr = 0
56 :
57 46269 : if (nx < 2) then
58 17 : return
59 : end if
60 :
61 46269 : if (nx == 2) then
62 1 : call mk_pmlinear(x, f1, slope_only, nwork, work1, str, ierr)
63 1 : return
64 : end if
65 :
66 46268 : if (nx == 3) then
67 7 : call mk_pmquad(x, f1, slope_only, nwork, work1, str, ierr)
68 7 : return
69 : end if
70 :
71 46261 : if (nx == 4) then
72 9 : call mk_pmcub4(x, f1, slope_only, nwork, work1, str, ierr)
73 9 : return
74 : end if
75 :
76 46252 : if (nwork < pm_work_size) then
77 0 : ierr = -1
78 0 : return
79 : end if
80 :
81 46252 : h(1:nx) => work1(1:nx)
82 46252 : s(1:nx) => work1(1+nx:2*nx)
83 46252 : p(1:nx) => work1(1+2*nx:3*nx)
84 :
85 : if (dbg) then
86 1911484 : h(:) = 0; s(:) = 0; p(:) = 0
87 : end if
88 :
89 621744 : do i=1,nx-1
90 621744 : h(i) = x(i+1) - x(i) ! width of interval
91 : end do
92 621744 : do i = 1, nx-1
93 621744 : 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 621744 : do i=1,nx-1
102 621744 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
103 : end do
104 :
105 575492 : do i=2,nx-1
106 575492 : 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 575492 : do i=2,nx-1
110 529240 : f(2,i) = (dsign(1d0, s(i-1))+dsign(1d0, s(i)))* &
111 1104732 : 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 46252 : 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 46252 : if (p(1)*s(1) <= 0) then
118 36660 : f(2, 1) = 0
119 9592 : else if (abs(p(1)) > 2*abs(s(1))) then
120 144 : f(2, 1) = 2*s(1)
121 : else
122 9448 : 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 46252 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
127 : ! slope at nx of parabola through last 3 points
128 46252 : if (p(nx)*s(nx-1) <= 0) then
129 31953 : f(2,nx) = 0
130 14299 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
131 660 : f(2,nx) = 2*s(nx-1)
132 : else
133 13639 : f(2,nx) = p(nx)
134 : end if
135 :
136 46252 : if (slope_only) return
137 :
138 621744 : do i=1,nx-1
139 575492 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
140 621744 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
141 : end do
142 46252 : f(3,nx) = 0
143 46252 : f(4,nx) = 0
144 :
145 46269 : end subroutine mk_pmcub
146 :
147 :
148 : ! optimize special case for nx = 4
149 9 : 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 9 : real(dp), dimension(:), pointer :: h, s, p
161 : integer :: i
162 9 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
163 9 : f(1:4,1:nx) => f1(1:4*nx)
164 :
165 9 : ierr = 0
166 :
167 9 : if (nwork < pm_work_size) then
168 0 : ierr = -1
169 0 : return
170 : end if
171 :
172 9 : h(1:nx) => work1(1:nx)
173 9 : s(1:nx) => work1(1+nx:2*nx)
174 9 : p(1:nx) => work1(1+2*nx:3*nx)
175 :
176 36 : do i=1,nx-1
177 36 : h(i) = x(i+1) - x(i) ! width of interval
178 : end do
179 36 : do i = 1, nx-1
180 36 : 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 36 : do i=1,nx-1
189 36 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
190 : end do
191 27 : do i=2,nx-1
192 27 : 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 27 : do i=2,nx-1
196 0 : f(2,i) = (dsign(1d0, s(i-1))+dsign(1d0, s(i)))* &
197 27 : 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 9 : 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 9 : if (p(1)*s(1) <= 0) then
204 0 : f(2, 1) = 0
205 9 : else if (abs(p(1)) > 2*abs(s(1))) then
206 0 : f(2, 1) = 2*s(1)
207 : else
208 9 : 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 9 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
213 : ! slope at nx of parabola through last 3 points
214 9 : if (p(nx)*s(nx-1) <= 0) then
215 1 : f(2, nx) = 0
216 8 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
217 0 : f(2, nx) = 2*s(nx-1)
218 : else
219 8 : f(2, nx) = p(nx)
220 : end if
221 :
222 9 : if (slope_only) return
223 :
224 36 : do i=1,nx-1
225 27 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
226 36 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
227 : end do
228 9 : f(3,nx) = 0
229 9 : f(4,nx) = 0
230 :
231 9 : end subroutine mk_pmcub4
232 :
233 :
234 : ! optimize special case for nx = 3
235 7 : 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 7 : real(dp), dimension(:), pointer :: h, s, p
248 : integer :: i
249 7 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
250 7 : f(1:4,1:nx) => f1(1:4*nx)
251 :
252 7 : if (nwork < pm_work_size) then
253 0 : ierr = -1
254 0 : return
255 : end if
256 7 : ierr = 0
257 :
258 7 : h(1:nx) => work1(1:nx)
259 7 : s(1:nx) => work1(1+nx:2*nx)
260 7 : p(1:nx) => work1(1+2*nx:3*nx)
261 :
262 21 : do i=1,nx-1
263 21 : h(i) = x(i+1) - x(i) ! width of interval
264 : end do
265 21 : do i = 1, nx-1
266 21 : 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 21 : do i=1,nx-1
275 21 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
276 : end do
277 14 : do i=2,nx-1
278 14 : 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 14 : do i=2,nx-1
282 0 : f(2,i) = (dsign(1d0, s(i-1))+dsign(1d0, s(i)))* &
283 14 : 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 7 : 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 7 : if (p(1)*s(1) <= 0) then
290 1 : f(2, 1) = 0
291 6 : else if (abs(p(1)) > 2*abs(s(1))) then
292 1 : f(2, 1) = 2*s(1)
293 : else
294 5 : 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 7 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
299 : ! slope at nx of parabola through last 3 points
300 7 : if (p(nx)*s(nx-1) <= 0) then
301 1 : f(2, nx) = 0
302 6 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
303 0 : f(2, nx) = 2*s(nx-1)
304 : else
305 6 : f(2, nx) = p(nx)
306 : end if
307 :
308 7 : if (slope_only) return
309 :
310 21 : do i=1,nx-1
311 14 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
312 21 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
313 : end do
314 7 : f(3,nx) = 0
315 7 : f(4,nx) = 0
316 :
317 7 : end subroutine mk_pmquad
318 :
319 :
320 : ! optimize special case for nx = 2
321 1 : 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 : real(dp) :: h, s
333 1 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
334 1 : f(1:4,1:nx) => f1(1:4*nx)
335 :
336 1 : ierr = 0
337 :
338 1 : if (nwork < pm_work_size) then
339 0 : ierr = -1
340 0 : return
341 : end if
342 :
343 1 : h = x(2) - x(1) ! width of interval
344 1 : 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 1 : s = (f(1, 2) - f(1, 1)) / h ! slope across interval
352 1 : f(2, 1) = s
353 1 : f(2, 2) = 0
354 :
355 1 : if (slope_only) return
356 :
357 3 : f(3, 1:2) = 0
358 3 : f(4, 1:2) = 0
359 :
360 1 : 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 : 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 0 : 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
|
