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