315 lines
9.4 KiB
Fortran
315 lines
9.4 KiB
Fortran
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
|