TABLE OF CONTENTS
ABINIT/cspint [ Functions ]
NAME
cspint
FUNCTION
CSPINT estimates the integral of a tabulated function.
INPUTS
OUTPUT
NOTES
The routine is given the value of a function F(X) at a set of
nodes XTAB, and estimates
Integral ( A <= X <= B ) F(X) DX
by computing the cubic natural spline S(X) that interpolates
F(X) at the nodes, and then computing
Integral ( A <= X <= B ) S(X) DX
exactly.
Other output from the program includes the definite integral
from X(1) to X(I) of S(X), and the coefficients necessary for
the user to evaluate the spline S(X) at any point.
Modified:
30 October 2000
Reference:
Philip Davis and Philip Rabinowitz,
Methods of Numerical Integration,
Blaisdell Publishing, 1967.
Parameters:
Input, real (dp) FTAB(NTAB), contains the tabulated values of
the function, FTAB(I) = F(XTAB(I)).
Input, real (dp) XTAB(NTAB), contains the points at which the
function was evaluated. The XTAB's must be distinct and
in ascending order.
Input, integer NTAB, the number of entries in FTAB and
XTAB. NTAB must be at least 3.
Input, real (dp) A, lower limit of integration.
Input, real (dp) B, upper limit of integration.
Output, real (dp) Y(3,NTAB), will contain the coefficients
of the interpolating natural spline over each subinterval.
For XTAB(I) <= X <= XTAB(I+1),
S(X) = FTAB(I) + Y(1,I)*(X-XTAB(I))
+ Y(2,I)*(X-XTAB(I))**2
+ Y(3,I)*(X-XTAB(I))**3
Output, real (dp) E(NTAB), E(I) = the definite integral from
XTAB(1) to XTAB(I) of S(X).
Workspace, real (dp) WORK(NTAB).
Output, real (dp) RESULT, the estimated value of the integral.
PARENTS
m_xc_vdw,mrgscr
CHILDREN
SOURCE
81 #if defined HAVE_CONFIG_H 82 #include "config.h" 83 #endif 84 85 #include "abi_common.h" 86 87 subroutine cspint ( ftab, xtab, ntab, a, b, y, e, work, result ) 88 89 use defs_basis 90 use m_profiling_abi 91 use m_errors 92 93 !This section has been created automatically by the script Abilint (TD). 94 !Do not modify the following lines by hand. 95 #undef ABI_FUNC 96 #define ABI_FUNC 'cspint' 97 !End of the abilint section 98 99 implicit none 100 101 integer, intent(in) :: ntab 102 103 real(dp), intent(in) :: a 104 real(dp), intent(in) :: b 105 real(dp), intent(inout) :: e(ntab) 106 real(dp), intent(in) :: ftab(ntab) 107 integer :: i 108 integer :: j 109 real(dp) :: r 110 real(dp), intent(out) :: result 111 real(dp) :: s 112 real(dp) :: term 113 real(dp) :: u 114 real(dp), intent(inout) :: work(ntab) 115 real(dp), intent(in) :: xtab(ntab) 116 real(dp), intent(inout) :: y(3,ntab) 117 118 if ( ntab < 3 ) then 119 write(std_out,'(a)' ) ' ' 120 write(std_out,'(a)' ) 'CSPINT - Fatal error!' 121 write(std_out,'(a,i6)' ) ' NTAB must be at least 3, but input NTAB = ',ntab 122 MSG_ERROR("Aborting now") 123 end if 124 125 do i = 1, ntab-1 126 127 if ( xtab(i+1) <= xtab(i) ) then 128 write(std_out,'(a)' ) ' ' 129 write(std_out,'(a)' ) 'CSPINT - Fatal error!' 130 write(std_out,'(a)' ) ' Nodes not in strict increasing order.' 131 write(std_out,'(a,i6)' ) ' XTAB(I) <= XTAB(I-1) for I=',i 132 write(std_out,'(a,g14.6)' ) ' XTAB(I) = ',xtab(i) 133 write(std_out,'(a,g14.6)' ) ' XTAB(I-1) = ',xtab(i-1) 134 MSG_ERROR("Aborting now") 135 end if 136 137 end do 138 139 s = zero 140 do i = 1, ntab-1 141 r = ( ftab(i+1) - ftab(i) ) / ( xtab(i+1) - xtab(i) ) 142 y(2,i) = r - s 143 s = r 144 end do 145 146 result = zero 147 s = zero 148 r = zero 149 y(2,1) = zero 150 y(2,ntab) = zero 151 152 do i = 2, ntab-1 153 y(2,i) = y(2,i) + r * y(2,i-1) 154 work(i) = two * ( xtab(i-1) - xtab(i+1) ) - r * s 155 s = xtab(i+1) - xtab(i) 156 r = s / work(i) 157 end do 158 159 do j = 2, ntab-1 160 i = ntab+1-j 161 y(2,i) = ( ( xtab(i+1) - xtab(i) ) * y(2,i+1) - y(2,i) ) / work(i) 162 end do 163 164 do i = 1, ntab-1 165 s = xtab(i+1) - xtab(i) 166 r = y(2,i+1) - y(2,i) 167 y(3,i) = r / s 168 y(2,i) = three * y(2,i) 169 y(1,i) = ( ftab(i+1) - ftab(i) ) / s - ( y(2,i) + r ) * s 170 end do 171 172 e(1) = 0.0D+00 173 do i = 1, ntab-1 174 s = xtab(i+1)-xtab(i) 175 term = ((( y(3,i) * quarter * s + y(2,i) * third ) * s & 176 + y(1,i) * half ) * s + ftab(i) ) * s 177 e(i+1) = e(i) + term 178 end do 179 ! 180 ! Determine where the endpoints A and B lie in the mesh of XTAB's. 181 ! 182 r = a 183 u = one 184 185 do j = 1, 2 186 ! 187 ! The endpoint is less than or equal to XTAB(1). 188 ! 189 if ( r <= xtab(1) ) then 190 result = result-u*((r-xtab(1))*y(1,1)*half +ftab(1))*(r-xtab(1)) 191 ! 192 ! The endpoint is greater than or equal to XTAB(NTAB). 193 ! 194 else if ( xtab(ntab) <= r ) then 195 196 result = result -u * ( e(ntab) + ( r - xtab(ntab) ) & 197 * ( ftab(ntab) + half * ( ftab(ntab-1) & 198 + ( xtab(ntab) - xtab(ntab-1) ) * y(1,ntab-1) ) & 199 * ( r - xtab(ntab) ))) 200 ! 201 ! The endpoint is strictly between XTAB(1) and XTAB(NTAB). 202 ! 203 else 204 205 do i = 1, ntab-1 206 207 if ( r <= xtab(i+1) ) then 208 r = r-xtab(i) 209 result = result-u*(e(i)+(((y(3,i)*quarter*r+y(2,i)*third)*r & 210 +y(1,i)*half )*r+ftab(i))*r) 211 go to 120 212 end if 213 214 end do 215 216 end if 217 218 120 continue 219 220 u = -one 221 r = b 222 223 end do 224 225 return 226 227 end subroutine cspint