TABLE OF CONTENTS


ABINIT/m_splines [ Modules ]

[ Top ] [ Modules ]

NAME

  m_splines

FUNCTION

  This module contains routines for spline interpolation.

COPYRIGHT

  Copyright (C) 2010-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_splines
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32 
33  implicit none
34 
35  public :: splfit
36  public :: spline
37  public :: spline_bicubic
38  public :: spline_c
39  public :: spline_complex
40  public :: spline_integrate
41  public :: splint
42  public :: splint_complex
43 
44  !FIXME deprecated
45  public :: intrpl
46 
47 ! *************************************************************************
48 
49 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.

PARENTS

CHILDREN

SOURCE

 994 SUBROUTINE INTRPL(L,X,Y,N,U,V,dv,dv2,ideriv)
 995 
 996 
 997 !This section has been created automatically by the script Abilint (TD).
 998 !Do not modify the following lines by hand.
 999 #undef ABI_FUNC
1000 #define ABI_FUNC 'INTRPL'
1001 !End of the abilint section
1002 
1003       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1004       IMPLICIT INTEGER(I-N)
1005 !
1006       PARAMETER (NQQ=12000)
1007 
1008       COMMON/QQ/ QQ(4,NQQ)
1009       DIMENSION X(L),Y(L),U(N),V(N),DV(NQQ),DV2(NQQ)
1010       EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3)
1011       REAL*8 M1,M2,M3,M4,M5
1012       EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5),&
1013      & (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
1014 !
1015 !     PRELIMINARY PROCESSING
1016 
1017       L0=L
1018       LM1=L0-1
1019       LM2=LM1-1
1020       LP1=L0+1
1021       N0=N
1022       IF(N0.GT.NQQ) THEN
1023           NQQV=NQQ
1024           write(std_out,2089) NQQV,N0
1025 !          CALL EXIT
1026       END IF
1027       IF(LM2.LT.0) GO TO 90
1028       IF(N0.LE.0) GO TO 91
1029       DO 11 I=2,L0
1030 
1031 !     IF(X(I-1)-X(I))11,95,96
1032       IF(X(I-1)-X(I).EQ.0.0D0) GO TO 95
1033       IF(X(I-1)-X(I).GT.0.0D0) GO TO 96
1034    11 CONTINUE
1035       IPV=0
1036 !
1037 !***  MAIN LOOP
1038       FINT=0.0D0
1039       DO 80 K=1,N0
1040       UK=U(K)
1041 !
1042 !***  ROUTINE TO LOCATE THE DESIRED POINT
1043        IF(UK.GE.X(L0)) GO TO 26
1044       IF(UK.LT.X(1)) GO TO 25
1045       IMN=2
1046       IMX=L0
1047    21 I=(IMN+IMX)/2
1048       IF(UK.GE.X(I)) GO TO 23
1049       IMX=I
1050       GO TO 24
1051    23 IMN=I+1
1052    24 IF(IMX.GT.IMN) GO TO 21
1053       I=IMX
1054       GO TO 30
1055    25 I=1
1056       GO TO 30
1057    26 I=LP1
1058       GO TO 30
1059 !
1060 !***  CHECK IF I=IPV
1061    30 IF(I.EQ.IPV) GO TO 70
1062       IPV=I
1063 !
1064 !***  ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND TO
1065 !***  ESTIMATE THEM IF NECESSARY
1066       J=I
1067       IF(J.EQ.1) J=2
1068       IF(J.EQ.LP1) J=L0
1069       X3=X(J-1)
1070       Y3=Y(J-1)
1071       X4=X(J)
1072       Y4=Y(J)
1073       A3=X4-X3
1074       M3=(Y4-Y3)/A3
1075       IF(LM2.EQ.0) GO TO 43
1076       IF(J.EQ.2) GO TO 41
1077       X2=X(J-2)
1078       Y2=Y(J-2)
1079       A2=X3-X2
1080       M2=(Y3-Y2)/A2
1081       IF(J.EQ.L0) GO TO 42
1082    41 X5=X(J+1)
1083       Y5=Y(J+1)
1084       A4=X5-X4
1085       M4=(Y5-Y4)/A4
1086       IF(J.EQ.2) M2=M3+M3-M4
1087       GO TO 45
1088    42 M4=M3+M3-M2
1089       GO TO 45
1090    43 M2=M3
1091    45 IF(J.LE.3) GO TO 46
1092       A1=X2-X(J-3)
1093       M1=(Y2-Y(J-3))/A1
1094       GO TO 47
1095    46 M1=M2+M2-M3
1096    47 IF(J.GE.LM1) GO TO 48
1097       A5=X(J+2)-X5
1098       M5=(Y(J+2)-Y5)/A5
1099       GO TO 50
1100    48 M5=M4+M4-M3
1101 !
1102 !***  NUMERICAL DIFFERENTIATION
1103    50 IF(I.EQ.LP1) GO TO 52
1104       W2=ABS(M4-M3)
1105       W3=ABS(M2-M1)
1106       SW=W2+W3
1107       IF(SW.NE.0.0) GO TO 51
1108       W2=0.5D0
1109       W3=0.5D0
1110       SW=1.0D0
1111    51 T3=(W2*M2+W3*M3)/SW
1112       IF(I.EQ.1) GO TO 54
1113    52 W3=ABS(M5-M4)
1114       W4=ABS(M3-M2)
1115       SW=W3+W4
1116       IF(SW.NE.0.0) GO TO 53
1117       W3=0.5D0
1118       W4=0.5D0
1119       SW=1.0D0
1120    53 T4=(W3*M3+W4*M4)/SW
1121       IF(I.NE.LP1) GO TO 60
1122       T3=T4
1123       SA=A2+A3
1124       T4=0.5D0*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA))
1125       X3=X4
1126       Y3=Y4
1127       A3=A2
1128       M3=M4
1129       GO TO 60
1130    54 T4=T3
1131       SA=A3+A4
1132       T3=0.5D0*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA))
1133       X3=X3-A4
1134       Y3=Y3-M2*A4
1135       A3=A4
1136       M3=M2
1137 !
1138 !***  COMPUTATION OF THE POLYNOMIAL
1139    60 Q2=(2.0D0*(M3-T3)+M3-T4)/A3
1140       Q3=(-M3-M3+T3+T4)/(A3*A3)
1141    70 DX=UK-P0
1142       V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3))
1143 
1144       IF(IDERIV.EQ.0) GO TO 80
1145       DV(K)=Q1+DX*(2.0D0*Q2+DX*3.0D0*Q3)
1146       DV2(k)=6.0D0*Q3*DX+2.d0*Q2
1147       QQ(1,K)=Q0
1148       QQ(2,K)=Q1
1149       QQ(3,K)=Q2
1150       QQ(4,K)=Q3
1151    80  CONTINUE
1152       RETURN
1153 !
1154 !***  ERROR EXIT
1155    90 write(std_out,2090)
1156       GO TO 99
1157    91 write(std_out,2091)
1158       GO TO 99
1159    95 write(std_out,2095)
1160       GO TO 97
1161    96 write(std_out,2096)
1162    97 write(std_out,2097)I,X(I)
1163    99 write(std_out,2099) L0,N0
1164       RETURN
1165 !
1166 !***  FORMAT STATEMENTS
1167  2089  FORMAT( 'WARNING ERROR IN INTRPL. MAX ALLOWED VALUE OF N0 IS',&
1168      & I3,' HERE N0 IS',I3)
1169  2090  FORMAT(1X/' N = 1 OR LESS.'/)
1170  2091  FORMAT(1X/' N = 0 OR LESS.'/)
1171  2095  FORMAT(1X/' IDENTICAL X VALUES.'/)
1172  2096  FORMAT(1X/' X VALUES OUT OF SEQUENCE.'/)
1173  2097  FORMAT(4X,'I =',I7,10X,6X,'X(I) =',E12.3)
1174  2099  FORMAT(4X,'L =',I7,10X,3X,'N =',I7/ &
1175      & ' ERROR DETECTED IN ROUTINE INTRPL')
1176 !
1177 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 delarg) 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)

PARENTS

      getnel,m_pawpwij,mkffnl,pawgylmg,psp8lo

CHILDREN

SOURCE

 92 subroutine splfit(arg,derfun,fun,ider,newarg,newfun,numarg,numnew)
 93 
 94 
 95 !This section has been created automatically by the script Abilint (TD).
 96 !Do not modify the following lines by hand.
 97 #undef ABI_FUNC
 98 #define ABI_FUNC 'splfit'
 99 !End of the abilint section
100 
101  implicit none
102 
103  integer, intent(in) :: ider,numarg,numnew
104  real(dp), intent(in) :: arg(numarg),fun(numarg,2),newarg(numnew)
105  real(dp), intent(out) :: derfun(numnew)
106  real(dp), intent(inout) :: newfun(numnew)
107 
108  integer :: i,jspl
109  real(dp) :: argmin,delarg,d,aa,bb,cc,dd
110  character(len=500) :: msg
111 
112 !argmin is smallest x value in spline fit; delarg is uniform spacing of spline argument
113  argmin=arg(1)
114  delarg=(arg(numarg)-argmin)/dble(numarg-1)
115 
116  if(delarg<tol12)then
117    write(msg,'(a,es16.8)') ' delarg should be strictly positive, while delarg= ',delarg
118    MSG_ERROR(msg)
119  endif
120 
121  jspl=-1
122 
123 !Do one loop for no grads, other for grads:
124  if (ider==0) then
125 ! Spline index loop for no grads:
126   do i=1,numnew
127    if (newarg(i).ge.arg(numarg)) then
128 ! MJV 6/4/2009 FIXME this message is never used
129     write(msg,1000)char(10),i,newarg(i), &
130 &     jspl,char(10),numarg,arg(numarg),char(10),char(10),char(10)
131 1000 format(a1,' splfit: for arg number',i8,2x,'of value', &
132 &    1p,e12.4,1x,'jspl=',i8,a1,' is >= numarg=',i8,  &
133 &    '  (max arg(numarg)=',e12.4,')',a1,             &
134 &    ' Means function values are being requested outside',       &
135 &    ' range of data.',a1,' Function and slope will be set to',  &
136 &    ' values at upper end of data.',a1)
137 
138     newfun(i)=fun(numarg,1)
139 
140    else if (newarg(i).le.arg(1)) then
141     newfun(i)=fun(1,1)
142 
143    else
144     jspl=1+int((newarg(i)-argmin)/delarg)
145     d=newarg(i)-arg(jspl)
146     bb = d/delarg
147     aa = 1.0d0-bb
148     cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0)
149     dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0)
150     newfun(i)=aa*fun(jspl,1)+bb*fun(jspl+1,1)+cc*fun(jspl,2)+dd*fun(jspl+1,2)
151    end if
152   enddo
153 
154  else if(ider==1)then
155 
156 ! Spline index loop includes grads:
157   do i=1,numnew
158 
159    if (newarg(i).ge.arg(numarg)) then
160     newfun(i)=fun(numarg,1)
161     derfun(i)=0.0d0
162 
163    else if (newarg(i).le.arg(1)) then
164     newfun(i)=fun(1,1)
165     derfun(i)=0.0d0
166 
167    else
168 
169 !   cubic spline interpolation:
170     jspl=1+int((newarg(i)-arg(1))/delarg)
171     d=newarg(i)-arg(jspl)
172     bb = d/delarg
173     aa = 1.0d0-bb
174     cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0)
175     dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0)
176     newfun(i)=aa*fun(jspl,1)+bb*fun(jspl+1,1)+cc*fun(jspl,2)+dd*fun(jspl+1,2)
177 !   spline fit to first derivative:
178 !   note correction of Numerical Recipes sign error
179     derfun(i) = (fun(jspl+1,1)-fun(jspl,1))/delarg +    &
180 &      (-(3.d0*aa**2-1.d0)*fun(jspl,2)+                 &
181 &        (3.d0*bb**2-1.d0)*fun(jspl+1,2)) * delarg/6.0d0
182 
183           end if
184   enddo
185 
186  else if (ider==2) then
187 
188   do i=1,numnew
189 
190    if (newarg(i).ge.arg(numarg)) then
191     derfun(i)=0.0d0
192 
193    else if (newarg(i).le.arg(1)) then
194     derfun(i)=0.0d0
195 
196    else
197 
198 !   cubic spline interpolation:
199     jspl=1+int((newarg(i)-argmin)/delarg)
200     d=newarg(i)-arg(jspl)
201     bb = d/delarg
202     aa = 1.0d0-bb
203 !   second derivative of spline (piecewise linear function)
204     derfun(i) = aa*fun(jspl,2)+bb*fun(jspl+1,2)
205 
206    end if
207   enddo
208 
209  end if
210 
211 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 ...

PARENTS

      atomden,calc_sigc_cd,calc_sigc_pole_cd,cc_derivatives,denfgr,get_tau_k
      init_occ_ent,integrho,m_atom,m_dens,m_entropyDMFT,m_ewald,m_paw_slater
      m_special_funcs,m_splines,outscfcv,pawinit,predict_string,psp10in
      psp10nl,psp11nl,psp1cc,psp1in,psp1nl,psp2in,psp2nl,psp3in,psp3nl,psp4cc
      psp5in,psp5nl,psp6cc,psp6in,psp8in,psp8lo,psp8nl,psp9in
      random_stopping_power,spline_paw_fncs,upf2abinit,vso_realspace_local

CHILDREN

SOURCE

253 subroutine spline( t, y, n, ybcbeg, ybcend, ypp )
254 
255 !*******************************************************************************
256 !
257 !  Discussion:
258 !
259 !    For data interpolation, the user must call SPLINE_CUBIC_SET to
260 !    determine the second derivative data, passing in the data to be
261 !    interpolated, and the desired boundary conditions.
262 !
263 !    The data to be interpolated, plus the SPLINE_CUBIC_SET output,
264 !    defines the spline.  The user may then call SPLINE_CUBIC_VAL to
265 !    evaluate the spline at any point.
266 !
267 !    The cubic spline is a piecewise cubic polynomial.  The intervals
268 !    are determined by the "knots" or abscissas of the data to be
269 !    interpolated.  The cubic spline has continous first and second
270 !    derivatives over the entire interval of interpolation.
271 !
272 !    For any point T in the interval T(IVAL), T(IVAL+1), the form of
273 !    the spline is
274 !
275 !      SPL(T) = A(IVAL)
276 !             + B(IVAL) * ( T - T(IVAL) )
277 !             + C(IVAL) * ( T - T(IVAL) )**2
278 !             + D(IVAL) * ( T - T(IVAL) )**3
279 !
280 !    If we assume that we know the values Y(*) and YPP(*), which represent
281 !    the values and second derivatives of the spline at each knot, then
282 !    the coefficients can be computed as:
283 !
284 !      A(IVAL) = Y(IVAL)
285 !      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
286 !        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
287 !      C(IVAL) = YPP(IVAL) / 2
288 !      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
289 !
290 !    Since the first derivative of the spline is
291 !
292 !      SPL'(T) =     B(IVAL)
293 !              + 2 * C(IVAL) * ( T - T(IVAL) )
294 !              + 3 * D(IVAL) * ( T - T(IVAL) )**2,
295 !
296 !    the requirement that the first derivative be continuous at interior
297 !    knot I results in a total of N-2 equations, of the form:
298 !
299 !      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
300 !      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
301 !
302 !    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
303 !
304 !      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
305 !      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
306 !      + YPP(IVAL-1) * H(IVAL-1)
307 !      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
308 !      =
309 !      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
310 !      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
311 !
312 !    or
313 !
314 !      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
315 !      + YPP(IVAL) * H(IVAL)
316 !      =
317 !      6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
318 !      - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
319 !
320 !    Boundary conditions must be applied at the first and last knots.
321 !    The resulting tridiagonal system can be solved for the YPP values.
322 !
323 !  Modified:
324 !
325 !    07 February 1999
326 !    28 November 2004 XGonze : double precision
327 !                              make arguments similar to the Numeric Recipes routine
328 !                              also use algorithmics similar to the Numeric Recipes routine
329 !
330 !  Author:
331 !
332 !    John Burkardt
333 !    (XGonze got it from http://www.psc.edu/~burkardt/src/spline/spline.html)
334 !
335 !  Parameters:
336 !
337 !    Input, integer N, the number of data points; N must be at least 2.
338 !    In the special case where N = 2 and IBCBEG = IBCEND = 0, the
339 !    spline will actually be linear.
340 !
341 !    Input, double precision T(N), the knot values, that is, the points where data
342 !    is specified.  The knot values should be distinct, and increasing.
343 !
344 !    Input, double precision Y(N), the data values to be interpolated.
345 !
346 !    Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
347 !    conditions if IBCBEG or IBCEND is equal to 1 or 2.
348 !
349 !    Output, double precision YPP(N), the second derivatives of the cubic spline.
350 !
351 !    Work space, double precision DIAG(N) - should be removed ...
352 !
353 !
354 !    XG041127 : In the initial implementation, one had the control on
355 !     IBCBEG and IBCEND. Now, they are determined by the values
356 !     of YBCBEG, YBCEND. Option 2 has been disabled.
357 !
358 !    Input, integer IBCBEG, left boundary condition flag:
359 !
360 !      0: the spline should be a quadratic over the first interval;
361 !      1: the first derivative at the left endpoint should be YBCBEG;
362 !      2: the second derivative at the left endpoint should be YBCBEG.
363 !
364 !    Input, integer IBCEND, right boundary condition flag:
365 !
366 !      0: the spline should be a quadratic over the last interval;
367 !      1: the first derivative at the right endpoint should be YBCEND;
368 !      2: the second derivative at the right endpoint should be YBCEND.
369 
370 !This section has been created automatically by the script Abilint (TD).
371 !Do not modify the following lines by hand.
372 #undef ABI_FUNC
373 #define ABI_FUNC 'spline'
374 !End of the abilint section
375 
376   implicit none
377 
378   integer, intent(in) :: n
379   real(dp), intent(in) :: t(n)
380   real(dp), intent(in) :: y(n)
381   real(dp), intent(in) :: ybcbeg
382   real(dp), intent(in) :: ybcend
383 
384   real(dp), intent(out) :: ypp(n)
385 
386   integer :: ibcbeg
387   integer :: ibcend
388   integer :: i,k
389   real(dp) :: ratio,pinv
390   real(dp), allocatable :: tmp(:)
391 !
392 !  Check.
393 !
394   if ( n <= 1 ) then
395     write(std_out,* ) ' '
396     write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!'
397     write(std_out,* ) '  The number of knots must be at least 2.'
398     write(std_out,* ) '  The input value of N = ', n
399     MSG_ERROR("Fatal error")
400   end if
401 
402   ABI_ALLOCATE(tmp,(n))
403 
404   do i = 1, n-1
405     if ( t(i) >= t(i+1) ) then
406       write(std_out,* ) ' '
407       write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!'
408       write(std_out,* ) '  The knots must be strictly increasing, but'
409       write(std_out,* ) '  T(',  i,') = ', t(i)
410       write(std_out,* ) '  T(',i+1,') = ', t(i+1)
411       MSG_ERROR("Fatal error")
412     end if
413   end do
414 !
415 !  XG041127
416   ibcbeg=1 ; ibcend=1
417   if(ybcbeg>1.0d+30)ibcbeg=0
418   if(ybcend>1.0d+30)ibcend=0
419 !
420 !  Set the first and last equations.
421 !
422   if ( ibcbeg == 0 ) then
423     ypp(1) = 0.d0
424     tmp(1) = 0.d0
425   else if ( ibcbeg == 1 ) then
426     ypp(1) = -0.5d0
427     tmp(1) = (3.d0/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg)
428   end if
429   if ( ibcend == 0 ) then
430     ypp(n) = 0.d0
431     tmp(n) = 0.d0
432   else if ( ibcend == 1 ) then
433     ypp(n) = 0.5d0
434     tmp(n) = (3.d0/(t(n)-t(n-1)))*(ybcend-(y(n)-y(n-1))/(t(n)-t(n-1)))
435   end if
436 
437 !
438 !  Set the intermediate equations.
439 !
440   do i=2,n-1
441    ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1))
442    pinv = 1.0d0/(ratio*ypp(i-1) + 2.0d0)
443    ypp(i) = (ratio-1.0d0)*pinv
444    tmp(i)=(6.0d0*((y(i+1)-y(i))/(t(i+1)-t(i))-(y(i)-y(i-1)) &
445 &    /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv
446    if (abs(tmp(i))<1.d5*tiny(0.d0)) tmp(i)=0.d0   !MT20050927
447   enddo
448 
449 ! Solve the equations
450   ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0d0)
451   do k=n-1,1,-1
452    ypp(k)=ypp(k)*ypp(k+1)+tmp(k)
453   enddo
454 
455   ABI_DEALLOCATE(tmp)
456 
457   return
458 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.

PARENTS

      m_xc_vdw

CHILDREN

SOURCE

493 subroutine spline_bicubic(n1,n2,x1,x2,y,der1_x1,der1_x2,der2_x1x2,spl_c)
494 
495 
496 !This section has been created automatically by the script Abilint (TD).
497 !Do not modify the following lines by hand.
498 #undef ABI_FUNC
499 #define ABI_FUNC 'spline_bicubic'
500 !End of the abilint section
501 
502   implicit none
503 
504   integer,intent(in)  :: n1,n2
505   real(dp),intent(in) :: x1(n1),x2(n2),y(n1,n2)
506   real(dp),intent(in) :: der1_x1(n1,n2),der1_x2(n1,n2),der2_x1x2(n1,n2)
507   real(dp),intent(out):: spl_c(4,4,n1,n2)
508 
509   integer :: i1,i2
510   real(dp) :: dx1,dx2,wt(16,16),z(16)
511 
512   data wt /1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4, &
513 &          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, &
514 &          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, &
515 &          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, &
516 &          -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, &
517 &          -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, &
518 &          -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/
519 
520   ! Set coefficients for i1<n1 and i2<n2
521   do i2 = 1,n2-1
522     do i1 = 1,n1-1
523       dx1 = x1(i1+1) - x1(i1)
524       dx2 = x2(i2+1) - x2(i2)
525       z(1)  = y(i1,i2)
526       z(2)  = y(i1+1,i2)
527       z(3)  = y(i1+1,i2+1)
528       z(4)  = y(i1,i2+1)
529       z(5)  = der1_x1(i1,i2) * dx1
530       z(6)  = der1_x1(i1+1,i2) * dx1
531       z(7)  = der1_x1(i1+1,i2+1) * dx1
532       z(8)  = der1_x1(i1,i2+1) * dx1
533       z(9)  = der1_x2(i1,i2) * dx2
534       z(10) = der1_x2(i1+1,i2) * dx2
535       z(11) = der1_x2(i1+1,i2+1) * dx2
536       z(12) = der1_x2(i1,i2+1) * dx2
537       z(13) = der2_x1x2(i1,i2) * dx1 * dx2
538       z(14) = der2_x1x2(i1+1,i2) * dx1 * dx2
539       z(15) = der2_x1x2(i1+1,i2+1) * dx1 * dx2
540       z(16) = der2_x1x2(i1,i2+1) * dx1 * dx2
541       z = matmul(wt,z)
542       spl_c(:,:,i1,i2) = reshape(z,(/4,4/),order=(/2,1/))
543     end do
544   end do
545 
546 ! Set coefficients for i1=n1 and i2=n2 (valid only at the border)
547   spl_c(:,:,n1,:) = 0
548   spl_c(:,:,:,n2) = 0
549   spl_c(1,1,n1,:) = y(n1,:)
550   spl_c(1,1,:,n2) = y(:,n2)
551 
552 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)

PARENTS

      m_green

CHILDREN

SOURCE

581 subroutine spline_c( nomega_lo, nomega_li, omega_lo, omega_li, splined_li, tospline_lo)
582 
583 
584 !This section has been created automatically by the script Abilint (TD).
585 !Do not modify the following lines by hand.
586 #undef ABI_FUNC
587 #define ABI_FUNC 'spline_c'
588 !End of the abilint section
589 
590  implicit none
591 
592 !Arguments --------------------------------------------
593 !scalars
594  integer, intent(in) :: nomega_lo, nomega_li
595  real(dp), intent(in) :: omega_lo(nomega_lo)
596  real(dp), intent(in) :: omega_li(nomega_li)
597  complex(dpc), intent(in) :: tospline_lo(nomega_lo)
598  complex(dpc), intent(out) :: splined_li(nomega_li)
599 
600 !Local variables---------------------------------------
601 !scalars
602  complex(dpc) :: ybcbeg, ybcend
603  complex(dpc), allocatable :: ysplin2_lo(:)
604 
605  ybcbeg=czero
606  ybcend=czero
607 
608  ABI_ALLOCATE(ysplin2_lo,(nomega_lo))
609  call spline_complex(omega_lo, tospline_lo, nomega_lo, ybcbeg, ybcend, ysplin2_lo)
610  call splint_complex( nomega_lo, omega_lo, tospline_lo,ysplin2_lo, nomega_li, omega_li, splined_li)
611  ABI_DEALLOCATE(ysplin2_lo)
612 
613 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.

PARENTS

      m_paw_dmft,m_splines

CHILDREN

SOURCE

649 subroutine spline_complex( t, y, n, ybcbeg, ybcend, ypp )
650 
651 
652 !This section has been created automatically by the script Abilint (TD).
653 !Do not modify the following lines by hand.
654 #undef ABI_FUNC
655 #define ABI_FUNC 'spline_complex'
656 !End of the abilint section
657 
658  implicit none
659 
660  integer, intent(in) :: n
661  real(dp), intent(in) :: t(n)
662  complex(dpc), intent(in) :: y(n)
663  complex(dpc), intent(in) :: ybcbeg
664  complex(dpc), intent(in) :: ybcend
665  complex(dpc), intent(out) :: ypp(n)
666 
667  real(dp), allocatable :: y_r(:)
668  real(dp) :: ybcbeg_r
669  real(dp) :: ybcend_r
670  real(dp), allocatable :: ypp_r(:)
671  real(dp), allocatable :: y_i(:)
672  real(dp) :: ybcbeg_i
673  real(dp) :: ybcend_i
674  real(dp), allocatable :: ypp_i(:)
675 
676  ABI_ALLOCATE(y_r,(n))
677  ABI_ALLOCATE(ypp_r,(n))
678  ABI_ALLOCATE(y_i,(n))
679  ABI_ALLOCATE(ypp_i,(n))
680  y_r=real(y)
681  y_i=aimag(y)    !vz_d
682  ybcbeg_r=real(ybcbeg)
683  ybcbeg_i=aimag(ybcbeg)    !vz_d
684  ybcend_r=real(ybcend)
685  ybcend_i=aimag(ybcend)    !vz_d
686  call spline( t, y_r, n, ybcbeg_r, ybcend_r, ypp_r )
687  call spline( t, y_i, n, ybcbeg_i, ybcend_i, ypp_i )
688  ypp=cmplx(ypp_r,ypp_i)
689  ABI_DEALLOCATE(y_r)
690  ABI_DEALLOCATE(ypp_r)
691  ABI_DEALLOCATE(y_i)
692  ABI_DEALLOCATE(ypp_i)
693 
694 end subroutine spline_complex

m_splines/spline_integrate [ Functions ]

[ Top ] [ m_splines ] [ Functions ]

NAME

  spline_integrate

FUNCTION

  Calculates an integral using cubic spline interpolation.

COPYRIGHT

  Copyright (C) 2010-2018 ABINIT Group (Yann Pouillon)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .
  For the initials of contributors, see
  ~abinit/doc/developers/contributors.txt .

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

PARENTS

      test_spline_integrate

CHILDREN

SOURCE

905 subroutine spline_integrate(integral,npts,dx,integrand)
906 
907 
908 !This section has been created automatically by the script Abilint (TD).
909 !Do not modify the following lines by hand.
910 #undef ABI_FUNC
911 #define ABI_FUNC 'spline_integrate'
912 !End of the abilint section
913 
914  implicit none
915 
916  integer,intent(in) :: npts
917  real(dp),intent(out) :: integral
918  real(dp),intent(in) :: dx,integrand(npts)
919 
920  integer :: ix
921  real(dp) :: ptmp,sf(npts),sf_der2(npts),sf_mesh(npts),utmp(npts)
922 
923  ! Prepare mesh
924  forall (ix=1:npts) sf_mesh(ix) = (ix - 1) * dx
925 
926  ! Calculate second derivative of integrand (adapted from Numercial Recipes)
927  sf_der2(1) = zero
928  sf_der2(npts) = zero
929  utmp(1) = zero
930 
931  do ix=2,npts-1
932   ptmp = half * sf_der2(ix-1) + two
933   sf_der2(ix) = (half - one) / ptmp
934   utmp(ix) = (three * (integrand(ix+1) + integrand(ix-1) - &
935 &  two*integrand(ix)) / (dx**2) - half * utmp(ix-1)) / ptmp
936  end do
937  do ix=npts-1,1,-1
938   sf_der2(ix) = sf_der2(ix) * sf_der2(ix+1) + utmp(ix)
939  end do
940 
941  ! Actually calculate integral
942  sf(:) = integrand(:) * dx
943  integral = (sf(1) + sf(npts)) / 2.0_dp - &
944 &           (sf_der2(1) + sf_der2(npts)) / 24.0_dp + &
945 &           sum(sf(2:npts-1)) - sum(sf_der2(2:npts-1)) / 12.0_dp
946 
947 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.

PARENTS

      atomden,calc_sigc_cd,calc_sigc_pole_cd,cc_derivatives,denfgr,get_tau_k
      m_atom,m_cut3d,m_entropyDMFT,m_paw_slater,m_special_funcs,m_splines
      outscfcv,partial_dos_fractions,predict_string,psp6cc
      random_stopping_power,spline_paw_fncs,vso_realspace_local,wvl_initro

CHILDREN

SOURCE

730 subroutine splint(nspline,xspline,yspline,ysplin2,nfit,xfit,yfit,ierr)
731 
732 
733 !This section has been created automatically by the script Abilint (TD).
734 !Do not modify the following lines by hand.
735 #undef ABI_FUNC
736 #define ABI_FUNC 'splint'
737 !End of the abilint section
738 
739  implicit none
740 
741  integer, intent(in) :: nfit, nspline
742  integer,optional,intent(out) :: ierr
743  real(dp), intent(in) :: xspline(nspline)
744  real(dp), intent(in) :: yspline(nspline)
745  real(dp), intent(in) :: ysplin2(nspline)
746  real(dp), intent(in) :: xfit(nfit)
747  real(dp), intent(out) :: yfit(nfit)
748 
749 !local
750  integer :: left,i,k,right,my_err
751  real(dp) :: delarg,invdelarg,aa,bb
752 
753 !source
754 
755  my_err=0
756 
757  left = 1
758  do i=1, nfit
759    yfit(i)=0.d0  ! Initialize for the unlikely event that rmax exceed r(mesh)
760    !
761    do k=left+1, nspline
762      if(xspline(k) >= xfit(i)) then
763        if(xspline(k-1) <= xfit(i)) then
764          right = k
765          left = k-1
766        else
767          if (k-1.eq.1 .and. i.eq.1) then
768            MSG_ERROR('xfit(1) < xspline(1)')
769            !my_err=my_err+1
770            !exit
771          else
772            MSG_ERROR('xfit not properly ordered')
773          end if
774        end if
775        delarg= xspline(right) - xspline(left)
776        invdelarg= 1.0d0/delarg
777        aa= (xspline(right)-xfit(i))*invdelarg
778        bb= (xfit(i)-xspline(left))*invdelarg
779 
780        yfit(i) = aa*yspline(left) + bb*yspline(right)    &
781 &               +( (aa*aa*aa-aa)*ysplin2(left) +         &
782 &                  (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0d0
783        exit
784      end if
785    end do ! k
786    !
787    if (k==nspline+1) my_err=my_err+1 ! xfit not found
788  end do ! i
789 
790  if (PRESENT(ierr)) ierr=my_err
791 
792 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

PARENTS

      m_paw_dmft,m_splines

CHILDREN

SOURCE

823 subroutine splint_complex (nspline,xspline,yspline,ysplin2,nfit,xfit,yfit)
824 
825 
826 !This section has been created automatically by the script Abilint (TD).
827 !Do not modify the following lines by hand.
828 #undef ABI_FUNC
829 #define ABI_FUNC 'splint_complex'
830 !End of the abilint section
831 
832  implicit none
833 
834  integer, intent(in) :: nfit, nspline
835  real(dp), intent(in) :: xspline(nspline)
836  complex(dpc), intent(in) :: yspline(nspline)
837  complex(dpc), intent(in) :: ysplin2(nspline)
838  real(dp), intent(in) :: xfit(nfit)
839  complex(dpc), intent(out) :: yfit(nfit)
840 
841  real(dp), allocatable :: ysplin2_r(:)
842  real(dp), allocatable :: ysplin2_i(:)
843  real(dp), allocatable :: yspline_r(:)
844  real(dp), allocatable :: yspline_i(:)
845  real(dp), allocatable :: yfit_r(:)
846  real(dp), allocatable :: yfit_i(:)
847 
848  ABI_ALLOCATE(yspline_r,(nspline))
849  ABI_ALLOCATE(yspline_i,(nspline))
850  ABI_ALLOCATE(ysplin2_r,(nspline))
851  ABI_ALLOCATE(ysplin2_i,(nspline))
852  ABI_ALLOCATE(yfit_r,(nfit))
853  ABI_ALLOCATE(yfit_i,(nfit))
854 
855 !local
856 
857 !source
858  yspline_r=real(yspline)
859  yspline_i=aimag(yspline)    !vz_d
860  ysplin2_r=real(ysplin2)
861  ysplin2_i=aimag(ysplin2)    !vz_d
862  call splint (nspline,xspline,yspline_r,ysplin2_r,nfit,xfit,yfit_r)
863  call splint (nspline,xspline,yspline_i,ysplin2_i,nfit,xfit,yfit_i)
864  yfit=cmplx(yfit_r,yfit_i)
865  ABI_DEALLOCATE(yspline_r)
866  ABI_DEALLOCATE(yspline_i)
867  ABI_DEALLOCATE(ysplin2_r)
868  ABI_DEALLOCATE(ysplin2_i)
869  ABI_DEALLOCATE(yfit_r)
870  ABI_DEALLOCATE(yfit_i)
871 
872 end subroutine splint_complex