Timur A. Fatkhullin
5279d1c41a
start developing of FITPACK C++ bindings mount_server.cpp: fix compilation error with GCC15
132 lines
3.7 KiB
Fortran
132 lines
3.7 KiB
Fortran
recursive subroutine fpintb(t,n,bint,nk1,x,y)
|
|
implicit none
|
|
c subroutine fpintb calculates integrals of the normalized b-splines
|
|
c nj,k+1(x) of degree k, defined on the set of knots t(j),j=1,2,...n.
|
|
c it makes use of the formulae of gaffney for the calculation of
|
|
c indefinite integrals of b-splines.
|
|
c
|
|
c calling sequence:
|
|
c call fpintb(t,n,bint,nk1,x,y)
|
|
c
|
|
c input parameters:
|
|
c t : real array,length n, containing the position of the knots.
|
|
c n : integer value, giving the number of knots.
|
|
c nk1 : integer value, giving the number of b-splines of degree k,
|
|
c defined on the set of knots ,i.e. nk1 = n-k-1.
|
|
c x,y : real values, containing the end points of the integration
|
|
c interval.
|
|
c output parameter:
|
|
c bint : array,length nk1, containing the integrals of the b-splines.
|
|
c ..
|
|
c ..scalars arguments..
|
|
integer n,nk1
|
|
real*8 x,y
|
|
c ..array arguments..
|
|
real*8 t(n),bint(nk1)
|
|
c ..local scalars..
|
|
integer i,ia,ib,it,j,j1,k,k1,l,li,lj,lk,l0,min
|
|
real*8 a,ak,arg,b,f,one
|
|
c ..local arrays..
|
|
real*8 aint(6),h(6),h1(6)
|
|
c initialization.
|
|
one = 0.1d+01
|
|
k1 = n-nk1
|
|
ak = k1
|
|
k = k1-1
|
|
do 10 i=1,nk1
|
|
bint(i) = 0.0d0
|
|
10 continue
|
|
c the integration limits are arranged in increasing order.
|
|
a = x
|
|
b = y
|
|
min = 0
|
|
if (a.lt.b) go to 30
|
|
if (a.eq.b) go to 160
|
|
go to 20
|
|
20 a = y
|
|
b = x
|
|
min = 1
|
|
30 if(a.lt.t(k1)) a = t(k1)
|
|
if(b.gt.t(nk1+1)) b = t(nk1+1)
|
|
if(a.gt.b) go to 160
|
|
c using the expression of gaffney for the indefinite integral of a
|
|
c b-spline we find that
|
|
c bint(j) = (t(j+k+1)-t(j))*(res(j,b)-res(j,a))/(k+1)
|
|
c where for t(l) <= x < t(l+1)
|
|
c res(j,x) = 0, j=1,2,...,l-k-1
|
|
c = 1, j=l+1,l+2,...,nk1
|
|
c = aint(j+k-l+1), j=l-k,l-k+1,...,l
|
|
c = sumi((x-t(j+i))*nj+i,k+1-i(x)/(t(j+k+1)-t(j+i)))
|
|
c i=0,1,...,k
|
|
l = k1
|
|
l0 = l+1
|
|
c set arg = a.
|
|
arg = a
|
|
do 90 it=1,2
|
|
c search for the knot interval t(l) <= arg < t(l+1).
|
|
40 if(arg.lt.t(l0) .or. l.eq.nk1) go to 50
|
|
l = l0
|
|
l0 = l+1
|
|
go to 40
|
|
c calculation of aint(j), j=1,2,...,k+1.
|
|
c initialization.
|
|
50 do 55 j=1,k1
|
|
aint(j) = 0.0d0
|
|
55 continue
|
|
aint(1) = (arg-t(l))/(t(l+1)-t(l))
|
|
h1(1) = one
|
|
do 70 j=1,k
|
|
c evaluation of the non-zero b-splines of degree j at arg,i.e.
|
|
c h(i+1) = nl-j+i,j(arg), i=0,1,...,j.
|
|
h(1) = 0.0d0
|
|
do 60 i=1,j
|
|
li = l+i
|
|
lj = li-j
|
|
f = h1(i)/(t(li)-t(lj))
|
|
h(i) = h(i)+f*(t(li)-arg)
|
|
h(i+1) = f*(arg-t(lj))
|
|
60 continue
|
|
c updating of the integrals aint.
|
|
j1 = j+1
|
|
do 70 i=1,j1
|
|
li = l+i
|
|
lj = li-j1
|
|
aint(i) = aint(i)+h(i)*(arg-t(lj))/(t(li)-t(lj))
|
|
h1(i) = h(i)
|
|
70 continue
|
|
if(it.eq.2) go to 100
|
|
c updating of the integrals bint
|
|
lk = l-k
|
|
ia = lk
|
|
do 80 i=1,k1
|
|
bint(lk) = -aint(i)
|
|
lk = lk+1
|
|
80 continue
|
|
c set arg = b.
|
|
arg = b
|
|
90 continue
|
|
c updating of the integrals bint.
|
|
100 lk = l-k
|
|
ib = lk-1
|
|
do 110 i=1,k1
|
|
bint(lk) = bint(lk)+aint(i)
|
|
lk = lk+1
|
|
110 continue
|
|
if(ib.lt.ia) go to 130
|
|
do 120 i=ia,ib
|
|
bint(i) = bint(i)+one
|
|
120 continue
|
|
c the scaling factors are taken into account.
|
|
130 f = one/ak
|
|
do 140 i=1,nk1
|
|
j = i+k1
|
|
bint(i) = bint(i)*(t(j)-t(i))*f
|
|
140 continue
|
|
c the order of the integration limits is taken into account.
|
|
if(min.eq.0) go to 160
|
|
do 150 i=1,nk1
|
|
bint(i) = -bint(i)
|
|
150 continue
|
|
160 return
|
|
end
|