Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/tautsp.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

tautsp.f90
View Options

!>
!> @file tautsp.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine tautsp ( tau, gtau, ntau, gamma, s, break, coef, l, k, iflag )
!*************************************************************************
!
!! TAUTSP constructs a cubic spline interpolant to given data.
!
! Discussion:
!
! If 0 < GAMMA, additional knots are introduced where needed to
! make the interpolant more flexible locally. This avoids extraneous
! inflection points typical of cubic spline interpolation at knots to
! rapidly changing data.
!
! Method:
!
! On the I-th interval, (TAU(I), TAU(I+1)), the interpolant is of the
! form:
!
! (*) f(u(x))=a+b*u + c*h(u,z) + d*h(1-u,1-z) ,
!
! with
!
! U = U(X) = ( X - TAU(I) ) / DTAU(I).
!
! Here,
! z=z(i) = addg(i+1)/(addg(i)+addg(i+1))
! (= .5, in case the denominator vanishes). with
! addg(j)=abs(ddg(j)), ddg(j) = dg(j+1)-dg(j),
! dg(j)=divdif ( j) = (gtau(j+1)-gtau(j))/dtau(j)
! and
! h(u,z)=alpha*u**3+(1-alpha)*(max(((u-zeta)/(1-zeta)),0)**3
! with
! alpha(z)=(1-gamma/3)/zeta
! zeta(z)=1-gamma*min((1 - z), 1/3)
! thus, for 1/3 <= z .le. 2/3, f is just a cubic polynomial on
! the interval i. otherwise, it has one additional knot, at
! tau(i)+zeta*dtau(i) .
! as z approaches 1, h(.,z) has an increasingly sharp bend near 1,
! thus allowing f to turn rapidly near the additional knot.
! in terms of f(j)=gtau(j) and
! fsecnd(j)= second derivative of f at tau(j),
! the coefficients for (*) are given as
! a=f(i)-d
! b=(f(i+1)-f(i)) - (c - d)
! c=fsecnd(i+1)*dtau(i)**2/hsecnd(1,z)
! d=fsecnd(i)*dtau(i)**2/hsecnd(1,1-z)
! hence can be computed once fsecnd(i),i=1,...,ntau, is fixed.
!
! F is automatically continuous and has a continuous second derivative
! (except when z=0 or 1 for some i). we determine fscnd(.) from
! the requirement that also the first derivative of F be continuous.
!
! In addition, we require that the third derivative be continuous
! across TAU(2) and across TAU(NTAU-1). This leads to a strictly
! diagonally dominant tridiagonal linear system for the fsecnd(i)
! which we solve by Gauss elimination without pivoting.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) TAU(NTAU), the sequence of data points.
! TAU must be strictly increasing.
!
! Input, real ( kind = 8 ) GTAU(NTAU), the corresponding sequence of
! function values.
!
! Input, integer NTAU, the number of data points. NTAU must be at least 4.
!
! Input, gamma indicates whether additional flexibility is desired.
! =0., no additional knots
! in (0.,3.), under certain conditions on the given data at
! points i-1, i, i+1, and i+2, a knot is added in the
! i-th interval, i=2,...,ntau-2. see description of meth-
! od below. the interpolant gets rounded with increasing
! gamma. a value of 2.5 for gamma is typical.
! in (3.,6.), same , except that knots might also be added in
! intervals in which an inflection point would be permit-
! ted. a value of 5.5 for gamma is typical.
!
! Output, break, coef, l, k give the pp-representation of the interpolant.
! specifically, for break(i) <= x .le. break(i+1), the
! interpolant has the form
! f(x)=coef(1,i) +dx(coef(2,i) +(dx/2)(coef(3,i) +(dx/3)coef(4,i)))
! with dx=x-break(i) and i=1,...,l .
!
! Output, iflag =1, ok
! =2, input was incorrect. a printout specifying the mistake
! was made.
! workspace
!
! Output, s is required, of size (ntau,6). the individual columns of this
! array contain the following quantities mentioned in the write-
! up and below.
! s(.,1)=dtau = tau(.+1)-tau
! s(.,2)=diag = diagonal in linear system
! s(.,3)=u = upper diagonal in linear system
! s(.,4)=r = right side for linear system (initially)
! = fsecnd = solution of linear system , namely the second
! derivatives of interpolant at tau
! s(.,5)=z = indicator of additional knots
! s(.,6)=1/hsecnd(1,x) with x = z or = 1-z. see below.
!
implicit none
integer ntau
real ( kind = 8 ) alph
real ( kind = 8 ) alpha
real ( kind = 8 ) break(*)
real ( kind = 8 ) c
real ( kind = 8 ) coef(4,*)
real ( kind = 8 ) d
real ( kind = 8 ) del
real ( kind = 8 ) denom
real ( kind = 8 ) divdif
real ( kind = 8 ) entry
real ( kind = 8 ) entry3
real ( kind = 8 ) factor
real ( kind = 8 ) factr2
real ( kind = 8 ) gam
real ( kind = 8 ) gamma
real ( kind = 8 ) gtau(ntau)
integer i
integer iflag
integer k
integer l
integer method
real ( kind = 8 ) onemg3
real ( kind = 8 ) onemzt
real ( kind = 8 ) ratio
real ( kind = 8 ) s(ntau,6)
real ( kind = 8 ) sixth
real ( kind = 8 ) tau(ntau)
real ( kind = 8 ) temp
real ( kind = 8 ) x
real ( kind = 8 ) z
real ( kind = 8 ) zeta
real ( kind = 8 ) zt2
alph(x) = min ( 1.0D+00, onemg3 / x )
!
! There must be at least 4 interpolation points.
!
if ( ntau < 4 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TAUTSP - Fatal error!'
write ( *, '(a)' ) ' Input NTAU must be at least 4.'
write ( *, '(a,i6)' ) ' NTAU = ', ntau
iflag = 2
stop
end if
!
! Construct delta tau and first and second (divided) differences of data.
!
do i = 1, ntau-1
s(i,1) = tau(i+1)-tau(i)
if ( s(i,1) <= 0.0D+00 ) then
write(*,30)i,tau(i),tau(i+1)
30 format (' point ',i3,' and the next',2e15.6,' are disordered')
iflag=2
return
end if
s(i+1,4) = ( gtau(i+1) - gtau(i) ) / s(i,1)
end do
do i = 2, ntau-1
s(i,4) = s(i+1,4)-s(i,4)
end do
!
! Construct system of equations for second derivatives at tau. at each
! interior data point, there is one continuity equation, at the first
! and the last interior data point there is an additional one for a
! total of NTAU equations in ntau unknowns.
!
i = 2
s(2,2) = s(1,1) / 3.0D+00
sixth = 1.0D+00 / 6.0D+00
method = 2
gam = gamma
if ( gam <= 0.0D+00 ) then
method=1
end if
if ( 3.0D+00 < gam ) then
method = 3
gam = gam - 3.0D+00
end if
onemg3 = 1.0D+00 - gam / 3.0D+00
!
! loop over i
!
70 continue
!
! Construct z(i) and zeta(i)
!
z = 0.5D+00
if ( method == 1) then
go to 100
end if
if ( method == 3) then
go to 90
end if
if ( s(i,4)*s(i+1,4) < 0.0D+00 ) then
go to 100
end if
90 continue
temp = abs ( s(i+1,4) )
denom = abs ( s(i,4) ) +temp
if ( denom /= 0.0D+00 ) then
z = temp/denom
if ( abs ( z - 0.5D+00 ) <= sixth ) then
z=0.5D+00
end if
end if
100 continue
s(i,5) = z
!
! Set up part of the i-th equation which depends on the i-th interval.
!
if ( z < 0.5D+00 ) then
zeta = gam*z
onemzt = 1.0D+00 - zeta
zt2 = zeta**2
alpha = alph(onemzt)
factor = zeta/(alpha*(zt2-1.0D+00 ) + 1.0D+00 )
s(i,6) = zeta*factor / 6.0D+00
s(i,2) = s(i,2) + s(i,1) &
* ( ( 1.0D+00 - alpha * onemzt ) * factor / 2.0D+00-s(i,6))
!
! If z=0 and the previous z = 1, then d(i) = 0. since then
! also u(i-1)=l(i+1) = 0, its value does not matter. reset
! d(i)=1 to insure nonzero pivot in elimination.
!
if ( s(i,2) <= 0.0D+00 ) then
s(i,2) = 1.0D+00
end if
s(i,3)=s(i,1) / 6.0D+00
else if ( z - 0.5D+00 == 0.0D+00 ) then
s(i,2)=s(i,2)+s(i,1) / 3.0D+00
s(i,3)=s(i,1) / 6.0D+00
else if ( 0.0D+00 < z - 0.5D+00 ) then
onemzt = gam*(1.0D+00 - z)
zeta = 1.0D+00 - onemzt
alpha = alph(zeta)
factor = onemzt/(1.0D+00 - alpha * zeta * ( 1.0D+00 + onemzt ) )
s(i,6) = onemzt*factor / 6.0D+00
s(i,2) = s(i,2)+s(i,1) / 3.0D+00
s(i,3) = s(i,6)*s(i,1)
end if
if ( 2 < i ) then
go to 190
end if
s(1,5) = 0.5D+00
!
! The first two equations enforce continuity of the first and of
! the third derivative across tau(2).
!
s(1,2)=s(1,1) / 6.0D+00
s(1,3)=s(2,2)
entry3=s(2,3)
if ( z-0.5D+00) 150, 160, 170
150 continue
factr2 = zeta * ( alpha * ( zt2 - 1.0D+00 ) + 1.0D+00 ) &
/ ( alpha * ( zeta * zt2 - 1.0D+00 ) + 1.0D+00 )
ratio=factr2*s(2,1)/s(1,2)
s(2,2)=factr2*s(2,1)+s(1,1)
s(2,3)=-factr2*s(1,1)
go to 180
160 continue
ratio=s(2,1)/s(1,2)
s(2,2)=s(2,1)+s(1,1)
s(2,3)=-s(1,1)
go to 180
170 continue
ratio=s(2,1)/s(1,2)
s(2,2)=s(2,1)+s(1,1)
s(2,3)=-s(1,1)*6.0D+00 * alpha * s(2,6)
!
! At this point, the first two equations read
! diag(1)*x1+u(1)*x2 + entry3*x3=r(2)
! -ratio*diag(1)*x1+diag(2)*x2 + u(2)*x3=0.0
! Eliminate first unknown from second equation
!
180 continue
s(2,2)=ratio*s(1,3)+s(2,2)
s(2,3)=ratio*entry3+s(2,3)
s(1,4)=s(2,4)
s(2,4)=ratio*s(1,4)
go to 200
190 continue
!
! The i-th equation enforces continuity of the first derivative
! across tau(i). it has been set up in statements 35 up to 40
! and 21 up to 25 and reads now
! -ratio*diag(i-1)*xi-1+diag(i)*xi + u(i)*xi+1=r(i) .
! eliminate (i-1)st unknown from this equation
!
s(i,2)=ratio*s(i-1,3)+s(i,2)
s(i,4)=ratio*s(i-1,4)+s(i,4)
!
! Set up the part of the next equation which depends on the
! i-th interval.
!
200 continue
if ( z- 0.5D+00 ) 210, 220, 230
210 continue
ratio = -s(i,6) * s(i,1) / s(i,2)
s(i+1,2)=s(i,1) / 3.0D+00
go to 240
220 continue
ratio=-(s(i,1) / 6.0D+00 ) / s(i,2)
s(i+1,2)=s(i,1) / 3.0D+00
go to 240
230 continue
ratio=-( s(i,1) / 6.0D+00 )/s(i,2)
s(i+1,2)=s(i,1)*((1.0D+00-zeta*alpha) * factor / 2.0D+00 - s(i,6) )
!
! end of i loop
!
240 continue
i=i+1
if ( i < ntau-1) then
go to 70
end if
s(i,5) = 0.5D+00
!
! The last two equations enforce continuity of third derivative and
! of first derivative across tau(ntau-1).
!
entry=ratio*s(i-1,3)+s(i,2)+s(i,1)/3.0D+00
s(i+1,2)=s(i,1)/6.0D+00
s(i+1,4)=ratio*s(i-1,4)+s(i,4)
if ( z- 0.5D+00 ) 250, 260, 270
250 continue
ratio = s(i,1) * 6.0D+00 * s(i-1,6) * alpha / s(i-1,2)
s(i,2)=ratio*s(i-1,3)+s(i,1)+s(i-1,1)
s(i,3)=-s(i-1,1)
go to 280
260 continue
ratio=s(i,1)/s(i-1,2)
s(i,2)=ratio*s(i-1,3)+s(i,1)+s(i-1,1)
s(i,3)=-s(i-1,1)
go to 280
270 continue
factr2=onemzt*(alpha*(onemzt**2-1.0D+00)+1.0D+00) &
/(alpha*(onemzt**3-1.0D+00)+1.0D+00)
ratio = factr2*s(i,1) / s(i-1,2)
s(i,2)=ratio*s(i-1,3)+factr2*s(i-1,1)+s(i,1)
s(i,3)=-factr2*s(i-1,1)
!
! At this point, the last two equations read:
!
! diag(i)*xi+ u(i)*xi+1=r(i)
! -ratio*diag(i)*xi+diag(i+1)*xi+1=r(i+1)
!
! Eliminate XI from the last equation.
!
280 continue
s(i,4)=ratio*s(i-1,4)
ratio=-entry/s(i,2)
s(i+1,2)=ratio*s(i,3)+s(i+1,2)
s(i+1,4)=ratio*s(i,4)+s(i+1,4)
!
! Back substitution.
!
s(ntau,4) = s(ntau,4) / s(ntau,2)
290 continue
s(i,4)=(s(i,4)-s(i,3)*s(i+1,4))/s(i,2)
i=i-1
if ( 1 < i ) then
go to 290
end if
s(1,4)=(s(1,4)-s(1,3)*s(2,4)-entry3*s(3,4))/s(1,2)
!
! Construct polynomial pieces.
!
break(1)=tau(1)
l=1
do i=1, ntau-1
coef(1,l)=gtau(i)
coef(3,l)=s(i,4)
divdif=(gtau(i+1)-gtau(i))/s(i,1)
z=s(i,5)
if ( z- 0.5D+00 ) 300, 310, 320
300 continue
if ( z == 0.0D+00 ) go to 330
zeta=gam*z
onemzt=1.0D+00-zeta
c=s(i+1,4) / 6.0D+00
d=s(i,4)*s(i,6)
l=l+1
del=zeta*s(i,1)
break(l)=tau(i)+del
zt2=zeta**2
alpha=alph(onemzt)
factor=onemzt**2*alpha
coef(1,l)=gtau(i)+divdif*del+s(i,1)**2*(d*onemzt*(factor-1.0D+00) &
+c*zeta*(zt2-1.0D+00))
coef(2,l)=divdif+s(i,1)*(d*(1.0D+00-3.0D+00*factor)+c*(3.0D+00*zt2-1.0D+00))
coef(3,l)=6.0D+00*(d*alpha*onemzt+c*zeta)
coef(4,l)=6.0D+00*(c-d*alpha)/s(i,1)
coef(4,l-1)=coef(4,l)-6.0D+00*d*(1.0D+00-alpha)/(del*zt2)
coef(2,l-1)=coef(2,l)-del*(coef(3,l)-(del/2.0D+00)*coef(4,l-1))
go to 340
310 continue
coef(2,l) = divdif - s(i,1) * ( 2.0D+00 * s(i,4) + s(i+1,4) ) / 6.0D+00
coef(4,l)=(s(i+1,4)-s(i,4))/s(i,1)
go to 340
320 continue
onemzt=gam*(1.0D+00-z)
if ( onemzt == 0.0D+00 ) then
go to 330
end if
zeta = 1.0D+00 - onemzt
alpha=alph(zeta)
c=s(i+1,4)*s(i,6)
d=s(i,4)/6.0D+00
del=zeta*s(i,1)
break(l+1)=tau(i)+del
coef(2,l)=divdif-s(i,1)*(2.0D+00*d+c)
coef(4,l)=6.0D+00*(c*alpha-d)/s(i,1)
l=l+1
coef(4,l)=coef(4,l-1)+6.0D+00*(1.0D+00-alpha)*c/(s(i,1)*onemzt**3)
coef(3,l)=coef(3,l-1)+del*coef(4,l-1)
coef(2,l)=coef(2,l-1)+del*(coef(3,l-1)+(del/2.0D+00)*coef(4,l-1))
coef(1,l)=coef(1,l-1)+del*(coef(2,l-1)+(del/2.0D+00)*(coef(3,l-1) &
+(del/3.0D+00)*coef(4,l-1)))
go to 340
330 continue
coef(2,l) = divdif
coef(3,l) = 0D+00
coef(4,l) = 0.0D+00
340 continue
l = l + 1
break(l) = tau(i+1)
end do
l = l - 1
k = 4
iflag = 1
return
end
c4science · Help