Line data Source code
1 : !!****m* ABINIT/m_splines
2 : !! NAME
3 : !! m_splines
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines for spline interpolation.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2010-2026 ABINIT group (YP, BAmadon)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_splines
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 :
28 : use m_fstrings, only : sjoin, itoa, ftoa
29 : !use m_time, only : timab
30 :
31 : implicit none
32 :
33 : public :: splfit
34 : public :: spline
35 : public :: spline_bicubic
36 : public :: spline_c
37 : public :: spline_complex
38 : public :: spline_integrate
39 : public :: splint
40 : public :: splint_complex
41 :
42 : !FIXME deprecated
43 : public :: intrpl
44 :
45 : ! *************************************************************************
46 :
47 : contains
48 : !!***
49 :
50 : !----------------------------------------------------------------------
51 :
52 : !!****f* m_splines/splfit
53 : !! NAME
54 : !! splfit
55 : !!
56 : !! FUNCTION
57 : !! Evaluate cubic spline fit to get function values on input set of ORDERED, UNFORMLY SPACED points.
58 : !! Optionally gives derivatives (first and second) at those points too.
59 : !! If point lies outside the range of arg, assign the extremal
60 : !! point values to these points, and zero derivative.
61 : !!
62 : !! INPUTS
63 : !! arg(numarg)=equally spaced arguments (spacing de) for data to which spline was fit.
64 : !! fun(numarg,2)=function values to which spline was fit and spline
65 : !! fit to second derivatives (from Numerical Recipes spline).
66 : !! ider= see above
67 : !! newarg(numnew)=new values of arguments at which function is desired.
68 : !! numarg=number of arguments at which spline was fit.
69 : !! numnew=number of arguments at which function values are desired.
70 : !!
71 : !! OUTPUT
72 : !! derfun(numnew)=(optional) values of first or second derivative of function.
73 : !! This is only computed for ider=1 or 2; otherwise derfun not used.
74 : !! newfun(numnew)=values of function at newarg(numnew).
75 : !! This is only computed for ider=0 or 1.
76 : !!
77 : !! NOTES
78 : !! if ider=0, compute only the function (contained in fun)
79 : !! if ider=1, compute the function (contained in fun) and its first derivative (in derfun)
80 : !! if ider=2, compute only the second derivative of the function (in derfun)
81 : !!
82 : !! SOURCE
83 :
84 17267854 : subroutine splfit(arg, derfun, fun, ider, newarg, newfun, numarg, numnew)
85 :
86 : integer, intent(in) :: ider, numarg, numnew
87 : real(dp), intent(in) :: arg(numarg), fun(numarg,2), newarg(numnew)
88 : real(dp), intent(out) :: derfun(numnew)
89 : real(dp), intent(inout) :: newfun(numnew)
90 :
91 : !Local variables---------------------------------------
92 : integer :: i,jspl
93 : real(dp) :: argmin,de,d,aa,bb,cc,dd,de2_dby_six,de_dby_six
94 : !real(dp) :: tsec(2)
95 :
96 : ! *************************************************************************
97 :
98 : ! argmin is smallest x value in spline fit; de is uniform spacing of spline argument
99 17267854 : argmin = arg(1)
100 17267854 : de = (arg(numarg) - argmin) / dble(numarg-1)
101 17267854 : de2_dby_six = de**2 / six
102 17267854 : de_dby_six = de / six
103 :
104 17267854 : if (de < tol12) then
105 0 : ABI_ERROR(sjoin('spacing should be strictly positive, while de is: ', ftoa(de)))
106 : endif
107 :
108 17267854 : jspl = -1
109 :
110 : ! Do one loop for no grads, other for grads
111 17267854 : select case (ider)
112 : case (0)
113 :
114 : ! Spline index loop for no grads:
115 2705894084 : do i=1,numnew
116 2705894084 : if (newarg(i) >= arg(numarg)) then
117 : ! function values are being requested outside range of data.',a1,'
118 : ! Function and slope will be set to values at upper end of data.
119 :
120 41689999 : newfun(i) = fun(numarg,1)
121 :
122 2653284947 : else if (newarg(i) <= arg(1)) then
123 60476326 : newfun(i) = fun(1,1)
124 :
125 : else
126 2592808621 : jspl = 1 + int((newarg(i) - argmin)/de)
127 2592808621 : d = newarg(i) - arg(jspl)
128 2592808621 : bb = d / de
129 2592808621 : aa = one - bb
130 2592808621 : cc = aa*(aa**2 -one) * de2_dby_six
131 2592808621 : dd = bb*(bb**2 -one) * de2_dby_six
132 2592808621 : newfun(i)= aa * fun(jspl,1) + bb*fun(jspl+1,1) + cc*fun(jspl,2) + dd*fun(jspl+1,2)
133 : end if
134 : enddo
135 :
136 : case (1)
137 :
138 : ! Spline index loop includes grads:
139 416714097 : do i=1,numnew
140 :
141 416714097 : if (newarg(i) >= arg(numarg)) then
142 28802920 : newfun(i) = fun(numarg,1)
143 28802920 : derfun(i) = zero
144 :
145 382718901 : else if (newarg(i) <= arg(1)) then
146 32356814 : newfun(i) = fun(1,1)
147 32356814 : derfun(i) = zero
148 :
149 : else
150 : ! cubic spline interpolation:
151 350362087 : jspl = 1 + int((newarg(i) - arg(1)) / de)
152 350362087 : d = newarg(i) - arg(jspl)
153 350362087 : bb = d / de
154 350362087 : aa = one - bb
155 350362087 : cc = aa*(aa**2 - one) * de2_dby_six
156 350362087 : dd = bb*(bb**2 - one) * de2_dby_six
157 350362087 : newfun(i) = aa*fun(jspl,1) + bb*fun(jspl+1,1) + cc*fun(jspl,2) + dd*fun(jspl+1,2)
158 : ! spline fit to first derivative:
159 : ! note correction of Numerical Recipes sign error
160 : derfun(i) = (fun(jspl+1,1)-fun(jspl,1)) / de + &
161 350362087 : (-(3.d0*aa**2 -one) * fun(jspl,2) + (3.d0*bb**2 -one) * fun(jspl+1,2)) * de_dby_six
162 :
163 : end if
164 : enddo
165 :
166 : case (2)
167 :
168 26829454 : do i=1,numnew
169 :
170 26829454 : if (newarg(i) >= arg(numarg)) then
171 0 : derfun(i) = zero
172 :
173 25673014 : else if (newarg(i) <= arg(1)) then
174 98 : derfun(i) = zero
175 :
176 : else
177 : ! cubic spline interpolation:
178 25672916 : jspl = 1 + int((newarg(i) - argmin) / de)
179 25672916 : d = newarg(i) - arg(jspl)
180 25672916 : bb = d / de
181 25672916 : aa = one - bb
182 : ! second derivative of spline (piecewise linear function)
183 25672916 : derfun(i) = aa*fun(jspl,2) + bb*fun(jspl+1,2)
184 :
185 : end if
186 : enddo
187 :
188 : case default
189 17267854 : ABI_ERROR(sjoin("Invalid ider:", itoa(ider)))
190 : end select
191 :
192 17267854 : end subroutine splfit
193 : !!***
194 :
195 : !----------------------------------------------------------------------
196 :
197 : !!****f* m_splines/spline
198 : !! NAME
199 : !! spline
200 : !!
201 : !! FUNCTION
202 : !! SPLINE (originally SPLINE_CUBIC_SET) computes the second derivatives
203 : !! of a cubic spline.
204 : !!
205 : !! INPUTS
206 : !! Input, integer N, the number of data points; N must be at least 2.
207 : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
208 : !! spline will actually be linear.
209 : !!
210 : !! Input, double precision T(N), the knot values, that is, the points where data
211 : !! is specified. The knot values should be distinct, and increasing.
212 : !!
213 : !! Input, double precision Y(N), the data values to be interpolated.
214 : !!
215 : !! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
216 : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
217 : !!
218 : !! OUTPUT
219 : !! Output, double precision YPP(N), the second derivatives of the cubic spline.
220 : !! Work space, double precision DIAG(N) - should be removed ...
221 : !!
222 : !! SOURCE
223 :
224 4453092 : subroutine spline( t, y, n, ybcbeg, ybcend, ypp )
225 :
226 : !*******************************************************************************
227 : !
228 : ! Discussion:
229 : !
230 : ! For data interpolation, the user must call SPLINE_CUBIC_SET to
231 : ! determine the second derivative data, passing in the data to be
232 : ! interpolated, and the desired boundary conditions.
233 : !
234 : ! The data to be interpolated, plus the SPLINE_CUBIC_SET output,
235 : ! defines the spline. The user may then call SPLINE_CUBIC_VAL to
236 : ! evaluate the spline at any point.
237 : !
238 : ! The cubic spline is a piecewise cubic polynomial. The intervals
239 : ! are determined by the "knots" or abscissas of the data to be
240 : ! interpolated. The cubic spline has continous first and second
241 : ! derivatives over the entire interval of interpolation.
242 : !
243 : ! For any point T in the interval T(IVAL), T(IVAL+1), the form of
244 : ! the spline is
245 : !
246 : ! SPL(T) = A(IVAL)
247 : ! + B(IVAL) * ( T - T(IVAL) )
248 : ! + C(IVAL) * ( T - T(IVAL) )**2
249 : ! + D(IVAL) * ( T - T(IVAL) )**3
250 : !
251 : ! If we assume that we know the values Y(*) and YPP(*), which represent
252 : ! the values and second derivatives of the spline at each knot, then
253 : ! the coefficients can be computed as:
254 : !
255 : ! A(IVAL) = Y(IVAL)
256 : ! B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
257 : ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
258 : ! C(IVAL) = YPP(IVAL) / 2
259 : ! D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
260 : !
261 : ! Since the first derivative of the spline is
262 : !
263 : ! SPL'(T) = B(IVAL)
264 : ! + 2 * C(IVAL) * ( T - T(IVAL) )
265 : ! + 3 * D(IVAL) * ( T - T(IVAL) )**2,
266 : !
267 : ! the requirement that the first derivative be continuous at interior
268 : ! knot I results in a total of N-2 equations, of the form:
269 : !
270 : ! B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
271 : ! + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
272 : !
273 : ! or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
274 : !
275 : ! ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
276 : ! - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
277 : ! + YPP(IVAL-1) * H(IVAL-1)
278 : ! + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
279 : ! =
280 : ! ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
281 : ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
282 : !
283 : ! or
284 : !
285 : ! YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
286 : ! + YPP(IVAL) * H(IVAL)
287 : ! =
288 : ! 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
289 : ! - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
290 : !
291 : ! Boundary conditions must be applied at the first and last knots.
292 : ! The resulting tridiagonal system can be solved for the YPP values.
293 : !
294 : ! Modified:
295 : !
296 : ! 07 February 1999
297 : ! 28 November 2004 XGonze : double precision
298 : ! make arguments similar to the Numeric Recipes routine
299 : ! also use algorithmics similar to the Numeric Recipes routine
300 : !
301 : ! Author:
302 : !
303 : ! John Burkardt
304 : ! (XGonze got it from http://www.psc.edu/~burkardt/src/spline/spline.html)
305 : !
306 : ! Parameters:
307 : !
308 : ! Input, integer N, the number of data points; N must be at least 2.
309 : ! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
310 : ! spline will actually be linear.
311 : !
312 : ! Input, double precision T(N), the knot values, that is, the points where data
313 : ! is specified. The knot values should be distinct, and increasing.
314 : !
315 : ! Input, double precision Y(N), the data values to be interpolated.
316 : !
317 : ! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
318 : ! conditions if IBCBEG or IBCEND is equal to 1 or 2.
319 : !
320 : ! Output, double precision YPP(N), the second derivatives of the cubic spline.
321 : !
322 : ! Work space, double precision DIAG(N) - should be removed ...
323 : !
324 : !
325 : ! XG041127 : In the initial implementation, one had the control on
326 : ! IBCBEG and IBCEND. Now, they are determined by the values
327 : ! of YBCBEG, YBCEND. Option 2 has been disabled.
328 : !
329 : ! Input, integer IBCBEG, left boundary condition flag:
330 : !
331 : ! 0: the spline should be a quadratic over the first interval;
332 : ! 1: the first derivative at the left endpoint should be YBCBEG;
333 : ! 2: the second derivative at the left endpoint should be YBCBEG.
334 : !
335 : ! Input, integer IBCEND, right boundary condition flag:
336 : !
337 : ! 0: the spline should be a quadratic over the last interval;
338 : ! 1: the first derivative at the right endpoint should be YBCEND;
339 : ! 2: the second derivative at the right endpoint should be YBCEND.
340 :
341 : integer, intent(in) :: n
342 : real(dp), intent(in) :: t(n)
343 : real(dp), intent(in) :: y(n)
344 : real(dp), intent(in) :: ybcbeg
345 : real(dp), intent(in) :: ybcend
346 :
347 : real(dp), intent(out) :: ypp(n)
348 :
349 : integer :: ibcbeg
350 : integer :: ibcend
351 : integer :: i,k
352 : real(dp) :: ratio,pinv
353 4453092 : real(dp), allocatable :: tmp(:)
354 : !
355 : ! Check.
356 : !
357 4453092 : if ( n <= 1 ) then
358 0 : write(std_out,* ) ' '
359 0 : write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!'
360 0 : write(std_out,* ) ' The number of knots must be at least 2.'
361 0 : write(std_out,* ) ' The input value of N = ', n
362 0 : ABI_ERROR("Fatal error")
363 : end if
364 :
365 13359276 : ABI_MALLOC(tmp,(n))
366 :
367 117126889 : do i = 1, n-1
368 117126889 : if ( t(i) >= t(i+1) ) then
369 0 : write(std_out,* ) ' '
370 0 : write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!'
371 0 : write(std_out,* ) ' The knots must be strictly increasing, but'
372 0 : write(std_out,* ) ' T(', i,') = ', t(i)
373 0 : write(std_out,* ) ' T(',i+1,') = ', t(i+1)
374 0 : ABI_ERROR("Fatal error")
375 : end if
376 : end do
377 : !
378 : ! XG041127
379 4453092 : ibcbeg=1 ; ibcend=1
380 4453092 : if(ybcbeg>1.0d+30)ibcbeg=0
381 4453092 : if(ybcend>1.0d+30)ibcend=0
382 : !
383 : ! Set the first and last equations.
384 : !
385 4453092 : if ( ibcbeg == 0 ) then
386 38916 : ypp(1) = 0.d0
387 38916 : tmp(1) = 0.d0
388 : else if ( ibcbeg == 1 ) then
389 4414176 : ypp(1) = -0.5d0
390 4414176 : tmp(1) = (3.d0/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg)
391 : end if
392 4453092 : if ( ibcend == 0 ) then
393 38916 : ypp(n) = 0.d0
394 38916 : tmp(n) = 0.d0
395 : else if ( ibcend == 1 ) then
396 4414176 : ypp(n) = 0.5d0
397 4414176 : tmp(n) = (3.d0/(t(n)-t(n-1)))*(ybcend-(y(n)-y(n-1))/(t(n)-t(n-1)))
398 : end if
399 :
400 : !
401 : ! Set the intermediate equations.
402 : !
403 112673797 : do i=2,n-1
404 108220705 : ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1))
405 108220705 : pinv = 1.0d0/(ratio*ypp(i-1) + 2.0d0)
406 108220705 : ypp(i) = (ratio-1.0d0)*pinv
407 : tmp(i)=(6.0d0*((y(i+1)-y(i))/(t(i+1)-t(i))-(y(i)-y(i-1)) &
408 108220705 : & /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv
409 112673797 : if (abs(tmp(i))<1.d5*tiny(0.d0)) tmp(i)=0.d0 !MT20050927
410 : enddo
411 :
412 : ! Solve the equations
413 4453092 : ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0d0)
414 117126889 : do k=n-1,1,-1
415 117126889 : ypp(k)=ypp(k)*ypp(k+1)+tmp(k)
416 : enddo
417 :
418 4453092 : ABI_FREE(tmp)
419 4453092 : end subroutine spline
420 : !!***
421 :
422 : !----------------------------------------------------------------------
423 :
424 : !!****f* m_splines/spline_bicubic
425 : !! NAME
426 : !! spline_bicubic
427 : !!
428 : !! FUNCTION
429 : !! Generates coefficients for bicubic spline interpolation.
430 : !!
431 : !! INPUTS
432 : !! n1 = length of first dimension
433 : !! n2 = length of second dimension
434 : !! x1 = positions on first dimension
435 : !! x2 = positions on second dimension
436 : !! y = function values on the (x1,x2) grid
437 : !! der1_x1 = first derivative of y wrt x1
438 : !! der1_x2 = first derivative of y wrt x2
439 : !! der2_x1x2 = second-order cross-derivative of y wrt x1x2
440 : !!
441 : !! OUTPUT
442 : !! spl_c = spline coefficients
443 : !!
444 : !! NOTES
445 : !! Adapted from Numerical Recipes and libbci.
446 : !!
447 : !! SOURCE
448 :
449 0 : subroutine spline_bicubic(n1,n2,x1,x2,y,der1_x1,der1_x2,der2_x1x2,spl_c)
450 :
451 : integer,intent(in) :: n1,n2
452 : real(dp),intent(in) :: x1(n1),x2(n2),y(n1,n2)
453 : real(dp),intent(in) :: der1_x1(n1,n2),der1_x2(n1,n2),der2_x1x2(n1,n2)
454 : real(dp),intent(out):: spl_c(4,4,n1,n2)
455 :
456 : integer :: i1,i2
457 : real(dp) :: dx1,dx2,wt(16,16),z(16)
458 :
459 : data wt /1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4, &
460 : & 8*0,3,0,-9,6,-2,0,6,-4,10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6, &
461 : & 2*0,6,-4,4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3, &
462 : & 2,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,0,1,-2,1,5*0, &
463 : & -3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,10*0,-3,3,2*0,2,-2,2*0, &
464 : & -1,1,6*0,3,-3,2*0,-2,2,5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2, &
465 : & -1,0,1,-2,1,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/
466 :
467 : ! Set coefficients for i1<n1 and i2<n2
468 0 : do i2 = 1,n2-1
469 0 : do i1 = 1,n1-1
470 0 : dx1 = x1(i1+1) - x1(i1)
471 0 : dx2 = x2(i2+1) - x2(i2)
472 0 : z(1) = y(i1,i2)
473 0 : z(2) = y(i1+1,i2)
474 0 : z(3) = y(i1+1,i2+1)
475 0 : z(4) = y(i1,i2+1)
476 0 : z(5) = der1_x1(i1,i2) * dx1
477 0 : z(6) = der1_x1(i1+1,i2) * dx1
478 0 : z(7) = der1_x1(i1+1,i2+1) * dx1
479 0 : z(8) = der1_x1(i1,i2+1) * dx1
480 0 : z(9) = der1_x2(i1,i2) * dx2
481 0 : z(10) = der1_x2(i1+1,i2) * dx2
482 0 : z(11) = der1_x2(i1+1,i2+1) * dx2
483 0 : z(12) = der1_x2(i1,i2+1) * dx2
484 0 : z(13) = der2_x1x2(i1,i2) * dx1 * dx2
485 0 : z(14) = der2_x1x2(i1+1,i2) * dx1 * dx2
486 0 : z(15) = der2_x1x2(i1+1,i2+1) * dx1 * dx2
487 0 : z(16) = der2_x1x2(i1,i2+1) * dx1 * dx2
488 0 : z = matmul(wt,z)
489 0 : spl_c(:,:,i1,i2) = reshape(z,(/4,4/),order=(/2,1/))
490 : end do
491 : end do
492 :
493 : ! Set coefficients for i1=n1 and i2=n2 (valid only at the border)
494 0 : spl_c(:,:,n1,:) = 0
495 0 : spl_c(:,:,:,n2) = 0
496 0 : spl_c(1,1,n1,:) = y(n1,:)
497 0 : spl_c(1,1,:,n2) = y(:,n2)
498 :
499 0 : end subroutine spline_bicubic
500 : !!***
501 :
502 : !----------------------------------------------------------------------
503 :
504 : !!****f* m_splines/spline_r
505 : !! NAME
506 : !! spline_r
507 : !!
508 : !! FUNCTION
509 : !! Computes the spline of a real function.
510 : !! If point lies outside the range of original grids, assign the extremal
511 : !! point values to either head or tail.
512 : !!
513 : !! INPUTS
514 : !! nomega_lo = number of point in the non regular grid (e.g. !logarithmic)
515 : !! nomega_li = number of point in the regular grid on which the spline is computed
516 : !! omega_lo = value of freq on the 1st grid
517 : !! omega_li = value of freq on the 2nd grid
518 : !! tospline_lo = function on the 1st grid
519 : !!
520 : !! OUTPUT
521 : !! splined_lo = spline (on the 2nd grid)
522 : !!
523 : !! SOURCE
524 :
525 0 : subroutine spline_r( nomega_lo, nomega_li, omega_lo, omega_li, splined_li, tospline_lo, extrapolate)
526 :
527 : !Arguments --------------------------------------------
528 : !scalars
529 : integer, intent(in) :: nomega_lo, nomega_li
530 : real(dp), intent(in) :: omega_lo(nomega_lo)
531 : real(dp), intent(in) :: omega_li(nomega_li)
532 : real(dp), intent(in) :: tospline_lo(nomega_lo)
533 : real(dp), intent(out) :: splined_li(nomega_li)
534 : logical, intent(in), optional :: extrapolate
535 :
536 : !Local variables---------------------------------------
537 : !scalars
538 : integer :: begin, end
539 : real(dp) :: ybcbeg, ybcend
540 0 : real(dp) :: ysplin2_lo(nomega_lo)
541 :
542 0 : ybcbeg=zero
543 0 : ybcend=zero
544 :
545 0 : begin = 1
546 0 : end = nomega_li
547 :
548 0 : call spline(omega_lo, tospline_lo, nomega_lo, ybcbeg, ybcend, ysplin2_lo)
549 0 : if (present(extrapolate)) then
550 0 : if (extrapolate) then
551 0 : do begin = 1, nomega_li
552 0 : if (omega_li(begin) >= omega_lo(1)) exit
553 : end do
554 0 : do end = nomega_li, 1, -1
555 0 : if (omega_li(end) <= omega_lo(nomega_lo)) exit
556 : end do
557 0 : ABI_CHECK(begin <= end, 'spline_c: omega_li not properly ordered')
558 : end if
559 : end if
560 0 : ABI_CHECK(begin <= end, 'spline_r: omega_li not properly ordered')
561 0 : call splint(nomega_lo, omega_lo, tospline_lo, ysplin2_lo, end-begin+1, omega_li(begin:end), splined_li(begin:end))
562 0 : if (present(extrapolate)) then
563 0 : if (extrapolate) then
564 0 : splined_li(1:begin-1) = tospline_lo(1)
565 0 : splined_li(end+1:nomega_li) = tospline_lo(nomega_lo)
566 : end if
567 : end if
568 :
569 0 : end subroutine spline_r
570 : !!***
571 :
572 : !----------------------------------------------------------------------
573 :
574 : !!****f* m_splines/spline_c
575 : !! NAME
576 : !! spline_c
577 : !!
578 : !! FUNCTION
579 : !! Computes the spline of a complex function.
580 : !!
581 : !! INPUTS
582 : !! nomega_lo = number of point in the non regular grid (e.g. !logarithmic)
583 : !! nomega_li = number of point in the regular grid on which the spline is computed
584 : !! omega_lo = value of freq on the 1st grid
585 : !! omega_li = value of freq on the 2nd grid
586 : !! tospline_lo = function on the 1st grid
587 : !!
588 : !! OUTPUT
589 : !! splined_lo = spline (on the 2nd grid)
590 : !!
591 : !! SOURCE
592 :
593 11570 : subroutine spline_c( nomega_lo, nomega_li, omega_lo, omega_li, splined_li, tospline_lo, extrapolate)
594 :
595 : !Arguments --------------------------------------------
596 : !scalars
597 : integer, intent(in) :: nomega_lo, nomega_li
598 : real(dp), intent(in) :: omega_lo(nomega_lo)
599 : real(dp), intent(in) :: omega_li(nomega_li)
600 : complex(dp), intent(in) :: tospline_lo(nomega_lo)
601 : complex(dp), intent(out) :: splined_li(nomega_li)
602 : logical, intent(in), optional :: extrapolate
603 :
604 : !Local variables---------------------------------------
605 : !scalars
606 : integer :: begin, end
607 : complex(dp) :: ybcbeg, ybcend
608 23140 : complex(dp) :: ysplin2_lo(nomega_lo)
609 :
610 11570 : ybcbeg=czero
611 11570 : ybcend=czero
612 :
613 11570 : begin = 1
614 11570 : end = nomega_li
615 :
616 11570 : call spline_complex(omega_lo, tospline_lo, nomega_lo, ybcbeg, ybcend, ysplin2_lo)
617 11570 : if (present(extrapolate)) then
618 11520 : if (extrapolate) then
619 23040 : do begin = 1, nomega_li
620 23040 : if (omega_li(begin) >= omega_lo(1)) exit
621 : end do
622 23040 : do end = nomega_li, 1, -1
623 23040 : if (omega_li(end) <= omega_lo(nomega_lo)) exit
624 : end do
625 11520 : ABI_CHECK(begin <= end, 'spline_c: omega_li not properly ordered')
626 : end if
627 : end if
628 11570 : call splint_complex( nomega_lo, omega_lo, tospline_lo,ysplin2_lo, end-begin+1, omega_li(begin:end), splined_li(begin:end))
629 11570 : if (present(extrapolate)) then
630 11520 : if (extrapolate) then
631 23040 : splined_li(1:begin-1) = tospline_lo(1)
632 23040 : splined_li(end+1:nomega_li) = tospline_lo(nomega_lo)
633 : end if
634 : end if
635 :
636 11570 : end subroutine spline_c
637 : !!***
638 :
639 : !----------------------------------------------------------------------
640 :
641 : !!****f* m_splines/spline_complex
642 : !! NAME
643 : !! spline_complex
644 : !!
645 : !! FUNCTION
646 : !! spline_complex interfaces the usual spline routine in a case of a
647 : !! complex function
648 : !!
649 : !! INPUTS
650 : !! Input, integer N, the number of data points; N must be at least 2.
651 : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
652 : !! spline will actually be linear.
653 : !!
654 : !! Input, double precision T(N), the knot values, that is, the points where data
655 : !! is specified. The knot values should be distinct, and increasing.
656 : !!
657 : !! Input, complex Y(N), the data values to be interpolated.
658 : !!
659 : !! Input, complex YBCBEG, YBCEND, the values to be used in the boundary
660 : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
661 : !!
662 : !! OUTPUT
663 : !! Output, complex YPP(N), the second derivatives of the cubic spline.
664 : !!
665 : !! SOURCE
666 :
667 29994 : subroutine spline_complex( t, y, n, ybcbeg, ybcend, ypp )
668 :
669 : integer, intent(in) :: n
670 : real(dp), intent(in) :: t(n)
671 : complex(dp), intent(in) :: y(n)
672 : complex(dp), intent(in) :: ybcbeg
673 : complex(dp), intent(in) :: ybcend
674 : complex(dp), intent(out) :: ypp(n)
675 :
676 29994 : real(dp), allocatable :: y_r(:)
677 : real(dp) :: ybcbeg_r
678 : real(dp) :: ybcend_r
679 29994 : real(dp), allocatable :: ypp_r(:)
680 29994 : real(dp), allocatable :: y_i(:)
681 : real(dp) :: ybcbeg_i
682 : real(dp) :: ybcend_i
683 29994 : real(dp), allocatable :: ypp_i(:)
684 :
685 89982 : ABI_MALLOC(y_r,(n))
686 59988 : ABI_MALLOC(ypp_r,(n))
687 59988 : ABI_MALLOC(y_i,(n))
688 59988 : ABI_MALLOC(ypp_i,(n))
689 8145078 : y_r=real(y)
690 8145078 : y_i=aimag(y) !vz_d
691 29994 : ybcbeg_r=real(ybcbeg)
692 29994 : ybcbeg_i=aimag(ybcbeg) !vz_d
693 29994 : ybcend_r=real(ybcend)
694 29994 : ybcend_i=aimag(ybcend) !vz_d
695 29994 : call spline( t, y_r, n, ybcbeg_r, ybcend_r, ypp_r )
696 29994 : call spline( t, y_i, n, ybcbeg_i, ybcend_i, ypp_i )
697 8115084 : ypp=cmplx(ypp_r,ypp_i)
698 29994 : ABI_FREE(y_r)
699 29994 : ABI_FREE(ypp_r)
700 29994 : ABI_FREE(y_i)
701 29994 : ABI_FREE(ypp_i)
702 :
703 29994 : end subroutine spline_complex
704 : !!***
705 :
706 : !----------------------------------------------------------------------
707 :
708 : !!****f* m_splines/splint
709 : !! NAME
710 : !! splint
711 : !!
712 : !! FUNCTION
713 : !! Compute spline interpolation. There is no hypothesis
714 : !! about the spacing of the input grid points.
715 : !!
716 : !! INPUTS
717 : !! nspline: number of grid points of input mesh
718 : !! xspline(nspline): input mesh
719 : !! yspline(nspline): function on input mesh
720 : !! ysplin2(nspline): second derivative of yspline on input mesh
721 : !! nfit: number of points of output mesh
722 : !! xfit(nfit): output mesh
723 : !!
724 : !! OUTPUT
725 : !! yfit(nfit): function on output mesh
726 : !! [ierr]=A non-zero value is used to signal that some points in xfit exceed xspline(nspline).
727 : !! The input value is incremented by the number of such points.
728 : !!
729 : !! SOURCE
730 :
731 50917729 : subroutine splint(nspline,xspline,yspline,ysplin2,nfit,xfit,yfit,ierr)
732 :
733 : integer, intent(in) :: nfit, nspline
734 : integer,optional,intent(out) :: ierr
735 : real(dp), intent(in) :: xspline(nspline)
736 : real(dp), intent(in) :: yspline(nspline)
737 : real(dp), intent(in) :: ysplin2(nspline)
738 : real(dp), intent(in) :: xfit(nfit)
739 : real(dp), intent(out) :: yfit(nfit)
740 :
741 : !local
742 : integer :: left,i,k,right,my_err
743 : real(dp) :: delarg,invdelarg,aa,bb
744 :
745 : !source
746 :
747 50917729 : my_err=0
748 :
749 50917729 : left = 1
750 12823376253 : do i=1, nfit
751 12772458524 : yfit(i)=0.d0 ! Initialize for the unlikely event that rmax exceed r(mesh)
752 : !
753 13051321540 : do k=left+1, nspline
754 13051321540 : if(xspline(k) >= xfit(i)) then
755 12757181822 : if(xspline(k-1) <= xfit(i)) then
756 : right = k
757 : left = k-1
758 : else
759 0 : if (k-1.eq.1 .and. i.eq.1) then
760 0 : ABI_ERROR('xfit(1) < xspline(1)')
761 : !my_err=my_err+1
762 : !exit
763 : else
764 0 : ABI_ERROR('xfit not properly ordered')
765 : end if
766 : end if
767 12757181822 : delarg= xspline(right) - xspline(left)
768 12757181822 : invdelarg= 1.0d0/delarg
769 12757181822 : aa= (xspline(right)-xfit(i))*invdelarg
770 12757181822 : bb= (xfit(i)-xspline(left))*invdelarg
771 :
772 : yfit(i) = aa*yspline(left) + bb*yspline(right) &
773 : & +( (aa*aa*aa-aa)*ysplin2(left) + &
774 12757181822 : & (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0d0
775 12757181822 : exit
776 : end if
777 : end do ! k
778 : !
779 12823376253 : if (k==nspline+1) my_err=my_err+1 ! xfit not found
780 : end do ! i
781 :
782 50917729 : if (PRESENT(ierr)) ierr=my_err
783 :
784 50917729 : end subroutine splint
785 : !!***
786 :
787 : !----------------------------------------------------------------------
788 :
789 : !!****f* m_splines/splint_complex
790 : !! NAME
791 : !! splint_complex
792 : !!
793 : !! FUNCTION
794 : !! Interface to the usual splint to compute *complex* spline interpolation. There is no hypothesis
795 : !! about the spacing of the input grid points.
796 : !!
797 : !! INPUTS
798 : !! nspline: number of grid points of input mesh
799 : !! xspline(nspline): input mesh
800 : !! yspline(nspline): complex function on input mesh
801 : !! ysplin2(nspline): second derivative of yspline on input mesh
802 : !! nfit: number of points of output mesh
803 : !! xfit(nfit): output mesh
804 : !!
805 : !! OUTPUT
806 : !! yfit(nfit): complex function on output mesh
807 : !!
808 : !! SOURCE
809 :
810 29994 : subroutine splint_complex (nspline,xspline,yspline,ysplin2,nfit,xfit,yfit)
811 :
812 : integer, intent(in) :: nfit, nspline
813 : real(dp), intent(in) :: xspline(nspline)
814 : complex(dp), intent(in) :: yspline(nspline)
815 : complex(dp), intent(in) :: ysplin2(nspline)
816 : real(dp), intent(in) :: xfit(nfit)
817 : complex(dp), intent(out) :: yfit(nfit)
818 :
819 29994 : real(dp), allocatable :: ysplin2_r(:)
820 29994 : real(dp), allocatable :: ysplin2_i(:)
821 29994 : real(dp), allocatable :: yspline_r(:)
822 29994 : real(dp), allocatable :: yspline_i(:)
823 29994 : real(dp), allocatable :: yfit_r(:)
824 29994 : real(dp), allocatable :: yfit_i(:)
825 :
826 89982 : ABI_MALLOC(yspline_r,(nspline))
827 59988 : ABI_MALLOC(yspline_i,(nspline))
828 59988 : ABI_MALLOC(ysplin2_r,(nspline))
829 59988 : ABI_MALLOC(ysplin2_i,(nspline))
830 89982 : ABI_MALLOC(yfit_r,(nfit))
831 59988 : ABI_MALLOC(yfit_i,(nfit))
832 :
833 : !local
834 :
835 : !source
836 8145078 : yspline_r=real(yspline)
837 8145078 : yspline_i=aimag(yspline) !vz_d
838 8145078 : ysplin2_r=real(ysplin2)
839 8145078 : ysplin2_i=aimag(ysplin2) !vz_d
840 29994 : call splint (nspline,xspline,yspline_r,ysplin2_r,nfit,xfit,yfit_r)
841 29994 : call splint (nspline,xspline,yspline_i,ysplin2_i,nfit,xfit,yfit_i)
842 6260473804 : yfit=cmplx(yfit_r,yfit_i)
843 29994 : ABI_FREE(yspline_r)
844 29994 : ABI_FREE(yspline_i)
845 29994 : ABI_FREE(ysplin2_r)
846 29994 : ABI_FREE(ysplin2_i)
847 29994 : ABI_FREE(yfit_r)
848 29994 : ABI_FREE(yfit_i)
849 :
850 29994 : end subroutine splint_complex
851 : !!***
852 :
853 : !!****f* m_splines/spline_integrate
854 : !! NAME
855 : !! spline_integrate
856 : !!
857 : !! FUNCTION
858 : !! Calculates an integral using cubic spline interpolation.
859 : !!
860 : !! INPUTS
861 : !! npts= number of grid points of input mesh
862 : !! dx= step of input mesh
863 : !! integrand= function on input mesh
864 : !!
865 : !! OUTPUT
866 : !! integral= integral of the input function
867 : !!
868 : !! SOURCE
869 :
870 0 : subroutine spline_integrate(integral,npts,dx,integrand)
871 :
872 : integer,intent(in) :: npts
873 : real(dp),intent(out) :: integral
874 : real(dp),intent(in) :: dx,integrand(npts)
875 :
876 : integer :: ix
877 0 : real(dp) :: ptmp,sf(npts),sf_der2(npts),sf_mesh(npts),utmp(npts)
878 :
879 : ! Prepare mesh
880 0 : forall (ix=1:npts) sf_mesh(ix) = (ix - 1) * dx
881 :
882 : ! Calculate second derivative of integrand (adapted from Numercial Recipes)
883 0 : sf_der2(1) = zero
884 0 : sf_der2(npts) = zero
885 0 : utmp(1) = zero
886 :
887 0 : do ix=2,npts-1
888 0 : ptmp = half * sf_der2(ix-1) + two
889 0 : sf_der2(ix) = (half - one) / ptmp
890 : utmp(ix) = (three * (integrand(ix+1) + integrand(ix-1) - &
891 0 : & two*integrand(ix)) / (dx**2) - half * utmp(ix-1)) / ptmp
892 : end do
893 0 : do ix=npts-1,1,-1
894 0 : sf_der2(ix) = sf_der2(ix) * sf_der2(ix+1) + utmp(ix)
895 : end do
896 :
897 : ! Actually calculate integral
898 0 : sf(:) = integrand(:) * dx
899 : integral = (sf(1) + sf(npts)) / 2.0_dp - &
900 : & (sf_der2(1) + sf_der2(npts)) / 24.0_dp + &
901 0 : & sum(sf(2:npts-1)) - sum(sf_der2(2:npts-1)) / 12.0_dp
902 :
903 0 : end subroutine spline_integrate
904 : !!***
905 :
906 : !!****f* m_splines/intrpl
907 : !! NAME
908 : !! intrpl
909 : !!
910 : !! FUNCTION
911 : !!
912 : !! DOUBLE PRECISION INTERPOLATION OF A SINGLE VALUED FUNCTION.
913 : !! THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION
914 : !! GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE
915 : !! AND FOR A GIVEN SET OF X VALUES(ABSCISSAE),THE VALUES OF
916 : !! A SINGLE VALUED FUNCTION Y=Y(X).
917 : !!
918 : !! THE SUBROUTINE ALSO CALCULATES FIRST DERIVATIVES DV(X) AND
919 : !! SECOND DERIVATIVE DV2(X)
920 : !
921 : !! THE INPUT PARAMETERS ARE;
922 : !!
923 : !! L=NUMBER OF DATA POINTS
924 : !! (MUST BE TWO OR GREATER)
925 : !! X=ARRAY OF DIMENSION L STORING THE X VALUES
926 : !! OF INPUT DATA POINTS (IN ASCENDING ORDER)
927 : !! Y=ARRAY OF DIMENSION L STORING THE Y VALUES OF INPUT DATA POINTS
928 : !! N=NUMBER OF POINTS AT WHICH INTERPOLATION OF THE Y-VALUES
929 : !! IS REQUIRED (MUST BE 1 OR GREATER)
930 : !! U=ARRAY OF DIMENSION N STORING THE X VALUES
931 : !! OF THE DESIRED POINTS
932 : !!
933 : !! THE OUTPUT PARAMETER IS V=ARRAY OF DIMENSION N WHERE THE
934 : !! INTERPOLATED Y VALUES ARE TO BE DISPLAYED
935 : !!
936 : !! INPUTS
937 : !! CUBIC SPLINE INTERPOLATION
938 : !!
939 : !! OUTPUT
940 : !!
941 : !! NOTES
942 : !! This routine is deprecated and will be replaced by the other routines of this module.
943 : !!
944 : !! SOURCE
945 :
946 0 : SUBROUTINE INTRPL(L,X,Y,N,U,V,dv,dv2,ideriv)
947 :
948 : IMPLICIT DOUBLE PRECISION (A-H,O-Z)
949 : IMPLICIT INTEGER(I-N)
950 : !
951 : PARAMETER (NQQ=12000)
952 :
953 : COMMON/QQ/ QQ(4,NQQ)
954 : DIMENSION X(L),Y(L),U(N),V(N),DV(NQQ),DV2(NQQ)
955 : EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3)
956 : REAL*8 M1,M2,M3,M4,M5
957 : EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5),&
958 : & (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
959 : !
960 : ! PRELIMINARY PROCESSING
961 :
962 0 : L0=L
963 0 : LM1=L0-1
964 0 : LM2=LM1-1
965 0 : LP1=L0+1
966 0 : N0=N
967 0 : IF(N0.GT.NQQ) THEN
968 0 : NQQV=NQQ
969 0 : write(std_out,2089) NQQV,N0
970 : ! CALL EXIT
971 : END IF
972 0 : IF(LM2.LT.0) GO TO 90
973 0 : IF(N0.LE.0) GO TO 91
974 0 : DO 11 I=2,L0
975 :
976 : ! IF(X(I-1)-X(I))11,95,96
977 0 : IF(X(I-1)-X(I).EQ.0.0D0) GO TO 95
978 0 : IF(X(I-1)-X(I).GT.0.0D0) GO TO 96
979 0 : 11 CONTINUE
980 : IPV=0
981 : !
982 : !*** MAIN LOOP
983 : FINT=0.0D0
984 0 : DO 80 K=1,N0
985 0 : UK=U(K)
986 : !
987 : !*** ROUTINE TO LOCATE THE DESIRED POINT
988 0 : IF(UK.GE.X(L0)) GO TO 26
989 0 : IF(UK.LT.X(1)) GO TO 25
990 0 : IMN=2
991 0 : IMX=L0
992 0 : 21 I=(IMN+IMX)/2
993 0 : IF(UK.GE.X(I)) GO TO 23
994 0 : IMX=I
995 0 : GO TO 24
996 0 : 23 IMN=I+1
997 0 : 24 IF(IMX.GT.IMN) GO TO 21
998 0 : I=IMX
999 0 : GO TO 30
1000 0 : 25 I=1
1001 0 : GO TO 30
1002 0 : 26 I=LP1
1003 0 : GO TO 30
1004 : !
1005 : !*** CHECK IF I=IPV
1006 0 : 30 IF(I.EQ.IPV) GO TO 70
1007 0 : IPV=I
1008 : !
1009 : !*** ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND TO
1010 : !*** ESTIMATE THEM IF NECESSARY
1011 0 : J=I
1012 0 : IF(J.EQ.1) J=2
1013 0 : IF(J.EQ.LP1) J=L0
1014 0 : X3=X(J-1)
1015 0 : Y3=Y(J-1)
1016 0 : X4=X(J)
1017 0 : Y4=Y(J)
1018 0 : A3=X4-X3
1019 0 : M3=(Y4-Y3)/A3
1020 0 : IF(LM2.EQ.0) GO TO 43
1021 0 : IF(J.EQ.2) GO TO 41
1022 0 : X2=X(J-2)
1023 0 : Y2=Y(J-2)
1024 0 : A2=X3-X2
1025 0 : M2=(Y3-Y2)/A2
1026 0 : IF(J.EQ.L0) GO TO 42
1027 0 : 41 X5=X(J+1)
1028 0 : Y5=Y(J+1)
1029 0 : A4=X5-X4
1030 0 : M4=(Y5-Y4)/A4
1031 0 : IF(J.EQ.2) M2=M3+M3-M4
1032 0 : GO TO 45
1033 0 : 42 M4=M3+M3-M2
1034 0 : GO TO 45
1035 : 43 M2=M3
1036 0 : 45 IF(J.LE.3) GO TO 46
1037 0 : A1=X2-X(J-3)
1038 0 : M1=(Y2-Y(J-3))/A1
1039 0 : GO TO 47
1040 0 : 46 M1=M2+M2-M3
1041 0 : 47 IF(J.GE.LM1) GO TO 48
1042 0 : A5=X(J+2)-X5
1043 0 : M5=(Y(J+2)-Y5)/A5
1044 0 : GO TO 50
1045 0 : 48 M5=M4+M4-M3
1046 : !
1047 : !*** NUMERICAL DIFFERENTIATION
1048 0 : 50 IF(I.EQ.LP1) GO TO 52
1049 0 : W2=ABS(M4-M3)
1050 0 : W3=ABS(M2-M1)
1051 0 : SW=W2+W3
1052 0 : IF(SW.NE.0.0) GO TO 51
1053 : W2=0.5D0
1054 : W3=0.5D0
1055 0 : SW=1.0D0
1056 0 : 51 T3=(W2*M2+W3*M3)/SW
1057 0 : IF(I.EQ.1) GO TO 54
1058 0 : 52 W3=ABS(M5-M4)
1059 0 : W4=ABS(M3-M2)
1060 0 : SW=W3+W4
1061 0 : IF(SW.NE.0.0) GO TO 53
1062 : W3=0.5D0
1063 : W4=0.5D0
1064 0 : SW=1.0D0
1065 0 : 53 T4=(W3*M3+W4*M4)/SW
1066 0 : IF(I.NE.LP1) GO TO 60
1067 : T3=T4
1068 0 : SA=A2+A3
1069 0 : T4=0.5D0*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA))
1070 : X3=X4
1071 : Y3=Y4
1072 0 : A3=A2
1073 0 : M3=M4
1074 0 : GO TO 60
1075 0 : 54 T4=T3
1076 0 : SA=A3+A4
1077 0 : T3=0.5D0*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA))
1078 0 : X3=X3-A4
1079 0 : Y3=Y3-M2*A4
1080 0 : A3=A4
1081 0 : M3=M2
1082 : !
1083 : !*** COMPUTATION OF THE POLYNOMIAL
1084 0 : 60 Q2=(2.0D0*(M3-T3)+M3-T4)/A3
1085 0 : Q3=(-M3-M3+T3+T4)/(A3*A3)
1086 0 : 70 DX=UK-P0
1087 0 : V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3))
1088 :
1089 0 : IF(IDERIV.EQ.0) GO TO 80
1090 0 : DV(K)=Q1+DX*(2.0D0*Q2+DX*3.0D0*Q3)
1091 0 : DV2(k)=6.0D0*Q3*DX+2.d0*Q2
1092 0 : QQ(1,K)=Q0
1093 0 : QQ(2,K)=Q1
1094 0 : QQ(3,K)=Q2
1095 0 : QQ(4,K)=Q3
1096 0 : 80 CONTINUE
1097 0 : RETURN
1098 : !
1099 : !*** ERROR EXIT
1100 0 : 90 write(std_out,2090)
1101 0 : GO TO 99
1102 0 : 91 write(std_out,2091)
1103 0 : GO TO 99
1104 0 : 95 write(std_out,2095)
1105 0 : GO TO 97
1106 0 : 96 write(std_out,2096)
1107 0 : 97 write(std_out,2097)I,X(I)
1108 0 : 99 write(std_out,2099) L0,N0
1109 0 : RETURN
1110 : !
1111 : !*** FORMAT STATEMENTS
1112 : 2089 FORMAT( 'WARNING ERROR IN INTRPL. MAX ALLOWED VALUE OF N0 IS',&
1113 : & I3,' HERE N0 IS',I3)
1114 : 2090 FORMAT(1X/' N = 1 OR LESS.'/)
1115 : 2091 FORMAT(1X/' N = 0 OR LESS.'/)
1116 : 2095 FORMAT(1X/' IDENTICAL X VALUES.'/)
1117 : 2096 FORMAT(1X/' X VALUES OUT OF SEQUENCE.'/)
1118 : 2097 FORMAT(4X,'I =',I7,10X,6X,'X(I) =',E12.3)
1119 : 2099 FORMAT(4X,'L =',I7,10X,3X,'N =',I7/ &
1120 : & ' ERROR DETECTED IN ROUTINE INTRPL')
1121 : !
1122 : END subroutine intrpl
1123 : !!***
1124 :
1125 : !!****f* m_splines/spline2
1126 : !! NAME
1127 : !! spline2
1128 : !!
1129 : !! FUNCTION
1130 : !! SPLINE2 computes the first derivatives of a cubic spline.
1131 : !!
1132 : !! INPUTS
1133 : !! Input, integer N, the number of data points; N must be at least 2.
1134 : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
1135 : !! spline will actually be linear.
1136 : !!
1137 : !! Input, double precision X(N), the knot values, that is, the points where data
1138 : !! is specified. The knot values should be distinct, and increasing.
1139 : !!
1140 : !! Input, double precision Y(N), the data values to be interpolated.
1141 : !!
1142 : !! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
1143 : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
1144 : !!
1145 : !! Input, integer IBCBEG, IBCEND, the type of boundary conditions at the first / last point
1146 : !! If 1, the value of the first derivative is constrained to YBCBEG / YBCEND
1147 : !! If 2, the value of the second derivative is constrained to YBCBEG / YBCEND
1148 : !! If 3, the third derivative is continuous at the second / second to last point (not-a-knot)
1149 : !!
1150 : !! OUTPUT
1151 : !! Output, double precision YP(N), the first derivatives of the cubic spline.
1152 : !! On [x(i),x(i+1)], the spline y(i)+yp(i)*(x-x(i))+c*(x-x(i))**2+d*(x-x(i))**3 can then
1153 : !! be reconstructed with c = (3*s-2*yp(i)-yp(i+1))/(x(i+1)-x(i))
1154 : !! d = (yp(i)+yp(i+1)-2*s)/(x(i+1)-x(i))**2
1155 : !! s = (y(i+1)-y(i))/(x(i+1)-x(i))
1156 : !!
1157 : !! Solving the tridiagonal system on the first derivatives rather than the
1158 : !! second derivatives as in the regular spline subroutine allows for the
1159 : !! fantastic "not-a-knot" boundary condition.
1160 : !!
1161 : !! SOURCE
1162 :
1163 0 : subroutine spline2(x,y,n,yp,ybcbeg,ybcend,ibcbeg,ibcend)
1164 :
1165 : integer,intent(in) :: n,ibcbeg,ibcend
1166 : real(dp),intent(in) :: ybcbeg,ybcend
1167 : real(dp),intent(in) :: x(n),y(n)
1168 : real(dp),intent(inout) :: yp(n)
1169 :
1170 : integer :: i
1171 : real(dp) :: pinv,ratio
1172 0 : real(dp),allocatable :: tmp(:)
1173 :
1174 0 : ABI_MALLOC(tmp,(n))
1175 :
1176 0 : if (ibcbeg==1) then
1177 0 : yp(1)=zero ; tmp(1)=ybcbeg
1178 0 : else if (ibcbeg==2) then
1179 0 : yp(1)=-half ; tmp(1)=(six*((y(2)-y(1))/(x(2)-x(1)))-ybcbeg*(x(2)-x(1)))/four
1180 0 : else if (ibcbeg==3) then
1181 0 : yp(1)=(x(1)-x(3))/(x(3)-x(2))
1182 : tmp(1)=((y(2)-y(1))*(x(3)-x(2))*(two*x(3)+x(2)-three*x(1))/(x(2)-x(1))+(y(3)-y(2))* &
1183 0 : & (x(2)-x(1))**2/(x(3)-x(2)))/((x(3)-x(2))*(x(3)-x(1)))
1184 : else
1185 0 : ABI_BUG("Option not recognized for ibcbeg !")
1186 : end if
1187 :
1188 0 : if (ibcend==1) then
1189 0 : yp(n)=zero ; tmp(n)=ybcend
1190 0 : else if (ibcend==2) then
1191 0 : yp(n)=half ; tmp(n)=(six*(y(n)-y(n-1))/(x(n)-x(n-1))+ybcend*(x(n)-x(n-1)))/four
1192 0 : else if (ibcend==3) then
1193 0 : yp(n)=(x(n)-x(n-2))/(x(n-1)-x(n-2))
1194 : tmp(n)=((y(n-1)-y(n-2))*(x(n)-x(n-1))**2/(x(n-1)-x(n-2))+ &
1195 : & (y(n)-y(n-1))*(x(n-1)-x(n-2))*(three*x(n)-x(n-1)-two*x(n-2)) &
1196 0 : & /(x(n)-x(n-1)))/((x(n-1)-x(n-2))*(x(n)-x(n-2)))
1197 : else
1198 0 : ABI_BUG("Option not recognized for ibcend !")
1199 : end if
1200 :
1201 0 : do i=2,n-1
1202 0 : ratio = (x(i+1)-x(i))/(x(i+1)-x(i-1))
1203 0 : pinv = one/(two+ratio*yp(i-1))
1204 0 : yp(i) = (ratio-one)*pinv
1205 : tmp(i) = (three*((y(i)-y(i-1))*(x(i+1)-x(i))/(x(i)-x(i-1))+ &
1206 : & (y(i+1)-y(i))*(x(i)-x(i-1))/(x(i+1)-x(i)))/(x(i+1)-x(i-1))-&
1207 0 : & ratio*tmp(i-1))*pinv
1208 0 : if (abs(tmp(i))<1.d5*tiny(zero)) tmp(i) = zero !MT20050927
1209 : end do
1210 :
1211 0 : yp(n) = (tmp(n)-yp(n)*tmp(n-1))/(yp(n)*yp(n-1)+one)
1212 0 : do i=n-1,1,-1
1213 0 : yp(i)=yp(i)*yp(i+1)+tmp(i)
1214 : end do
1215 :
1216 0 : ABI_FREE(tmp)
1217 :
1218 0 : end subroutine spline2
1219 : !!***
1220 :
1221 : !!****f* m_splines/spline2_complex
1222 : !! NAME
1223 : !! spline2_complex
1224 : !!
1225 : !! FUNCTION
1226 : !! SPLINE2_COMPLEX computes the first derivatives of a cubic spline for a complex function.
1227 : !!
1228 : !! INPUTS
1229 : !! Input, integer N, the number of data points; N must be at least 2.
1230 : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
1231 : !! spline will actually be linear.
1232 : !!
1233 : !! Input, double precision X(N), the knot values, that is, the points where data
1234 : !! is specified. The knot values should be distinct, and increasing.
1235 : !!
1236 : !! Input, double precision Y(N), the data values to be interpolated.
1237 : !!
1238 : !! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
1239 : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
1240 : !!
1241 : !! Input, integer IBCBEG, IBCEND, the type of boundary conditions at the first / last point
1242 : !! If 1, the value of the first derivative is constrained to YBCBEG / YBCEND
1243 : !! If 2, the value of the second derivative is constrained to YBCBEG / YBCEND
1244 : !! If 3, the third derivative is continuous at the second / second to last point (not-a-knot)
1245 : !!
1246 : !! OUTPUT
1247 : !! Output, double precision YP(N), the first derivatives of the cubic spline.
1248 : !! On [x(i),x(i+1)], the spline y(i)+yp(i)*(x-x(i))+c*(x-x(i))**2+d*(x-x(i))**3 can then
1249 : !! be reconstructed with c = (3*s-2*yp(i)-yp(i+1))/(x(i+1)-x(i))
1250 : !! d = (yp(i)+yp(i+1)-2*s)/(x(i+1)-x(i))**2
1251 : !! s = (y(i+1)-y(i))/(x(i+1)-x(i))
1252 : !!
1253 : !! Solving the tridiagonal system on the first derivatives rather than the
1254 : !! second derivatives as in the regular spline subroutine allows for the
1255 : !! fantastic "not-a-knot" boundary condition.
1256 : !!
1257 : !! SOURCE
1258 :
1259 0 : subroutine spline2_complex(x,y,n,yp,ybcbeg,ybcend,ibcbeg,ibcend)
1260 :
1261 : integer,intent(in) :: n,ibcbeg,ibcend
1262 : complex(dp),intent(in) :: ybcbeg,ybcend
1263 : real(dp),intent(in) :: x(n)
1264 : complex(dp),intent(in) :: y(n)
1265 : complex(dp),intent(inout) :: yp(n)
1266 :
1267 0 : real(dp),allocatable :: y_tmp(:),yp_i(:),yp_r(:)
1268 :
1269 0 : ABI_MALLOC(y_tmp,(n))
1270 0 : ABI_MALLOC(yp_i,(n))
1271 0 : ABI_MALLOC(yp_r,(n))
1272 :
1273 0 : y_tmp(:) = dble(y(:))
1274 0 : call spline2(x(:),y_tmp(:),n,yp_r(:),dble(ybcbeg),dble(ybcend),ibcbeg,ibcend)
1275 0 : y_tmp(:) = aimag(y(:))
1276 0 : call spline2(x(:),y_tmp(:),n,yp_i(:),aimag(ybcbeg),aimag(ybcend),ibcbeg,ibcend)
1277 0 : yp(:) = cmplx(yp_r(:),yp_i(:),kind=dp)
1278 :
1279 0 : ABI_FREE(y_tmp)
1280 0 : ABI_FREE(yp_i)
1281 0 : ABI_FREE(yp_r)
1282 :
1283 0 : end subroutine spline2_complex
1284 : !!***
1285 :
1286 : end module m_splines
1287 : !!***
|
