eddy/mountcontrol
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
mountcontrol/cxx/fitpack/fpgrpa.f
Timur A. Fatkhullin 5279d1c41a add FITPACK Fortran library 
start developing of FITPACK C++ bindings
mount_server.cpp: fix compilation error with GCC15
2025-05-05 17:24:21 +03:00

315 lines
9.4 KiB
Fortran
Raw Blame History

recursive subroutine fpgrpa(ifsu,ifsv,ifbu,ifbv,idim,ipar,u,mu,
* v,mv,z,mz,tu,nu,tv,nv,p,c,nc,fp,fpu,fpv,mm,mvnu,spu,spv,
* right,q,au,au1,av,av1,bu,bv,nru,nrv)
implicit none
c ..
c ..scalar arguments..
real*8 p,fp
integer ifsu,ifsv,ifbu,ifbv,idim,mu,mv,mz,nu,nv,nc,mm,mvnu
c ..array arguments..
real*8 u(mu),v(mv),z(mz*idim),tu(nu),tv(nv),c(nc*idim),fpu(nu),
* fpv(nv),spu(mu,4),spv(mv,4),right(mm*idim),q(mvnu),au(nu,5),
* au1(nu,4),av(nv,5),av1(nv,4),bu(nu,5),bv(nv,5)
integer ipar(2),nru(mu),nrv(mv)
c ..local scalars..
real*8 arg,fac,term,one,half,value
integer i,id,ii,it,iz,i1,i2,j,jz,k,k1,k2,l,l1,l2,mvv,k0,muu,
* ncof,nroldu,nroldv,number,nmd,numu,numu1,numv,numv1,nuu,nvv,
* nu4,nu7,nu8,nv4,nv7,nv8, n33
c ..local arrays..
real*8 h(5)
c ..subroutine references..
c fpback,fpbspl,fpdisc,fpbacp,fptrnp,fptrpe
c ..
c let
c | (spu) | | (spv) |
c (au) = | ---------- | (av) = | ---------- |
c | (1/p) (bu) | | (1/p) (bv) |
c
c | z ' 0 |
c q = | ------ |
c | 0 ' 0 |
c
c with c : the (nu-4) x (nv-4) matrix which contains the b-spline
c coefficients.
c z : the mu x mv matrix which contains the function values.
c spu,spv: the mu x (nu-4), resp. mv x (nv-4) observation matrices
c according to the least-squares problems in the u-,resp.
c v-direction.
c bu,bv : the (nu-7) x (nu-4),resp. (nv-7) x (nv-4) matrices
c containing the discontinuity jumps of the derivatives
c of the b-splines in the u-,resp.v-variable at the knots
c the b-spline coefficients of the smoothing spline are then calculated
c as the least-squares solution of the following over-determined linear
c system of equations
c
c (1) (av) c (au)' = q
c
c subject to the constraints
c
c (2) c(nu-3+i,j) = c(i,j), i=1,2,3 ; j=1,2,...,nv-4
c if(ipar(1).ne.0)
c
c (3) c(i,nv-3+j) = c(i,j), j=1,2,3 ; i=1,2,...,nu-4
c if(ipar(2).ne.0)
c
c set constants
one = 1
half = 0.5
c initialization
nu4 = nu-4
nu7 = nu-7
nu8 = nu-8
nv4 = nv-4
nv7 = nv-7
nv8 = nv-8
muu = mu
if(ipar(1).ne.0) muu = mu-1
mvv = mv
if(ipar(2).ne.0) mvv = mv-1
c it depends on the value of the flags ifsu,ifsv,ifbu and ibvand
c on the value of p whether the matrices (spu), (spv), (bu) and (bv)
c still must be determined.
if(ifsu.ne.0) go to 50
c calculate the non-zero elements of the matrix (spu) which is the ob-
c servation matrix according to the least-squares spline approximation
c problem in the u-direction.
l = 4
l1 = 5
number = 0
do 40 it=1,muu
arg = u(it)
10 if(arg.lt.tu(l1) .or. l.eq.nu4) go to 20
l = l1
l1 = l+1
number = number+1
go to 10
20 call fpbspl(tu,nu,3,arg,l,h)
do 30 i=1,4
spu(it,i) = h(i)
30 continue
nru(it) = number
40 continue
ifsu = 1
c calculate the non-zero elements of the matrix (spv) which is the ob-
c servation matrix according to the least-squares spline approximation
c problem in the v-direction.
50 if(ifsv.ne.0) go to 100
l = 4
l1 = 5
number = 0
do 90 it=1,mvv
arg = v(it)
60 if(arg.lt.tv(l1) .or. l.eq.nv4) go to 70
l = l1
l1 = l+1
number = number+1
go to 60
70 call fpbspl(tv,nv,3,arg,l,h)
do 80 i=1,4
spv(it,i) = h(i)
80 continue
nrv(it) = number
90 continue
ifsv = 1
100 if(p.le.0.) go to 150
c calculate the non-zero elements of the matrix (bu).
if(ifbu.ne.0 .or. nu8.eq.0) go to 110
call fpdisc(tu,nu,5,bu,nu)
ifbu = 1
c calculate the non-zero elements of the matrix (bv).
110 if(ifbv.ne.0 .or. nv8.eq.0) go to 150
call fpdisc(tv,nv,5,bv,nv)
ifbv = 1
c substituting (2) and (3) into (1), we obtain the overdetermined
c system
c (4) (avv) (cr) (auu)' = (qq)
c from which the nuu*nvv remaining coefficients
c c(i,j) , i=1,...,nu-4-3*ipar(1) ; j=1,...,nv-4-3*ipar(2) ,
c the elements of (cr), are then determined in the least-squares sense.
c we first determine the matrices (auu) and (qq). then we reduce the
c matrix (auu) to upper triangular form (ru) using givens rotations.
c we apply the same transformations to the rows of matrix qq to obtain
c the (mv) x nuu matrix g.
c we store matrix (ru) into au (and au1 if ipar(1)=1) and g into q.
150 if(ipar(1).ne.0) go to 160
nuu = nu4
call fptrnp(mu,mv,idim,nu,nru,spu,p,bu,z,au,q,right)
go to 180
160 nuu = nu7
call fptrpe(mu,mv,idim,nu,nru,spu,p,bu,z,au,au1,q,right)
c we determine the matrix (avv) and then we reduce this matrix to
c upper triangular form (rv) using givens rotations.
c we apply the same transformations to the columns of matrix
c g to obtain the (nvv) x (nuu) matrix h.
c we store matrix (rv) into av (and av1 if ipar(2)=1) and h into c.
180 if(ipar(2).ne.0) go to 190
nvv = nv4
call fptrnp(mv,nuu,idim,nv,nrv,spv,p,bv,q,av,c,right)
go to 200
190 nvv = nv7
call fptrpe(mv,nuu,idim,nv,nrv,spv,p,bv,q,av,av1,c,right)
c backward substitution to obtain the b-spline coefficients as the
c solution of the linear system (rv) (cr) (ru)' = h.
c first step: solve the system (rv) (c1) = h.
200 ncof = nuu*nvv
k = 1
if(ipar(2).ne.0) go to 240
do 220 ii=1,idim
do 220 i=1,nuu
call fpback(av,c(k),nvv,5,c(k),nv)
k = k+nvv
220 continue
go to 300
240 do 260 ii=1,idim
do 260 i=1,nuu
call fpbacp(av,av1,c(k),nvv,4,c(k),5,nv)
k = k+nvv
260 continue
c second step: solve the system (cr) (ru)' = (c1).
300 if(ipar(1).ne.0) go to 400
do 360 ii=1,idim
k = (ii-1)*ncof
do 360 j=1,nvv
k = k+1
l = k
do 320 i=1,nuu
right(i) = c(l)
l = l+nvv
320 continue
call fpback(au,right,nuu,5,right,nu)
l = k
do 340 i=1,nuu
c(l) = right(i)
l = l+nvv
340 continue
360 continue
go to 500
400 do 460 ii=1,idim
k = (ii-1)*ncof
do 460 j=1,nvv
k = k+1
l = k
do 420 i=1,nuu
right(i) = c(l)
l = l+nvv
420 continue
call fpbacp(au,au1,right,nuu,4,right,5,nu)
l = k
do 440 i=1,nuu
c(l) = right(i)
l = l+nvv
440 continue
460 continue
c calculate from the conditions (2)-(3), the remaining b-spline
c coefficients.
500 if(ipar(2).eq.0) go to 600
i = 0
j = 0
do 560 id=1,idim
do 560 l=1,nuu
ii = i
do 520 k=1,nvv
i = i+1
j = j+1
q(i) = c(j)
520 continue
do 540 k=1,3
ii = ii+1
i = i+1
q(i) = q(ii)
540 continue
560 continue
ncof = nv4*nuu
nmd = ncof*idim
do 580 i=1,nmd
c(i) = q(i)
580 continue
600 if(ipar(1).eq.0) go to 700
i = 0
j = 0
n33 = 3*nv4
do 660 id=1,idim
ii = i
do 620 k=1,ncof
i = i+1
j = j+1
q(i) = c(j)
620 continue
do 640 k=1,n33
ii = ii+1
i = i+1
q(i) = q(ii)
640 continue
660 continue
ncof = nv4*nu4
nmd = ncof*idim
do 680 i=1,nmd
c(i) = q(i)
680 continue
c calculate the quantities
c res(i,j) = (z(i,j) - s(u(i),v(j)))**2 , i=1,2,..,mu;j=1,2,..,mv
c fp = sumi=1,mu(sumj=1,mv(res(i,j)))
c fpu(r) = sum''i(sumj=1,mv(res(i,j))) , r=1,2,...,nu-7
c tu(r+3) <= u(i) <= tu(r+4)
c fpv(r) = sumi=1,mu(sum''j(res(i,j))) , r=1,2,...,nv-7
c tv(r+3) <= v(j) <= tv(r+4)
700 fp = 0.
do 720 i=1,nu
fpu(i) = 0.
720 continue
do 740 i=1,nv
fpv(i) = 0.
740 continue
nroldu = 0
c main loop for the different grid points.
do 860 i1=1,muu
numu = nru(i1)
numu1 = numu+1
nroldv = 0
iz = (i1-1)*mv
do 840 i2=1,mvv
numv = nrv(i2)
numv1 = numv+1
iz = iz+1
c evaluate s(u,v) at the current grid point by making the sum of the
c cross products of the non-zero b-splines at (u,v), multiplied with
c the appropriate b-spline coefficients.
term = 0.
k0 = numu*nv4+numv
jz = iz
do 800 id=1,idim
k1 = k0
value = 0.
do 780 l1=1,4
k2 = k1
fac = spu(i1,l1)
do 760 l2=1,4
k2 = k2+1
value = value+fac*spv(i2,l2)*c(k2)
760 continue
k1 = k1+nv4
780 continue
c calculate the squared residual at the current grid point.
term = term+(z(jz)-value)**2
jz = jz+mz
k0 = k0+ncof
800 continue
c adjust the different parameters.
fp = fp+term
fpu(numu1) = fpu(numu1)+term
fpv(numv1) = fpv(numv1)+term
fac = term*half
if(numv.eq.nroldv) go to 820
fpv(numv1) = fpv(numv1)-fac
fpv(numv) = fpv(numv)+fac
820 nroldv = numv
if(numu.eq.nroldu) go to 840
fpu(numu1) = fpu(numu1)-fac
fpu(numu) = fpu(numu)+fac
840 continue
nroldu = numu
860 continue
return
end
Reference in New Issue View Git Blame Copy Permalink