Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/newnot.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

newnot.f90
View Options

!>
!> @file newnot.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 newnot ( break, coef, l, k, brknew, lnew, coefg )
!*************************************************************************
!
!! NEWNOT returns LNEW+1 knots which are equidistributed on (A,B).
!
! Discussion:
!
! (a,b) = (break(1),break(l+1)) with respect to a certain monotone
! function G related to the K-th root of the K-th derivative of the
! piecewise polynomial function F whose piecewise polynomial
! representation is contained in break, coef, l, k .
!
! method
!
! The K-th derivative of the given piecewise polynomial function F does
! not exist, except perhaps as a linear combination of delta functions.
! Nevertheless, we construct a piecewise constant function H with
! breakpoint sequence BREAK which is approximately proportional
! to abs(d**k(f)).
!
! Specifically, on (break(i), break(i+1)),
!
! abs(jump at break(i) of pc) abs(jump at break(i+1) of pc)
! h=-------------- + ----------------------------
! break(i+1)-break(i-1) break(i+2) - break(i)
!
! with pc the p.constant (k-1)st derivative of f .
! then, the p.linear function g is constructed as
!
! g(x) = integral of h(y)**(1/k) for y from a to x
!
! and its pp coefficients stored in coefg .
!
! then brknew is determined by
!
! brknew(i) = a+g**(-1)((i-1)*step) , i=1,...,lnew+1
!
! where step=g(b)/lnew and (a,b) = (break(1),break(l+1)).
!
! In the event that pc=d**(k-1)(f) is constant in (a,b) and
! therefore h=0 identically, brknew is chosen uniformly spaced.
!
! optional printed output
! coefg.....the pp coefficients of g are printed out if iprint is set
! positive in data statement below.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, break, coef, l, k.....contains the pp-representation of a certain
! function f of order k . specifically,
! d**(k-1)f(x)=coef(k,i) for break(i) <= x < break(i+1)
!
! Input, lnew.....number of intervals into which the interval (a,b) is to be
! sectioned by the new breakpoint sequence brknew .
!
! Output, real ( kind = 8 ) BRKNEW(LNEW+1), the new breakpoint sequence.
!
! Output, coefg.....the coefficient part of the pp-repr. break, coefg, l, 2
! for the monotone p.linear function g with respect to which brknew will
! be equidistributed.
!
implicit none
integer k
integer l
integer lnew
real ( kind = 8 ) break(l+1)
real ( kind = 8 ) brknew(lnew+1)
real ( kind = 8 ) coef(k,l)
real ( kind = 8 ) coefg(2,l)
real ( kind = 8 ) dif
real ( kind = 8 ) difprv
integer i
integer, save :: iprint = 0
integer j
real ( kind = 8 ) oneovk
real ( kind = 8 ) step
real ( kind = 8 ) stepi
brknew(1) = break(1)
brknew(lnew+1) = break(l+1)
!
! If G is constant, BRKNEW is uniform.
!
if ( l <= 1) then
step = (break(l+1)-break(1))/ real ( lnew, kind = 8 )
do i = 2, lnew
brknew(i) = break(1) + real ( i - 1, kind = 8 ) * step
end do
return
end if
!
! Construct the continuous piecewise linear function G.
!
oneovk = 1.0D+00 / real ( k, kind = 8 )
coefg(1,1) = 0.0D+00
difprv = abs(coef(k,2)-coef(k,1))/(break(3)-break(1))
do i = 2, l
dif = abs(coef(k,i)-coef(k,i-1))/(break(i+1)-break(i-1))
coefg(2,i-1) = (dif+difprv)**oneovk
coefg(1,i) = coefg(1,i-1)+coefg(2,i-1)*(break(i)-break(i-1))
difprv = dif
end do
coefg(2,l) = ( 2.0D+00 * difprv )**oneovk
!
! step = g(b)/lnew
!
step=(coefg(1,l)+coefg(2,l)*(break(l+1)-break(l))) / real ( lnew, kind = 8 )
if ( 0 < iprint ) then
write(*,20)step,(i,coefg(1,i),coefg(2,i),i=1,l)
end if
20 format (' step =',e16.7/(i5,2e16.5))
!
! if G is constant, BRKNEW is uniform.
!
if ( step <= 0.0D+00 ) then
step = (break(l+1)-break(1)) / real ( lnew, kind = 8 )
do i = 2, lnew
brknew(i) = break(1) + real ( i - 1, kind = 8 ) * step
end do
return
end if
!
! for i=2,...,lnew, construct brknew(i)=a+g**(-1)(stepi),
! with stepi=(i-1)*step . this requires inversion of the p.lin-
! ear function g . for this, j is found so that
! g(break(j)) <= stepi .le. g(break(j+1))
! and then
! brknew(i) = break(j)+(stepi-g(break(j)))/dg(break(j)) .
! the midpoint is chosen if dg(break(j))=0 .
!
j = 1
do i = 2, lnew
stepi = real ( i - 1, kind = 8 ) * step
do
if ( j == l ) then
exit
end if
if ( stepi <= coefg(1,j+1) ) then
exit
end if
j = j + 1
end do
if ( coefg(2,j) /= 0.0D+00 ) then
brknew(i) = break(j)+(stepi-coefg(1,j))/coefg(2,j)
else
brknew(i) = ( break(j) + break(j+1) ) / 2.0D+00
end if
end do
return
end
c4science · Help