Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Simple splines
10 : !> Splines are fully specified by the interpolation points, except that
11 : !> at the ends, we have the freedom to prescribe the second derivatives.
12 : !> If we know a derivative at an end (exactly), then best is to impose that.
13 : !> Otherwise, it is better to use the "consistent" end conditions: the second
14 : !> derivative is determined such that it is smooth.
15 : !>
16 : !> High level API: spline3, spline3ders.
17 : !> Low level API: the rest of public soubroutines.
18 : !>
19 : !> Use the high level API to obtain cubic spline fit with consistent boundary
20 : !> conditions and optionally the derivatives. Use the low level API if more fine
21 : !> grained control is needed.
22 : !>
23 : !> This module is based on a code written by John E. Pask, LLNL.
24 : !> \par History
25 : !> Adapted for use in CP2K (30.12.2016,JGH)
26 : !> \author JGH
27 : ! **************************************************************************************************
28 : MODULE splines
29 :
30 : USE kinds, ONLY: dp
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'splines'
38 :
39 : PUBLIC :: spline3, spline3ders
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief ...
45 : !> \param x ...
46 : !> \param y ...
47 : !> \param xnew ...
48 : !> \return ...
49 : ! **************************************************************************************************
50 0 : FUNCTION spline3(x, y, xnew) RESULT(ynew)
51 : ! Takes the function values 'y' on the grid 'x' and returns new values 'ynew'
52 : ! at the given grid 'xnew' using cubic splines interpolation with such
53 : ! boundary conditions so that the 2nd derivative is consistent with the
54 : ! interpolating cubic.
55 : REAL(dp), INTENT(in) :: x(:), y(:), xnew(:)
56 : REAL(dp) :: ynew(SIZE(xnew))
57 :
58 : INTEGER :: i, ip
59 0 : REAL(dp) :: c(0:4, SIZE(x) - 1)
60 :
61 : ! get spline parameters: 2nd derivs at ends determined by cubic interpolation
62 0 : CALL spline3pars(x, y, [2, 2], [0._dp, 0._dp], c)
63 :
64 0 : ip = 0
65 0 : DO i = 1, SIZE(xnew)
66 0 : ip = iixmin(xnew(i), x, ip)
67 0 : ynew(i) = poly3(xnew(i), c(:, ip))
68 : END DO
69 0 : END FUNCTION
70 :
71 : ! **************************************************************************************************
72 : !> \brief ...
73 : !> \param x ...
74 : !> \param y ...
75 : !> \param xnew ...
76 : !> \param ynew ...
77 : !> \param dynew ...
78 : !> \param d2ynew ...
79 : ! **************************************************************************************************
80 0 : SUBROUTINE spline3ders(x, y, xnew, ynew, dynew, d2ynew)
81 : ! Just like 'spline', but also calculate 1st and 2nd derivatives
82 : REAL(dp), INTENT(in) :: x(:), y(:), xnew(:)
83 : REAL(dp), INTENT(out), OPTIONAL :: ynew(:), dynew(:), d2ynew(:)
84 :
85 : INTEGER :: i, ip
86 0 : REAL(dp) :: c(0:4, SIZE(x) - 1)
87 :
88 0 : CALL spline3pars(x, y, [2, 2], [0._dp, 0._dp], c)
89 :
90 0 : ip = 0
91 0 : DO i = 1, SIZE(xnew)
92 0 : ip = iixmin(xnew(i), x, ip)
93 0 : IF (PRESENT(ynew)) ynew(i) = poly3(xnew(i), c(:, ip))
94 0 : IF (PRESENT(dynew)) dynew(i) = dpoly3(xnew(i), c(:, ip))
95 0 : IF (PRESENT(d2ynew)) d2ynew(i) = d2poly3(xnew(i), c(:, ip))
96 : END DO
97 0 : END SUBROUTINE
98 :
99 : ! **************************************************************************************************
100 : !> \brief ...
101 : !> \param xi ...
102 : !> \param yi ...
103 : !> \param bctype ...
104 : !> \param bcval ...
105 : !> \param c ...
106 : ! **************************************************************************************************
107 0 : SUBROUTINE spline3pars(xi, yi, bctype, bcval, c)
108 : ! Returns parameters c defining cubic spline interpolating x-y data xi, yi, with
109 : ! boundary conditions specified by bcytpe, bcvals
110 : REAL(dp), INTENT(in) :: xi(:), yi(:)
111 : INTEGER, INTENT(in) :: bctype(2)
112 : REAL(dp), INTENT(in) :: bcval(2)
113 : REAL(dp), INTENT(out) :: c(0:, :)
114 :
115 : INTEGER :: i, i2, info, ipiv(4), j, n
116 : REAL(dp) :: Ae(4, 4), be(4), bemat(4, 1), c1, c2, c3, c4, ce(4), d2p1, d2pn, x0, xe(4), &
117 0 : ye(4), hi(SIZE(c, 2)), cs(2*SIZE(c, 2)), bs(2*SIZE(c, 2)), bmat(2*SIZE(c, 2), 1), &
118 0 : As(5, 2*SIZE(c, 2))
119 0 : INTEGER :: ipiv2(2*SIZE(c, 2))
120 :
121 : ! x values of data
122 : ! y values of data
123 : ! type of boundary condition at each end:
124 : ! bctype(1) = type at left end, bctype(2) = type at right end.
125 : ! 1 = specified 2nd derivative, 2 = 2nd derivative consistent with interpolating cubic.
126 : ! boundary condition values at each end:
127 : ! bcval(1) = value at left end, bcval(2) = value at right end
128 : ! parameters defining spline: c(i,j) = ith parameter of jth
129 : ! spline polynomial, p_j = sum_{i=1}^4 c(i,j) (x-c(0,j))^(i-1), j = 1..n-1, n = # of data pts.
130 : ! dimensions: c(0:4,1:n-1)
131 : ! spline eq. matrix -- LAPACK band form
132 : ! spline eq. rhs vector
133 : ! spline eq. solution vector
134 : ! spline intervals
135 : ! end-cubic eq. matrix
136 : ! end-cubic eq. rhs vector
137 : ! end-cubic eq. solution vector
138 : ! x,y values at ends
139 : ! 2nd derivatives at ends
140 : ! expansion center
141 : ! expansion coefficients
142 : ! number of data points
143 : ! lapack variables
144 :
145 : ! check input parameters
146 0 : IF (bctype(1) < 1 .OR. bctype(1) > 2) CALL stop_error("spline3pars error: bctype /= 1 or 2.")
147 0 : IF (bctype(2) < 1 .OR. bctype(2) > 2) CALL stop_error("spline3pars error: bctype /= 1 or 2.")
148 0 : IF (SIZE(c, 1) /= 5) CALL stop_error("spline3pars error: size(c,1) /= 5.")
149 0 : IF (SIZE(c, 2) /= SIZE(xi) - 1) CALL stop_error("spline3pars error: size(c,2) /= size(xi)-1.")
150 0 : IF (SIZE(xi) /= SIZE(yi)) CALL stop_error("spline3pars error: size(xi) /= size(yi)")
151 :
152 : ! To get rid of compiler warnings:
153 0 : d2p1 = 0
154 0 : d2pn = 0
155 :
156 : ! initializations
157 0 : n = SIZE(xi)
158 0 : DO i = 1, n - 1
159 0 : hi(i) = xi(i + 1) - xi(i)
160 : END DO
161 :
162 : ! compute interpolating-cubic 2nd derivs at ends, if required
163 : ! left end
164 0 : IF (bctype(1) == 2) THEN
165 0 : IF (n < 4) CALL stop_error("spline3pars error: n < 4")
166 0 : xe = xi(1:4)
167 0 : ye = yi(1:4)
168 0 : x0 = xe(1) ! center at end
169 0 : DO i = 1, 4
170 0 : DO j = 1, 4
171 0 : Ae(i, j) = (xe(i) - x0)**(j - 1)
172 : END DO
173 : END DO
174 0 : Ae(:, 1) = 1 ! set 0^0 = 1
175 0 : be = ye; bemat(:, 1) = be
176 0 : CALL dgesv(4, 1, Ae, 4, ipiv, bemat, 4, info)
177 0 : IF (info /= 0) CALL stop_error("spline3pars error: dgesv error.")
178 0 : ce = bemat(:, 1)
179 0 : d2p1 = 2*ce(3)
180 : END IF
181 : ! right end
182 0 : IF (bctype(2) == 2) THEN
183 0 : IF (n < 4) CALL stop_error("spline3pars error: n < 4")
184 0 : xe = xi(n - 3:n)
185 0 : ye = yi(n - 3:n)
186 0 : x0 = xe(4) ! center at end
187 0 : DO i = 1, 4
188 0 : DO j = 1, 4
189 0 : Ae(i, j) = (xe(i) - x0)**(j - 1)
190 : END DO
191 : END DO
192 0 : Ae(:, 1) = 1 ! set 0^0 = 1
193 0 : be = ye; bemat(:, 1) = be
194 0 : CALL dgesv(4, 1, Ae, 4, ipiv, bemat, 4, info)
195 0 : IF (info /= 0) CALL stop_error("spline3pars error: dgesv error.")
196 0 : ce = bemat(:, 1)
197 0 : d2pn = 2*ce(3)
198 : END IF
199 :
200 : ! set 2nd derivs at ends
201 0 : IF (bctype(1) == 1) d2p1 = bcval(1)
202 0 : IF (bctype(2) == 1) d2pn = bcval(2)
203 :
204 : ! construct spline equations -- LAPACK band form
205 : ! basis: phi1 = -(x-x_i)/h_i, phi2 = (x-x_{i+1})/h_i, phi3 = phi1^3-phi1, phi4 = phi2^3-phi2
206 : ! on interval [x_i,x_{i+1}] of length h_i = x_{i+1}-x_i
207 : !A=0 ! full matrix
208 0 : As = 0
209 : ! left end condition
210 0 : As(4, 1) = 6/hi(1)**2
211 0 : bs(1) = d2p1
212 : ! internal knot conditions
213 0 : DO i = 2, n - 1
214 0 : i2 = 2*(i - 1)
215 0 : As(5, i2 - 1) = 1/hi(i - 1)
216 0 : As(4, i2) = 2/hi(i - 1)
217 0 : As(3, i2 + 1) = 2/hi(i)
218 0 : As(2, i2 + 2) = 1/hi(i)
219 0 : As(5, i2) = 1/hi(i - 1)**2
220 0 : As(4, i2 + 1) = -1/hi(i)**2
221 0 : bs(i2) = (yi(i + 1) - yi(i))/hi(i) - (yi(i) - yi(i - 1))/hi(i - 1)
222 0 : bs(i2 + 1) = 0
223 : END DO
224 : ! right end condition
225 0 : As(4, 2*(n - 1)) = 6/hi(n - 1)**2
226 0 : bs(2*(n - 1)) = d2pn
227 :
228 : ! solve spline equations -- LAPACK band form
229 0 : bmat(:, 1) = bs
230 0 : CALL dgbsv(2*(n - 1), 1, 2, 1, As, 5, ipiv2, bmat, 2*(n - 1), info)
231 0 : IF (info /= 0) CALL stop_error("spline3pars error: dgbsv error.")
232 0 : cs = bmat(:, 1)
233 :
234 : ! transform to (x-x0)^(i-1) basis and return
235 0 : DO i = 1, n - 1
236 : ! coefficients in spline basis:
237 0 : c1 = yi(i)
238 0 : c2 = yi(i + 1)
239 0 : c3 = cs(2*i - 1)
240 0 : c4 = cs(2*i)
241 : ! coefficients in (x-x0)^(i-1) basis
242 0 : c(0, i) = xi(i)
243 0 : c(1, i) = c1
244 0 : c(2, i) = -(c1 - c2 + 2*c3 + c4)/hi(i)
245 0 : c(3, i) = 3*c3/hi(i)**2
246 0 : c(4, i) = (-c3 + c4)/hi(i)**3
247 : END DO
248 0 : END SUBROUTINE
249 :
250 : !--------------------------------------------------------------------------------------------------!
251 :
252 : ! **************************************************************************************************
253 : !> \brief ...
254 : !> \param x ...
255 : !> \param xi ...
256 : !> \param c ...
257 : !> \param val ...
258 : !> \param der ...
259 : ! **************************************************************************************************
260 0 : SUBROUTINE spline3valder(x, xi, c, val, der)
261 : ! Returns value and 1st derivative of spline defined by knots xi and parameters c
262 : ! returned by spline3pars
263 : REAL(dp), INTENT(in) :: x, xi(:), c(0:, :)
264 : REAL(dp), INTENT(out) :: val, der
265 :
266 : INTEGER :: i1, n
267 :
268 : ! point at which to evaluate spline
269 : ! spline knots (x values of data)
270 : ! spline parameters: c(i,j) = ith parameter of jth
271 : ! spline polynomial, p_j = sum_{i=1}^4 c(i,j) (x-c(0,j))^(i-1), j = 1..n-1, n = # of data pts.
272 : ! dimensions: c(0:4,1:n-1)
273 : ! value of spline at x
274 : ! 1st derivative of spline at x
275 : ! number of knots
276 :
277 : ! initialize, check input parameters
278 0 : n = SIZE(xi)
279 0 : IF (SIZE(c, 1) /= 5) CALL stop_error("spline3 error: size(c,1) /= 5.")
280 0 : IF (SIZE(c, 2) /= SIZE(xi) - 1) CALL stop_error("spline3 error: size(c,2) /= size(xi)-1.")
281 : ! find interval containing x
282 0 : i1 = iix(x, xi)
283 : ! return value and derivative
284 0 : val = poly3(x, c(:, i1))
285 0 : der = dpoly3(x, c(:, i1))
286 0 : END SUBROUTINE
287 :
288 : !--------------------------------------------------------------------------------------------------!
289 :
290 : ! **************************************************************************************************
291 : !> \brief ...
292 : !> \param x ...
293 : !> \param xi ...
294 : !> \return ...
295 : ! **************************************************************************************************
296 0 : INTEGER FUNCTION iix(x, xi) RESULT(i1)
297 : ! Returns index i of interval [xi(i),xi(i+1)] containing x in mesh xi,
298 : ! with intervals indexed by left-most points.
299 : ! N.B.: x outside [x1,xn] are indexed to nearest end.
300 : ! Uses bisection, except if "x" lies in the first or second elements (which is
301 : ! often the case)
302 : REAL(dp), INTENT(in) :: x, xi(:)
303 :
304 : INTEGER :: i2, ic, n
305 :
306 : ! target value
307 : ! mesh, xi(i) < xi(i+1)
308 : ! number of mesh points
309 :
310 0 : n = SIZE(xi)
311 0 : i1 = 1
312 0 : IF (n < 2) THEN
313 0 : CALL stop_error("error in iix: n < 2")
314 0 : ELSEIF (n == 2) THEN
315 : i1 = 1
316 0 : ELSEIF (n == 3) THEN
317 0 : IF (x <= xi(2)) THEN ! first element
318 : i1 = 1
319 : ELSE
320 0 : i1 = 2
321 : END IF
322 0 : ELSEIF (x <= xi(1)) THEN ! left end
323 : i1 = 1
324 0 : ELSEIF (x <= xi(2)) THEN ! first element
325 : i1 = 1
326 0 : ELSEIF (x <= xi(3)) THEN ! second element
327 : i1 = 2
328 0 : ELSEIF (x >= xi(n)) THEN ! right end
329 0 : i1 = n - 1
330 : ELSE
331 : ! bisection: xi(i1) <= x < xi(i2)
332 : i1 = 3; i2 = n
333 : DO
334 0 : IF (i2 - i1 == 1) EXIT
335 0 : ic = i1 + (i2 - i1)/2
336 0 : IF (x >= xi(ic)) THEN
337 : i1 = ic
338 : ELSE
339 0 : i2 = ic
340 : END IF
341 : END DO
342 : END IF
343 0 : END FUNCTION
344 :
345 : ! **************************************************************************************************
346 : !> \brief ...
347 : !> \param x ...
348 : !> \param xi ...
349 : !> \param i_min ...
350 : !> \return ...
351 : ! **************************************************************************************************
352 0 : INTEGER FUNCTION iixmin(x, xi, i_min) RESULT(ip)
353 : ! Just like iix, but assumes that x >= xi(i_min)
354 : REAL(dp), INTENT(in) :: x, xi(:)
355 : INTEGER, INTENT(in) :: i_min
356 :
357 0 : IF (i_min >= 1 .AND. i_min <= SIZE(xi) - 1) THEN
358 0 : ip = iix(x, xi(i_min:)) + i_min - 1
359 : ELSE
360 0 : ip = iix(x, xi)
361 : END IF
362 0 : END FUNCTION
363 :
364 : !--------------------------------------------------------------------------------------------------!
365 :
366 : ! **************************************************************************************************
367 : !> \brief ...
368 : !> \param x ...
369 : !> \param n ...
370 : !> \param x1 ...
371 : !> \param xn ...
372 : !> \return ...
373 : ! **************************************************************************************************
374 0 : FUNCTION iixun(x, n, x1, xn)
375 : ! Returns index i of interval [x(i),x(i+1)] containing x in uniform mesh defined by
376 : ! x(i) = x1 + (i-1)/(n-1)*(xn-x1), i = 1 .. n,
377 : ! with intervals indexed by left-most points.
378 : ! N.B.: x outside [x1,xn] are indexed to nearest end.
379 : REAL(dp), INTENT(in) :: x
380 : INTEGER, INTENT(in) :: n
381 : REAL(dp), INTENT(in) :: x1, xn
382 : INTEGER :: iixun
383 :
384 : INTEGER :: i
385 :
386 : ! index i of interval [x(i),x(i+1)] containing x
387 : ! target value
388 : ! number of mesh points
389 : ! initial point of mesh
390 : ! final point of mesh
391 :
392 : ! compute index
393 0 : i = INT((x - x1)/(xn - x1)*(n - 1)) + 1
394 : ! reset if outside 1..n
395 0 : IF (i < 1) i = 1
396 0 : IF (i > n - 1) i = n - 1
397 0 : iixun = i
398 0 : END FUNCTION
399 :
400 : !--------------------------------------------------------------------------------------------------!
401 :
402 : ! **************************************************************************************************
403 : !> \brief ...
404 : !> \param x ...
405 : !> \param n ...
406 : !> \param x1 ...
407 : !> \param alpha ...
408 : !> \param beta ...
409 : !> \return ...
410 : ! **************************************************************************************************
411 0 : FUNCTION iixexp(x, n, x1, alpha, beta)
412 : ! Returns index i of interval [x(i),x(i+1)] containing x in exponential mesh defined by
413 : ! x(i) = x1 + alpha [ exp(beta(i-1)) - 1 ], i = 1 .. n,
414 : ! where alpha = (x(n) - x(1))/[ exp(beta(n-1)) - 1 ],
415 : ! beta = log(r)/(n-2), r = (x(n)-x(n-1))/(x(2)-x(1)) = ratio of last to first interval,
416 : ! and intervals indexed by left-most points.
417 : ! N.B.: x outside [x1,xn] are indexed to nearest end.
418 : REAL(dp), INTENT(in) :: x
419 : INTEGER, INTENT(in) :: n
420 : REAL(dp), INTENT(in) :: x1, alpha, beta
421 : INTEGER :: iixexp
422 :
423 : INTEGER :: i
424 :
425 : ! index i of interval [x(i),x(i+1)] containing x
426 : ! target value
427 : ! number of mesh points
428 : ! initial point of mesh
429 : ! mesh parameter:
430 : ! x(i) = x1 + alpha [ exp(beta(i-1)) - 1 ], i = 1 .. n,
431 : ! where alpha = (x(n) - x(1))/[ exp(beta(n-1)) - 1 ],
432 : ! beta = log(r)/(n-2), r = (x(n)-x(n-1))/(x(2)-x(1)) = ratio of last to first interval,
433 : ! mesh parameter
434 :
435 : ! compute index
436 0 : i = INT(LOG((x - x1)/alpha + 1)/beta) + 1
437 : ! reset if outside 1..n
438 0 : IF (i < 1) i = 1
439 0 : IF (i > n - 1) i = n - 1
440 0 : iixexp = i
441 0 : END FUNCTION
442 :
443 : !--------------------------------------------------------------------------------------------------!
444 :
445 : ! **************************************************************************************************
446 : !> \brief ...
447 : !> \param x ...
448 : !> \param c ...
449 : !> \return ...
450 : ! **************************************************************************************************
451 0 : FUNCTION poly3(x, c)
452 : ! returns value of polynomial \sum_{i=1}^4 c(i) (x-c(0))^(i-1)
453 : REAL(dp), INTENT(in) :: x, c(0:)
454 : REAL(dp) :: poly3
455 :
456 : REAL(dp) :: dx
457 :
458 : ! point at which to evaluate polynomial
459 : ! coefficients: poly = \sum_{i=1}^4 c(i) (x-c(0))^(i-1)
460 :
461 0 : dx = x - c(0)
462 0 : poly3 = c(1) + c(2)*dx + c(3)*dx**2 + c(4)*dx**3
463 0 : END FUNCTION
464 :
465 : !--------------------------------------------------------------------------------------------------!
466 :
467 : ! **************************************************************************************************
468 : !> \brief ...
469 : !> \param x ...
470 : !> \param c ...
471 : !> \return ...
472 : ! **************************************************************************************************
473 0 : FUNCTION dpoly3(x, c)
474 : ! returns 1st derivative of polynomial \sum_{i=1}^4 c(i) (x-c(0))^(i-1)
475 : REAL(dp), INTENT(in) :: x, c(0:)
476 : REAL(dp) :: dpoly3
477 :
478 : REAL(dp) :: dx
479 :
480 : ! point at which to evaluate polynomial
481 : ! coefficients: poly = \sum_{i=1}^4 c(i) (x-c(0))^(i-1)
482 :
483 0 : dx = x - c(0)
484 0 : dpoly3 = c(2) + 2*c(3)*dx + 3*c(4)*dx**2
485 0 : END FUNCTION
486 :
487 : !--------------------------------------------------------------------------------------------------!
488 :
489 : ! **************************************************************************************************
490 : !> \brief ...
491 : !> \param x ...
492 : !> \param c ...
493 : !> \return ...
494 : ! **************************************************************************************************
495 0 : FUNCTION d2poly3(x, c)
496 : ! returns 2nd derivative of polynomial \sum_{i=1}^4 c(i) (x-c(0))^(i-1)
497 : REAL(dp), INTENT(in) :: x, c(0:)
498 : REAL(dp) :: d2poly3
499 :
500 : REAL(dp) :: dx
501 :
502 : ! point at which to evaluate polynomial
503 : ! coefficients: poly = \sum_{i=1}^4 c(i) (x-c(0))^(i-1)
504 :
505 0 : dx = x - c(0)
506 0 : d2poly3 = 2*c(3) + 6*c(4)*dx
507 0 : END FUNCTION
508 :
509 : !--------------------------------------------------------------------------------------------------!
510 :
511 : ! **************************************************************************************************
512 : !> \brief ...
513 : !> \param msg ...
514 : ! **************************************************************************************************
515 0 : SUBROUTINE stop_error(msg)
516 : ! Aborts the program
517 : CHARACTER(LEN=*) :: msg
518 :
519 : ! Message to print on stdout
520 0 : CPABORT(msg)
521 0 : END SUBROUTINE
522 :
523 : END MODULE splines
|
