TABLE OF CONTENTS
ABINIT/intrpl [ Functions ]
NAME
intrpl
FUNCTION
COPYRIGHT
Copyright (C) 2014-2017 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
28 #if defined HAVE_CONFIG_H 29 #include "config.h" 30 #endif 31 32 #include "abi_common.h" 33 34 35 !-------------------------------------------------------------------------------- 36 SUBROUTINE INTRPL(L,X,Y,N,U,V,dv,dv2,ideriv) 37 38 ! 39 ! CUBIC SPLINE INTERPOLATION 40 ! 41 42 USE DEFS_BASIS, ONLY : STD_OUT 43 44 !This section has been created automatically by the script Abilint (TD). 45 !Do not modify the following lines by hand. 46 #undef ABI_FUNC 47 #define ABI_FUNC 'INTRPL' 48 !End of the abilint section 49 50 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 51 IMPLICIT INTEGER(I-N) 52 ! 53 PARAMETER (NQQ=12000) 54 ! 55 ! DOUBLE PRECISION INTERPOLATION OF A SINGLE VALUED FUNCTION. 56 ! THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION 57 ! GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE 58 ! AND FOR A GIVEN SET OF X VALUES(ABSCISSAE),THE VALUES OF 59 ! A SINGLE VALUED FUNCTION Y=Y(X). 60 ! 61 ! THE SUBROUTINE ALSO CALCULATES FIRST DERIVATIVES DV(X) AND 62 ! SECOND DERIVATIVE DV2(X) 63 64 ! THE INPUT PARAMETERS ARE; 65 ! 66 ! L=NUMBER OF DATA POINTS 67 ! (MUST BE TWO OR GREATER) 68 ! X=ARRAY OF DIMENSION L STORING THE X VALUES 69 ! OF INPUT DATA POINTS (IN ASCENDING ORDER) 70 ! Y=ARRAY OF DIMENSION L STORING THE Y VALUES OF INPUT DATA POINTS 71 ! N=NUMBER OF POINTS AT WHICH INTERPOLATION OF THE Y-VALUES 72 ! IS REQUIRED (MUST BE 1 OR GREATER) 73 ! U=ARRAY OF DIMENSION N STORING THE X VALUES 74 ! OF THE DESIRED POINTS 75 ! 76 ! THE OUTPUT PARAMETER IS V=ARRAY OF DIMENSION N WHERE THE 77 ! INTERPOLATED Y VALUES ARE TO BE DISPLAYED 78 ! 79 80 COMMON/QQ/ QQ(4,NQQ) 81 DIMENSION X(L),Y(L),U(N),V(N),DV(NQQ),DV2(NQQ) 82 EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3) 83 REAL*8 M1,M2,M3,M4,M5 84 EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5),& 85 & (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3) 86 ! 87 ! PRELIMINARY PROCESSING 88 89 L0=L 90 LM1=L0-1 91 LM2=LM1-1 92 LP1=L0+1 93 N0=N 94 IF(N0.GT.NQQ) THEN 95 NQQV=NQQ 96 write(std_out,2089) NQQV,N0 97 ! CALL EXIT 98 END IF 99 IF(LM2.LT.0) GO TO 90 100 IF(N0.LE.0) GO TO 91 101 DO 11 I=2,L0 102 103 ! IF(X(I-1)-X(I))11,95,96 104 IF(X(I-1)-X(I).EQ.0.0D0) GO TO 95 105 IF(X(I-1)-X(I).GT.0.0D0) GO TO 96 106 11 CONTINUE 107 IPV=0 108 ! 109 !*** MAIN LOOP 110 FINT=0.0D0 111 DO 80 K=1,N0 112 UK=U(K) 113 ! 114 !*** ROUTINE TO LOCATE THE DESIRED POINT 115 IF(UK.GE.X(L0)) GO TO 26 116 IF(UK.LT.X(1)) GO TO 25 117 IMN=2 118 IMX=L0 119 21 I=(IMN+IMX)/2 120 IF(UK.GE.X(I)) GO TO 23 121 IMX=I 122 GO TO 24 123 23 IMN=I+1 124 24 IF(IMX.GT.IMN) GO TO 21 125 I=IMX 126 GO TO 30 127 25 I=1 128 GO TO 30 129 26 I=LP1 130 GO TO 30 131 ! 132 !*** CHECK IF I=IPV 133 30 IF(I.EQ.IPV) GO TO 70 134 IPV=I 135 ! 136 !*** ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND TO 137 !*** ESTIMATE THEM IF NECESSARY 138 J=I 139 IF(J.EQ.1) J=2 140 IF(J.EQ.LP1) J=L0 141 X3=X(J-1) 142 Y3=Y(J-1) 143 X4=X(J) 144 Y4=Y(J) 145 A3=X4-X3 146 M3=(Y4-Y3)/A3 147 IF(LM2.EQ.0) GO TO 43 148 IF(J.EQ.2) GO TO 41 149 X2=X(J-2) 150 Y2=Y(J-2) 151 A2=X3-X2 152 M2=(Y3-Y2)/A2 153 IF(J.EQ.L0) GO TO 42 154 41 X5=X(J+1) 155 Y5=Y(J+1) 156 A4=X5-X4 157 M4=(Y5-Y4)/A4 158 IF(J.EQ.2) M2=M3+M3-M4 159 GO TO 45 160 42 M4=M3+M3-M2 161 GO TO 45 162 43 M2=M3 163 45 IF(J.LE.3) GO TO 46 164 A1=X2-X(J-3) 165 M1=(Y2-Y(J-3))/A1 166 GO TO 47 167 46 M1=M2+M2-M3 168 47 IF(J.GE.LM1) GO TO 48 169 A5=X(J+2)-X5 170 M5=(Y(J+2)-Y5)/A5 171 GO TO 50 172 48 M5=M4+M4-M3 173 ! 174 !*** NUMERICAL DIFFERENTIATION 175 50 IF(I.EQ.LP1) GO TO 52 176 W2=ABS(M4-M3) 177 W3=ABS(M2-M1) 178 SW=W2+W3 179 IF(SW.NE.0.0) GO TO 51 180 W2=0.5D0 181 W3=0.5D0 182 SW=1.0D0 183 51 T3=(W2*M2+W3*M3)/SW 184 IF(I.EQ.1) GO TO 54 185 52 W3=ABS(M5-M4) 186 W4=ABS(M3-M2) 187 SW=W3+W4 188 IF(SW.NE.0.0) GO TO 53 189 W3=0.5D0 190 W4=0.5D0 191 SW=1.0D0 192 53 T4=(W3*M3+W4*M4)/SW 193 IF(I.NE.LP1) GO TO 60 194 T3=T4 195 SA=A2+A3 196 T4=0.5D0*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA)) 197 X3=X4 198 Y3=Y4 199 A3=A2 200 M3=M4 201 GO TO 60 202 54 T4=T3 203 SA=A3+A4 204 T3=0.5D0*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA)) 205 X3=X3-A4 206 Y3=Y3-M2*A4 207 A3=A4 208 M3=M2 209 ! 210 !*** COMPUTATION OF THE POLYNOMIAL 211 60 Q2=(2.0D0*(M3-T3)+M3-T4)/A3 212 Q3=(-M3-M3+T3+T4)/(A3*A3) 213 70 DX=UK-P0 214 V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3)) 215 216 IF(IDERIV.EQ.0) GO TO 80 217 DV(K)=Q1+DX*(2.0D0*Q2+DX*3.0D0*Q3) 218 DV2(k)=6.0D0*Q3*DX+2.d0*Q2 219 QQ(1,K)=Q0 220 QQ(2,K)=Q1 221 QQ(3,K)=Q2 222 QQ(4,K)=Q3 223 80 CONTINUE 224 RETURN 225 ! 226 !*** ERROR EXIT 227 90 write(std_out,2090) 228 GO TO 99 229 91 write(std_out,2091) 230 GO TO 99 231 95 write(std_out,2095) 232 GO TO 97 233 96 write(std_out,2096) 234 97 write(std_out,2097)I,X(I) 235 99 write(std_out,2099) L0,N0 236 RETURN 237 ! 238 !*** FORMAT STATEMENTS 239 2089 FORMAT( 'WARNING ERROR IN INTRPL. MAX ALLOWED VALUE OF N0 IS',& 240 & I3,' HERE N0 IS',I3) 241 2090 FORMAT(1X/' N = 1 OR LESS.'/) 242 2091 FORMAT(1X/' N = 0 OR LESS.'/) 243 2095 FORMAT(1X/' IDENTICAL X VALUES.'/) 244 2096 FORMAT(1X/' X VALUES OUT OF SEQUENCE.'/) 245 2097 FORMAT(4X,'I =',I7,10X,6X,'X(I) =',E12.3) 246 2099 FORMAT(4X,'L =',I7,10X,3X,'N =',I7/ & 247 & ' ERROR DETECTED IN ROUTINE INTRPL') 248 ! 249 END