TABLE OF CONTENTS
- ABINIT/m_splines
- m_splines/intrpl
- m_splines/splfit
- m_splines/spline
- m_splines/spline_bicubic
- m_splines/spline_c
- m_splines/spline_complex
- m_splines/spline_integrate
- m_splines/splint
- m_splines/splint_complex
ABINIT/m_splines [ Modules ]
NAME
m_splines
FUNCTION
This module contains routines for spline interpolation.
COPYRIGHT
Copyright (C) 2010-2024 ABINIT group (YP, BAmadon) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
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
m_splines/intrpl [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
intrpl
FUNCTION
DOUBLE PRECISION INTERPOLATION OF A SINGLE VALUED FUNCTION. THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE AND FOR A GIVEN SET OF X VALUES(ABSCISSAE),THE VALUES OF A SINGLE VALUED FUNCTION Y=Y(X). THE SUBROUTINE ALSO CALCULATES FIRST DERIVATIVES DV(X) AND SECOND DERIVATIVE DV2(X) THE INPUT PARAMETERS ARE; L=NUMBER OF DATA POINTS (MUST BE TWO OR GREATER) X=ARRAY OF DIMENSION L STORING THE X VALUES OF INPUT DATA POINTS (IN ASCENDING ORDER) Y=ARRAY OF DIMENSION L STORING THE Y VALUES OF INPUT DATA POINTS N=NUMBER OF POINTS AT WHICH INTERPOLATION OF THE Y-VALUES IS REQUIRED (MUST BE 1 OR GREATER) U=ARRAY OF DIMENSION N STORING THE X VALUES OF THE DESIRED POINTS THE OUTPUT PARAMETER IS V=ARRAY OF DIMENSION N WHERE THE INTERPOLATED Y VALUES ARE TO BE DISPLAYED
INPUTS
CUBIC SPLINE INTERPOLATION
OUTPUT
NOTES
This routine is deprecated and will be replaced by the other routines of this module.
SOURCE
856 SUBROUTINE INTRPL(L,X,Y,N,U,V,dv,dv2,ideriv) 857 858 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 859 IMPLICIT INTEGER(I-N) 860 ! 861 PARAMETER (NQQ=12000) 862 863 COMMON/QQ/ QQ(4,NQQ) 864 DIMENSION X(L),Y(L),U(N),V(N),DV(NQQ),DV2(NQQ) 865 EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3) 866 REAL*8 M1,M2,M3,M4,M5 867 EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5),& 868 & (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3) 869 ! 870 ! PRELIMINARY PROCESSING 871 872 L0=L 873 LM1=L0-1 874 LM2=LM1-1 875 LP1=L0+1 876 N0=N 877 IF(N0.GT.NQQ) THEN 878 NQQV=NQQ 879 write(std_out,2089) NQQV,N0 880 ! CALL EXIT 881 END IF 882 IF(LM2.LT.0) GO TO 90 883 IF(N0.LE.0) GO TO 91 884 DO 11 I=2,L0 885 886 ! IF(X(I-1)-X(I))11,95,96 887 IF(X(I-1)-X(I).EQ.0.0D0) GO TO 95 888 IF(X(I-1)-X(I).GT.0.0D0) GO TO 96 889 11 CONTINUE 890 IPV=0 891 ! 892 !*** MAIN LOOP 893 FINT=0.0D0 894 DO 80 K=1,N0 895 UK=U(K) 896 ! 897 !*** ROUTINE TO LOCATE THE DESIRED POINT 898 IF(UK.GE.X(L0)) GO TO 26 899 IF(UK.LT.X(1)) GO TO 25 900 IMN=2 901 IMX=L0 902 21 I=(IMN+IMX)/2 903 IF(UK.GE.X(I)) GO TO 23 904 IMX=I 905 GO TO 24 906 23 IMN=I+1 907 24 IF(IMX.GT.IMN) GO TO 21 908 I=IMX 909 GO TO 30 910 25 I=1 911 GO TO 30 912 26 I=LP1 913 GO TO 30 914 ! 915 !*** CHECK IF I=IPV 916 30 IF(I.EQ.IPV) GO TO 70 917 IPV=I 918 ! 919 !*** ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND TO 920 !*** ESTIMATE THEM IF NECESSARY 921 J=I 922 IF(J.EQ.1) J=2 923 IF(J.EQ.LP1) J=L0 924 X3=X(J-1) 925 Y3=Y(J-1) 926 X4=X(J) 927 Y4=Y(J) 928 A3=X4-X3 929 M3=(Y4-Y3)/A3 930 IF(LM2.EQ.0) GO TO 43 931 IF(J.EQ.2) GO TO 41 932 X2=X(J-2) 933 Y2=Y(J-2) 934 A2=X3-X2 935 M2=(Y3-Y2)/A2 936 IF(J.EQ.L0) GO TO 42 937 41 X5=X(J+1) 938 Y5=Y(J+1) 939 A4=X5-X4 940 M4=(Y5-Y4)/A4 941 IF(J.EQ.2) M2=M3+M3-M4 942 GO TO 45 943 42 M4=M3+M3-M2 944 GO TO 45 945 43 M2=M3 946 45 IF(J.LE.3) GO TO 46 947 A1=X2-X(J-3) 948 M1=(Y2-Y(J-3))/A1 949 GO TO 47 950 46 M1=M2+M2-M3 951 47 IF(J.GE.LM1) GO TO 48 952 A5=X(J+2)-X5 953 M5=(Y(J+2)-Y5)/A5 954 GO TO 50 955 48 M5=M4+M4-M3 956 ! 957 !*** NUMERICAL DIFFERENTIATION 958 50 IF(I.EQ.LP1) GO TO 52 959 W2=ABS(M4-M3) 960 W3=ABS(M2-M1) 961 SW=W2+W3 962 IF(SW.NE.0.0) GO TO 51 963 W2=0.5D0 964 W3=0.5D0 965 SW=1.0D0 966 51 T3=(W2*M2+W3*M3)/SW 967 IF(I.EQ.1) GO TO 54 968 52 W3=ABS(M5-M4) 969 W4=ABS(M3-M2) 970 SW=W3+W4 971 IF(SW.NE.0.0) GO TO 53 972 W3=0.5D0 973 W4=0.5D0 974 SW=1.0D0 975 53 T4=(W3*M3+W4*M4)/SW 976 IF(I.NE.LP1) GO TO 60 977 T3=T4 978 SA=A2+A3 979 T4=0.5D0*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA)) 980 X3=X4 981 Y3=Y4 982 A3=A2 983 M3=M4 984 GO TO 60 985 54 T4=T3 986 SA=A3+A4 987 T3=0.5D0*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA)) 988 X3=X3-A4 989 Y3=Y3-M2*A4 990 A3=A4 991 M3=M2 992 ! 993 !*** COMPUTATION OF THE POLYNOMIAL 994 60 Q2=(2.0D0*(M3-T3)+M3-T4)/A3 995 Q3=(-M3-M3+T3+T4)/(A3*A3) 996 70 DX=UK-P0 997 V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3)) 998 999 IF(IDERIV.EQ.0) GO TO 80 1000 DV(K)=Q1+DX*(2.0D0*Q2+DX*3.0D0*Q3) 1001 DV2(k)=6.0D0*Q3*DX+2.d0*Q2 1002 QQ(1,K)=Q0 1003 QQ(2,K)=Q1 1004 QQ(3,K)=Q2 1005 QQ(4,K)=Q3 1006 80 CONTINUE 1007 RETURN 1008 ! 1009 !*** ERROR EXIT 1010 90 write(std_out,2090) 1011 GO TO 99 1012 91 write(std_out,2091) 1013 GO TO 99 1014 95 write(std_out,2095) 1015 GO TO 97 1016 96 write(std_out,2096) 1017 97 write(std_out,2097)I,X(I) 1018 99 write(std_out,2099) L0,N0 1019 RETURN 1020 ! 1021 !*** FORMAT STATEMENTS 1022 2089 FORMAT( 'WARNING ERROR IN INTRPL. MAX ALLOWED VALUE OF N0 IS',& 1023 & I3,' HERE N0 IS',I3) 1024 2090 FORMAT(1X/' N = 1 OR LESS.'/) 1025 2091 FORMAT(1X/' N = 0 OR LESS.'/) 1026 2095 FORMAT(1X/' IDENTICAL X VALUES.'/) 1027 2096 FORMAT(1X/' X VALUES OUT OF SEQUENCE.'/) 1028 2097 FORMAT(4X,'I =',I7,10X,6X,'X(I) =',E12.3) 1029 2099 FORMAT(4X,'L =',I7,10X,3X,'N =',I7/ & 1030 & ' ERROR DETECTED IN ROUTINE INTRPL') 1031 ! 1032 END subroutine intrpl
m_splines/splfit [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
splfit
FUNCTION
Evaluate cubic spline fit to get function values on input set of ORDERED, UNFORMLY SPACED points. Optionally gives derivatives (first and second) at those points too. If point lies outside the range of arg, assign the extremal point values to these points, and zero derivative.
INPUTS
arg(numarg)=equally spaced arguments (spacing de) for data to which spline was fit. fun(numarg,2)=function values to which spline was fit and spline fit to second derivatives (from Numerical Recipes spline). ider= see above newarg(numnew)=new values of arguments at which function is desired. numarg=number of arguments at which spline was fit. numnew=number of arguments at which function values are desired.
OUTPUT
derfun(numnew)=(optional) values of first or second derivative of function. This is only computed for ider=1 or 2; otherwise derfun not used. newfun(numnew)=values of function at newarg(numnew). This is only computed for ider=0 or 1.
NOTES
if ider=0, compute only the function (contained in fun) if ider=1, compute the function (contained in fun) and its first derivative (in derfun) if ider=2, compute only the second derivative of the function (in derfun)
SOURCE
84 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 argmin = arg(1) 100 de = (arg(numarg) - argmin) / dble(numarg-1) 101 de2_dby_six = de**2 / six 102 de_dby_six = de / six 103 104 if (de < tol12) then 105 ABI_ERROR(sjoin('spacing should be strictly positive, while de is: ', ftoa(de))) 106 endif 107 108 jspl = -1 109 110 ! Do one loop for no grads, other for grads 111 select case (ider) 112 case (0) 113 114 ! Spline index loop for no grads: 115 do i=1,numnew 116 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 newfun(i) = fun(numarg,1) 121 122 else if (newarg(i) <= arg(1)) then 123 newfun(i) = fun(1,1) 124 125 else 126 jspl = 1 + int((newarg(i) - argmin)/de) 127 d = newarg(i) - arg(jspl) 128 bb = d / de 129 aa = one - bb 130 cc = aa*(aa**2 -one) * de2_dby_six 131 dd = bb*(bb**2 -one) * de2_dby_six 132 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 do i=1,numnew 140 141 if (newarg(i) >= arg(numarg)) then 142 newfun(i) = fun(numarg,1) 143 derfun(i) = zero 144 145 else if (newarg(i) <= arg(1)) then 146 newfun(i) = fun(1,1) 147 derfun(i) = zero 148 149 else 150 ! cubic spline interpolation: 151 jspl = 1 + int((newarg(i) - arg(1)) / de) 152 d = newarg(i) - arg(jspl) 153 bb = d / de 154 aa = one - bb 155 cc = aa*(aa**2 - one) * de2_dby_six 156 dd = bb*(bb**2 - one) * de2_dby_six 157 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 (-(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 do i=1,numnew 169 170 if (newarg(i) >= arg(numarg)) then 171 derfun(i) = zero 172 173 else if (newarg(i) <= arg(1)) then 174 derfun(i) = zero 175 176 else 177 ! cubic spline interpolation: 178 jspl = 1 + int((newarg(i) - argmin) / de) 179 d = newarg(i) - arg(jspl) 180 bb = d / de 181 aa = one - bb 182 ! second derivative of spline (piecewise linear function) 183 derfun(i) = aa*fun(jspl,2) + bb*fun(jspl+1,2) 184 185 end if 186 enddo 187 188 case default 189 ABI_ERROR(sjoin("Invalid ider:", itoa(ider))) 190 end select 191 192 end subroutine splfit
m_splines/spline [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
spline
FUNCTION
SPLINE (originally SPLINE_CUBIC_SET) computes the second derivatives of a cubic spline.
INPUTS
Input, integer N, the number of data points; N must be at least 2.
In the special case where N = 2 and IBCBEG = IBCEND = 0, the
spline will actually be linear.
Input, double precision T(N), the knot values, that is, the points where data
is specified. The knot values should be distinct, and increasing.
Input, double precision Y(N), the data values to be interpolated.
Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
conditions if IBCBEG or IBCEND is equal to 1 or 2.
OUTPUT
Output, double precision YPP(N), the second derivatives of the cubic spline.
Work space, double precision DIAG(N) - should be removed ...
SOURCE
224 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 real(dp), allocatable :: tmp(:) 354 ! 355 ! Check. 356 ! 357 if ( n <= 1 ) then 358 write(std_out,* ) ' ' 359 write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!' 360 write(std_out,* ) ' The number of knots must be at least 2.' 361 write(std_out,* ) ' The input value of N = ', n 362 ABI_ERROR("Fatal error") 363 end if 364 365 ABI_MALLOC(tmp,(n)) 366 367 do i = 1, n-1 368 if ( t(i) >= t(i+1) ) then 369 write(std_out,* ) ' ' 370 write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!' 371 write(std_out,* ) ' The knots must be strictly increasing, but' 372 write(std_out,* ) ' T(', i,') = ', t(i) 373 write(std_out,* ) ' T(',i+1,') = ', t(i+1) 374 ABI_ERROR("Fatal error") 375 end if 376 end do 377 ! 378 ! XG041127 379 ibcbeg=1 ; ibcend=1 380 if(ybcbeg>1.0d+30)ibcbeg=0 381 if(ybcend>1.0d+30)ibcend=0 382 ! 383 ! Set the first and last equations. 384 ! 385 if ( ibcbeg == 0 ) then 386 ypp(1) = 0.d0 387 tmp(1) = 0.d0 388 else if ( ibcbeg == 1 ) then 389 ypp(1) = -0.5d0 390 tmp(1) = (3.d0/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg) 391 end if 392 if ( ibcend == 0 ) then 393 ypp(n) = 0.d0 394 tmp(n) = 0.d0 395 else if ( ibcend == 1 ) then 396 ypp(n) = 0.5d0 397 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 do i=2,n-1 404 ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1)) 405 pinv = 1.0d0/(ratio*ypp(i-1) + 2.0d0) 406 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 & /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv 409 if (abs(tmp(i))<1.d5*tiny(0.d0)) tmp(i)=0.d0 !MT20050927 410 enddo 411 412 ! Solve the equations 413 ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0d0) 414 do k=n-1,1,-1 415 ypp(k)=ypp(k)*ypp(k+1)+tmp(k) 416 enddo 417 418 ABI_FREE(tmp) 419 end subroutine spline
m_splines/spline_bicubic [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
spline_bicubic
FUNCTION
Generates coefficients for bicubic spline interpolation.
INPUTS
n1 = length of first dimension n2 = length of second dimension x1 = positions on first dimension x2 = positions on second dimension y = function values on the (x1,x2) grid der1_x1 = first derivative of y wrt x1 der1_x2 = first derivative of y wrt x2 der2_x1x2 = second-order cross-derivative of y wrt x1x2
OUTPUT
spl_c = spline coefficients
NOTES
Adapted from Numerical Recipes and libbci.
SOURCE
449 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 do i2 = 1,n2-1 469 do i1 = 1,n1-1 470 dx1 = x1(i1+1) - x1(i1) 471 dx2 = x2(i2+1) - x2(i2) 472 z(1) = y(i1,i2) 473 z(2) = y(i1+1,i2) 474 z(3) = y(i1+1,i2+1) 475 z(4) = y(i1,i2+1) 476 z(5) = der1_x1(i1,i2) * dx1 477 z(6) = der1_x1(i1+1,i2) * dx1 478 z(7) = der1_x1(i1+1,i2+1) * dx1 479 z(8) = der1_x1(i1,i2+1) * dx1 480 z(9) = der1_x2(i1,i2) * dx2 481 z(10) = der1_x2(i1+1,i2) * dx2 482 z(11) = der1_x2(i1+1,i2+1) * dx2 483 z(12) = der1_x2(i1,i2+1) * dx2 484 z(13) = der2_x1x2(i1,i2) * dx1 * dx2 485 z(14) = der2_x1x2(i1+1,i2) * dx1 * dx2 486 z(15) = der2_x1x2(i1+1,i2+1) * dx1 * dx2 487 z(16) = der2_x1x2(i1,i2+1) * dx1 * dx2 488 z = matmul(wt,z) 489 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 spl_c(:,:,n1,:) = 0 495 spl_c(:,:,:,n2) = 0 496 spl_c(1,1,n1,:) = y(n1,:) 497 spl_c(1,1,:,n2) = y(:,n2) 498 499 end subroutine spline_bicubic
m_splines/spline_c [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
spline_c
FUNCTION
Computes the spline of a complex function.
INPUTS
nomega_lo = number of point in the non regular grid (e.g. !logarithmic) nomega_li = number of point in the regular grid on which the spline is computed omega_lo = value of freq on the 1st grid omega_li = value of freq on the 2nd grid tospline_lo = function on the 1st grid
OUTPUT
splined_lo = spline (on the 2nd grid)
SOURCE
523 subroutine spline_c( nomega_lo, nomega_li, omega_lo, omega_li, splined_li, tospline_lo) 524 525 !Arguments -------------------------------------------- 526 !scalars 527 integer, intent(in) :: nomega_lo, nomega_li 528 real(dp), intent(in) :: omega_lo(nomega_lo) 529 real(dp), intent(in) :: omega_li(nomega_li) 530 complex(dpc), intent(in) :: tospline_lo(nomega_lo) 531 complex(dpc), intent(out) :: splined_li(nomega_li) 532 533 !Local variables--------------------------------------- 534 !scalars 535 complex(dpc) :: ybcbeg, ybcend 536 complex(dpc), allocatable :: ysplin2_lo(:) 537 538 ybcbeg=czero 539 ybcend=czero 540 541 ABI_MALLOC(ysplin2_lo,(nomega_lo)) 542 call spline_complex(omega_lo, tospline_lo, nomega_lo, ybcbeg, ybcend, ysplin2_lo) 543 call splint_complex( nomega_lo, omega_lo, tospline_lo,ysplin2_lo, nomega_li, omega_li, splined_li) 544 ABI_FREE(ysplin2_lo) 545 546 end subroutine spline_c
m_splines/spline_complex [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
spline_complex
FUNCTION
spline_complex interfaces the usual spline routine in a case of a complex function
INPUTS
Input, integer N, the number of data points; N must be at least 2.
In the special case where N = 2 and IBCBEG = IBCEND = 0, the
spline will actually be linear.
Input, double precision T(N), the knot values, that is, the points where data
is specified. The knot values should be distinct, and increasing.
Input, complex Y(N), the data values to be interpolated.
Input, complex YBCBEG, YBCEND, the values to be used in the boundary
conditions if IBCBEG or IBCEND is equal to 1 or 2.
OUTPUT
Output, complex YPP(N), the second derivatives of the cubic spline.
SOURCE
577 subroutine spline_complex( t, y, n, ybcbeg, ybcend, ypp ) 578 579 integer, intent(in) :: n 580 real(dp), intent(in) :: t(n) 581 complex(dpc), intent(in) :: y(n) 582 complex(dpc), intent(in) :: ybcbeg 583 complex(dpc), intent(in) :: ybcend 584 complex(dpc), intent(out) :: ypp(n) 585 586 real(dp), allocatable :: y_r(:) 587 real(dp) :: ybcbeg_r 588 real(dp) :: ybcend_r 589 real(dp), allocatable :: ypp_r(:) 590 real(dp), allocatable :: y_i(:) 591 real(dp) :: ybcbeg_i 592 real(dp) :: ybcend_i 593 real(dp), allocatable :: ypp_i(:) 594 595 ABI_MALLOC(y_r,(n)) 596 ABI_MALLOC(ypp_r,(n)) 597 ABI_MALLOC(y_i,(n)) 598 ABI_MALLOC(ypp_i,(n)) 599 y_r=real(y) 600 y_i=aimag(y) !vz_d 601 ybcbeg_r=real(ybcbeg) 602 ybcbeg_i=aimag(ybcbeg) !vz_d 603 ybcend_r=real(ybcend) 604 ybcend_i=aimag(ybcend) !vz_d 605 call spline( t, y_r, n, ybcbeg_r, ybcend_r, ypp_r ) 606 call spline( t, y_i, n, ybcbeg_i, ybcend_i, ypp_i ) 607 ypp=cmplx(ypp_r,ypp_i) 608 ABI_FREE(y_r) 609 ABI_FREE(ypp_r) 610 ABI_FREE(y_i) 611 ABI_FREE(ypp_i) 612 613 end subroutine spline_complex
m_splines/spline_integrate [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
spline_integrate
FUNCTION
Calculates an integral using cubic spline interpolation.
INPUTS
npts= number of grid points of input mesh dx= step of input mesh integrand= function on input mesh
OUTPUT
integral= integral of the input function
SOURCE
780 subroutine spline_integrate(integral,npts,dx,integrand) 781 782 integer,intent(in) :: npts 783 real(dp),intent(out) :: integral 784 real(dp),intent(in) :: dx,integrand(npts) 785 786 integer :: ix 787 real(dp) :: ptmp,sf(npts),sf_der2(npts),sf_mesh(npts),utmp(npts) 788 789 ! Prepare mesh 790 forall (ix=1:npts) sf_mesh(ix) = (ix - 1) * dx 791 792 ! Calculate second derivative of integrand (adapted from Numercial Recipes) 793 sf_der2(1) = zero 794 sf_der2(npts) = zero 795 utmp(1) = zero 796 797 do ix=2,npts-1 798 ptmp = half * sf_der2(ix-1) + two 799 sf_der2(ix) = (half - one) / ptmp 800 utmp(ix) = (three * (integrand(ix+1) + integrand(ix-1) - & 801 & two*integrand(ix)) / (dx**2) - half * utmp(ix-1)) / ptmp 802 end do 803 do ix=npts-1,1,-1 804 sf_der2(ix) = sf_der2(ix) * sf_der2(ix+1) + utmp(ix) 805 end do 806 807 ! Actually calculate integral 808 sf(:) = integrand(:) * dx 809 integral = (sf(1) + sf(npts)) / 2.0_dp - & 810 & (sf_der2(1) + sf_der2(npts)) / 24.0_dp + & 811 & sum(sf(2:npts-1)) - sum(sf_der2(2:npts-1)) / 12.0_dp 812 813 end subroutine spline_integrate
m_splines/splint [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
splint
FUNCTION
Compute spline interpolation. There is no hypothesis about the spacing of the input grid points.
INPUTS
nspline: number of grid points of input mesh xspline(nspline): input mesh yspline(nspline): function on input mesh ysplin2(nspline): second derivative of yspline on input mesh nfit: number of points of output mesh xfit(nfit): output mesh
OUTPUT
yfit(nfit): function on output mesh
[ierr]=A non-zero value is used to signal that some points in xfit exceed xspline(nspline).
The input value is incremented by the number of such points.
SOURCE
641 subroutine splint(nspline,xspline,yspline,ysplin2,nfit,xfit,yfit,ierr) 642 643 integer, intent(in) :: nfit, nspline 644 integer,optional,intent(out) :: ierr 645 real(dp), intent(in) :: xspline(nspline) 646 real(dp), intent(in) :: yspline(nspline) 647 real(dp), intent(in) :: ysplin2(nspline) 648 real(dp), intent(in) :: xfit(nfit) 649 real(dp), intent(out) :: yfit(nfit) 650 651 !local 652 integer :: left,i,k,right,my_err 653 real(dp) :: delarg,invdelarg,aa,bb 654 655 !source 656 657 my_err=0 658 659 left = 1 660 do i=1, nfit 661 yfit(i)=0.d0 ! Initialize for the unlikely event that rmax exceed r(mesh) 662 ! 663 do k=left+1, nspline 664 if(xspline(k) >= xfit(i)) then 665 if(xspline(k-1) <= xfit(i)) then 666 right = k 667 left = k-1 668 else 669 if (k-1.eq.1 .and. i.eq.1) then 670 ABI_ERROR('xfit(1) < xspline(1)') 671 !my_err=my_err+1 672 !exit 673 else 674 ABI_ERROR('xfit not properly ordered') 675 end if 676 end if 677 delarg= xspline(right) - xspline(left) 678 invdelarg= 1.0d0/delarg 679 aa= (xspline(right)-xfit(i))*invdelarg 680 bb= (xfit(i)-xspline(left))*invdelarg 681 682 yfit(i) = aa*yspline(left) + bb*yspline(right) & 683 & +( (aa*aa*aa-aa)*ysplin2(left) + & 684 & (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0d0 685 exit 686 end if 687 end do ! k 688 ! 689 if (k==nspline+1) my_err=my_err+1 ! xfit not found 690 end do ! i 691 692 if (PRESENT(ierr)) ierr=my_err 693 694 end subroutine splint
m_splines/splint_complex [ Functions ]
[ Top ] [ m_splines ] [ Functions ]
NAME
splint_complex
FUNCTION
Interface to the usual splint to compute *complex* spline interpolation. There is no hypothesis about the spacing of the input grid points.
INPUTS
nspline: number of grid points of input mesh xspline(nspline): input mesh yspline(nspline): complex function on input mesh ysplin2(nspline): second derivative of yspline on input mesh nfit: number of points of output mesh xfit(nfit): output mesh
OUTPUT
yfit(nfit): complex function on output mesh
SOURCE
720 subroutine splint_complex (nspline,xspline,yspline,ysplin2,nfit,xfit,yfit) 721 722 integer, intent(in) :: nfit, nspline 723 real(dp), intent(in) :: xspline(nspline) 724 complex(dpc), intent(in) :: yspline(nspline) 725 complex(dpc), intent(in) :: ysplin2(nspline) 726 real(dp), intent(in) :: xfit(nfit) 727 complex(dpc), intent(out) :: yfit(nfit) 728 729 real(dp), allocatable :: ysplin2_r(:) 730 real(dp), allocatable :: ysplin2_i(:) 731 real(dp), allocatable :: yspline_r(:) 732 real(dp), allocatable :: yspline_i(:) 733 real(dp), allocatable :: yfit_r(:) 734 real(dp), allocatable :: yfit_i(:) 735 736 ABI_MALLOC(yspline_r,(nspline)) 737 ABI_MALLOC(yspline_i,(nspline)) 738 ABI_MALLOC(ysplin2_r,(nspline)) 739 ABI_MALLOC(ysplin2_i,(nspline)) 740 ABI_MALLOC(yfit_r,(nfit)) 741 ABI_MALLOC(yfit_i,(nfit)) 742 743 !local 744 745 !source 746 yspline_r=real(yspline) 747 yspline_i=aimag(yspline) !vz_d 748 ysplin2_r=real(ysplin2) 749 ysplin2_i=aimag(ysplin2) !vz_d 750 call splint (nspline,xspline,yspline_r,ysplin2_r,nfit,xfit,yfit_r) 751 call splint (nspline,xspline,yspline_i,ysplin2_i,nfit,xfit,yfit_i) 752 yfit=cmplx(yfit_r,yfit_i) 753 ABI_FREE(yspline_r) 754 ABI_FREE(yspline_i) 755 ABI_FREE(ysplin2_r) 756 ABI_FREE(ysplin2_i) 757 ABI_FREE(yfit_r) 758 ABI_FREE(yfit_i) 759 760 end subroutine splint_complex