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