659 lines
19 KiB
Fortran
659 lines
19 KiB
Fortran
recursive subroutine fpgrsp(ifsu,ifsv,ifbu,ifbv,iback,u,mu,v,
|
|
* mv,r,mr,dr,iop0,iop1,tu,nu,tv,nv,p,c,nc,sq,fp,fpu,fpv,mm,
|
|
* mvnu,spu,spv,right,q,au,av1,av2,bu,bv,a0,a1,b0,b1,c0,c1,
|
|
* cosi,nru,nrv)
|
|
implicit none
|
|
c ..
|
|
c ..scalar arguments..
|
|
real*8 p,sq,fp
|
|
integer ifsu,ifsv,ifbu,ifbv,iback,mu,mv,mr,iop0,iop1,nu,nv,nc,
|
|
* mm,mvnu
|
|
c ..array arguments..
|
|
real*8 u(mu),v(mv),r(mr),dr(6),tu(nu),tv(nv),c(nc),fpu(nu),fpv(nv)
|
|
*,
|
|
* spu(mu,4),spv(mv,4),right(mm),q(mvnu),au(nu,5),av1(nv,6),c0(nv),
|
|
* av2(nv,4),a0(2,mv),b0(2,nv),cosi(2,nv),bu(nu,5),bv(nv,5),c1(nv),
|
|
* a1(2,mv),b1(2,nv)
|
|
integer nru(mu),nrv(mv)
|
|
c ..local scalars..
|
|
real*8 arg,co,dr01,dr02,dr03,dr11,dr12,dr13,fac,fac0,fac1,pinv,piv
|
|
*,
|
|
* si,term,one,three,half
|
|
integer i,ic,ii,ij,ik,iq,irot,it,ir,i0,i1,i2,i3,j,jj,jk,jper,
|
|
* j0,j1,k,k1,k2,l,l0,l1,l2,mvv,ncof,nrold,nroldu,nroldv,number,
|
|
* numu,numu1,numv,numv1,nuu,nu4,nu7,nu8,nu9,nv11,nv4,nv7,nv8,n1
|
|
c ..local arrays..
|
|
real*8 h(5),h1(5),h2(4)
|
|
c ..function references..
|
|
integer min0
|
|
real*8 cos,sin
|
|
c ..subroutine references..
|
|
c fpback,fpbspl,fpgivs,fpcyt1,fpcyt2,fpdisc,fpbacp,fprota
|
|
c ..
|
|
c let
|
|
c | (spu) | | (spv) |
|
|
c (au) = | -------------- | (av) = | -------------- |
|
|
c | sqrt(1/p) (bu) | | sqrt(1/p) (bv) |
|
|
c
|
|
c | r ' 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 r : 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(i,nv-3+j) = c(i,j), j=1,2,3 ; i=1,2,...,nu-4
|
|
c
|
|
c (3) if iop0 = 0 c(1,j) = dr(1)
|
|
c iop0 = 1 c(1,j) = dr(1)
|
|
c c(2,j) = dr(1)+(dr(2)*cosi(1,j)+dr(3)*cosi(2,j))*
|
|
c tu(5)/3. = c0(j) , j=1,2,...nv-4
|
|
c
|
|
c (4) if iop1 = 0 c(nu-4,j) = dr(4)
|
|
c iop1 = 1 c(nu-4,j) = dr(4)
|
|
c c(nu-5,j) = dr(4)+(dr(5)*cosi(1,j)+dr(6)*cosi(2,j))
|
|
c *(tu(nu-4)-tu(nu-3))/3. = c1(j)
|
|
c
|
|
c set constants
|
|
one = 1
|
|
three = 3
|
|
half = 0.5
|
|
c initialization
|
|
nu4 = nu-4
|
|
nu7 = nu-7
|
|
nu8 = nu-8
|
|
nu9 = nu-9
|
|
nv4 = nv-4
|
|
nv7 = nv-7
|
|
nv8 = nv-8
|
|
nv11 = nv-11
|
|
nuu = nu4-iop0-iop1-2
|
|
if(p.gt.0.) pinv = one/p
|
|
c it depends on the value of the flags ifsu,ifsv,ifbu,ifbv,iop0,iop1
|
|
c and on the value of p whether the matrices (spu), (spv), (bu), (bv),
|
|
c (cosi) still must be determined.
|
|
if(ifsu.ne.0) go to 30
|
|
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 25 it=1,mu
|
|
arg = u(it)
|
|
10 if(arg.lt.tu(l1) .or. l.eq.nu4) go to 15
|
|
l = l1
|
|
l1 = l+1
|
|
number = number+1
|
|
go to 10
|
|
15 call fpbspl(tu,nu,3,arg,l,h)
|
|
do 20 i=1,4
|
|
spu(it,i) = h(i)
|
|
20 continue
|
|
nru(it) = number
|
|
25 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.
|
|
30 if(ifsv.ne.0) go to 85
|
|
l = 4
|
|
l1 = 5
|
|
number = 0
|
|
do 50 it=1,mv
|
|
arg = v(it)
|
|
35 if(arg.lt.tv(l1) .or. l.eq.nv4) go to 40
|
|
l = l1
|
|
l1 = l+1
|
|
number = number+1
|
|
go to 35
|
|
40 call fpbspl(tv,nv,3,arg,l,h)
|
|
do 45 i=1,4
|
|
spv(it,i) = h(i)
|
|
45 continue
|
|
nrv(it) = number
|
|
50 continue
|
|
ifsv = 1
|
|
if(iop0.eq.0 .and. iop1.eq.0) go to 85
|
|
c calculate the coefficients of the interpolating splines for cos(v)
|
|
c and sin(v).
|
|
do 55 i=1,nv4
|
|
cosi(1,i) = 0.
|
|
cosi(2,i) = 0.
|
|
55 continue
|
|
if(nv7.lt.4) go to 85
|
|
do 65 i=1,nv7
|
|
l = i+3
|
|
arg = tv(l)
|
|
call fpbspl(tv,nv,3,arg,l,h)
|
|
do 60 j=1,3
|
|
av1(i,j) = h(j)
|
|
60 continue
|
|
cosi(1,i) = cos(arg)
|
|
cosi(2,i) = sin(arg)
|
|
65 continue
|
|
call fpcyt1(av1,nv7,nv)
|
|
do 80 j=1,2
|
|
do 70 i=1,nv7
|
|
right(i) = cosi(j,i)
|
|
70 continue
|
|
call fpcyt2(av1,nv7,right,right,nv)
|
|
do 75 i=1,nv7
|
|
cosi(j,i+1) = right(i)
|
|
75 continue
|
|
cosi(j,1) = cosi(j,nv7+1)
|
|
cosi(j,nv7+2) = cosi(j,2)
|
|
cosi(j,nv4) = cosi(j,3)
|
|
80 continue
|
|
85 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 90
|
|
call fpdisc(tu,nu,5,bu,nu)
|
|
ifbu = 1
|
|
c calculate the non-zero elements of the matrix (bv).
|
|
90 if(ifbv.ne.0 .or. nv8.eq.0) go to 150
|
|
call fpdisc(tv,nv,5,bv,nv)
|
|
ifbv = 1
|
|
c substituting (2),(3) and (4) into (1), we obtain the overdetermined
|
|
c system
|
|
c (5) (avv) (cc) (auu)' = (qq)
|
|
c from which the nuu*nv7 remaining coefficients
|
|
c c(i,j) , i=2+iop0,3+iop0,...,nu-5-iop1,j=1,2,...,nv-7.
|
|
c the elements of (cc), are then determined in the least-squares sense.
|
|
c simultaneously, we compute the resulting sum of squared residuals sq.
|
|
150 dr01 = dr(1)
|
|
dr11 = dr(4)
|
|
do 155 i=1,mv
|
|
a0(1,i) = dr01
|
|
a1(1,i) = dr11
|
|
155 continue
|
|
if(nv8.eq.0 .or. p.le.0.) go to 165
|
|
do 160 i=1,nv8
|
|
b0(1,i) = 0.
|
|
b1(1,i) = 0.
|
|
160 continue
|
|
165 mvv = mv
|
|
if(iop0.eq.0) go to 195
|
|
fac = (tu(5)-tu(4))/three
|
|
dr02 = dr(2)*fac
|
|
dr03 = dr(3)*fac
|
|
do 170 i=1,nv4
|
|
c0(i) = dr01+dr02*cosi(1,i)+dr03*cosi(2,i)
|
|
170 continue
|
|
do 180 i=1,mv
|
|
number = nrv(i)
|
|
fac = 0.
|
|
do 175 j=1,4
|
|
number = number+1
|
|
fac = fac+c0(number)*spv(i,j)
|
|
175 continue
|
|
a0(2,i) = fac
|
|
180 continue
|
|
if(nv8.eq.0 .or. p.le.0.) go to 195
|
|
do 190 i=1,nv8
|
|
number = i
|
|
fac = 0.
|
|
do 185 j=1,5
|
|
fac = fac+c0(number)*bv(i,j)
|
|
number = number+1
|
|
185 continue
|
|
b0(2,i) = fac*pinv
|
|
190 continue
|
|
mvv = mv+nv8
|
|
195 if(iop1.eq.0) go to 225
|
|
fac = (tu(nu4)-tu(nu4+1))/three
|
|
dr12 = dr(5)*fac
|
|
dr13 = dr(6)*fac
|
|
do 200 i=1,nv4
|
|
c1(i) = dr11+dr12*cosi(1,i)+dr13*cosi(2,i)
|
|
200 continue
|
|
do 210 i=1,mv
|
|
number = nrv(i)
|
|
fac = 0.
|
|
do 205 j=1,4
|
|
number = number+1
|
|
fac = fac+c1(number)*spv(i,j)
|
|
205 continue
|
|
a1(2,i) = fac
|
|
210 continue
|
|
if(nv8.eq.0 .or. p.le.0.) go to 225
|
|
do 220 i=1,nv8
|
|
number = i
|
|
fac = 0.
|
|
do 215 j=1,5
|
|
fac = fac+c1(number)*bv(i,j)
|
|
number = number+1
|
|
215 continue
|
|
b1(2,i) = fac*pinv
|
|
220 continue
|
|
mvv = mv+nv8
|
|
c we first determine the matrices (auu) and (qq). then we reduce the
|
|
c matrix (auu) to an unit upper triangular form (ru) using givens
|
|
c rotations without square roots. we apply the same transformations to
|
|
c the rows of matrix qq to obtain the mv x nuu matrix g.
|
|
c we store matrix (ru) into au and g into q.
|
|
225 l = mvv*nuu
|
|
c initialization.
|
|
sq = 0.
|
|
if(l.eq.0) go to 245
|
|
do 230 i=1,l
|
|
q(i) = 0.
|
|
230 continue
|
|
do 240 i=1,nuu
|
|
do 240 j=1,5
|
|
au(i,j) = 0.
|
|
240 continue
|
|
l = 0
|
|
245 nrold = 0
|
|
n1 = nrold+1
|
|
do 420 it=1,mu
|
|
number = nru(it)
|
|
c find the appropriate column of q.
|
|
250 do 260 j=1,mvv
|
|
right(j) = 0.
|
|
260 continue
|
|
if(nrold.eq.number) go to 280
|
|
if(p.le.0.) go to 410
|
|
c fetch a new row of matrix (bu).
|
|
do 270 j=1,5
|
|
h(j) = bu(n1,j)*pinv
|
|
270 continue
|
|
i0 = 1
|
|
i1 = 5
|
|
go to 310
|
|
c fetch a new row of matrix (spu).
|
|
280 do 290 j=1,4
|
|
h(j) = spu(it,j)
|
|
290 continue
|
|
c find the appropriate column of q.
|
|
do 300 j=1,mv
|
|
l = l+1
|
|
right(j) = r(l)
|
|
300 continue
|
|
i0 = 1
|
|
i1 = 4
|
|
310 j0 = n1
|
|
j1 = nu7-number
|
|
c take into account that we eliminate the constraints (3)
|
|
315 if(j0-1.gt.iop0) go to 335
|
|
fac0 = h(i0)
|
|
do 320 j=1,mv
|
|
right(j) = right(j)-fac0*a0(j0,j)
|
|
320 continue
|
|
if(mv.eq.mvv) go to 330
|
|
j = mv
|
|
do 325 jj=1,nv8
|
|
j = j+1
|
|
right(j) = right(j)-fac0*b0(j0,jj)
|
|
325 continue
|
|
330 j0 = j0+1
|
|
i0 = i0+1
|
|
go to 315
|
|
c take into account that we eliminate the constraints (4)
|
|
335 if(j1-1.gt.iop1) go to 360
|
|
fac1 = h(i1)
|
|
do 340 j=1,mv
|
|
right(j) = right(j)-fac1*a1(j1,j)
|
|
340 continue
|
|
if(mv.eq.mvv) go to 350
|
|
j = mv
|
|
do 345 jj=1,nv8
|
|
j = j+1
|
|
right(j) = right(j)-fac1*b1(j1,jj)
|
|
345 continue
|
|
350 j1 = j1+1
|
|
i1 = i1-1
|
|
go to 335
|
|
360 irot = nrold-iop0-1
|
|
if(irot.lt.0) irot = 0
|
|
c rotate the new row of matrix (auu) into triangle.
|
|
if(i0.gt.i1) go to 390
|
|
do 385 i=i0,i1
|
|
irot = irot+1
|
|
piv = h(i)
|
|
if(piv.eq.0.) go to 385
|
|
c calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,au(irot,1),co,si)
|
|
c apply that transformation to the rows of matrix (qq).
|
|
iq = (irot-1)*mvv
|
|
do 370 j=1,mvv
|
|
iq = iq+1
|
|
call fprota(co,si,right(j),q(iq))
|
|
370 continue
|
|
c apply that transformation to the columns of (auu).
|
|
if(i.eq.i1) go to 385
|
|
i2 = 1
|
|
i3 = i+1
|
|
do 380 j=i3,i1
|
|
i2 = i2+1
|
|
call fprota(co,si,h(j),au(irot,i2))
|
|
380 continue
|
|
385 continue
|
|
c we update the sum of squared residuals.
|
|
390 do 395 j=1,mvv
|
|
sq = sq+right(j)**2
|
|
395 continue
|
|
if(nrold.eq.number) go to 420
|
|
410 nrold = n1
|
|
n1 = n1+1
|
|
go to 250
|
|
420 continue
|
|
if(nuu.eq.0) go to 800
|
|
c we determine the matrix (avv) and then we reduce her to an unit
|
|
c upper triangular form (rv) using givens rotations without square
|
|
c roots. we apply the same transformations to the columns of matrix
|
|
c g to obtain the (nv-7) x (nu-6-iop0-iop1) matrix h.
|
|
c we store matrix (rv) into av1 and av2, h into c.
|
|
c the nv7 x nv7 triangular unit upper matrix (rv) has the form
|
|
c | av1 ' |
|
|
c (rv) = | ' av2 |
|
|
c | 0 ' |
|
|
c with (av2) a nv7 x 4 matrix and (av1) a nv11 x nv11 unit upper
|
|
c triangular matrix of bandwidth 5.
|
|
ncof = nuu*nv7
|
|
c initialization.
|
|
do 430 i=1,ncof
|
|
c(i) = 0.
|
|
430 continue
|
|
do 440 i=1,nv4
|
|
av1(i,5) = 0.
|
|
do 440 j=1,4
|
|
av1(i,j) = 0.
|
|
av2(i,j) = 0.
|
|
440 continue
|
|
jper = 0
|
|
nrold = 0
|
|
do 770 it=1,mv
|
|
number = nrv(it)
|
|
450 if(nrold.eq.number) go to 480
|
|
if(p.le.0.) go to 760
|
|
c fetch a new row of matrix (bv).
|
|
n1 = nrold+1
|
|
do 460 j=1,5
|
|
h(j) = bv(n1,j)*pinv
|
|
460 continue
|
|
c find the appropriate row of g.
|
|
do 465 j=1,nuu
|
|
right(j) = 0.
|
|
465 continue
|
|
if(mv.eq.mvv) go to 510
|
|
l = mv+n1
|
|
do 470 j=1,nuu
|
|
right(j) = q(l)
|
|
l = l+mvv
|
|
470 continue
|
|
go to 510
|
|
c fetch a new row of matrix (spv)
|
|
480 h(5) = 0.
|
|
do 490 j=1,4
|
|
h(j) = spv(it,j)
|
|
490 continue
|
|
c find the appropriate row of g.
|
|
l = it
|
|
do 500 j=1,nuu
|
|
right(j) = q(l)
|
|
l = l+mvv
|
|
500 continue
|
|
c test whether there are non-zero values in the new row of (avv)
|
|
c corresponding to the b-splines n(j;v),j=nv7+1,...,nv4.
|
|
510 if(nrold.lt.nv11) go to 710
|
|
if(jper.ne.0) go to 550
|
|
c initialize the matrix (av2).
|
|
jk = nv11+1
|
|
do 540 i=1,4
|
|
ik = jk
|
|
do 520 j=1,5
|
|
if(ik.le.0) go to 530
|
|
av2(ik,i) = av1(ik,j)
|
|
ik = ik-1
|
|
520 continue
|
|
530 jk = jk+1
|
|
540 continue
|
|
jper = 1
|
|
c if one of the non-zero elements of the new row corresponds to one of
|
|
c the b-splines n(j;v),j=nv7+1,...,nv4, we take account of condition
|
|
c (2) for setting up this row of (avv). the row is stored in h1( the
|
|
c part with respect to av1) and h2 (the part with respect to av2).
|
|
550 do 560 i=1,4
|
|
h1(i) = 0.
|
|
h2(i) = 0.
|
|
560 continue
|
|
h1(5) = 0.
|
|
j = nrold-nv11
|
|
do 600 i=1,5
|
|
j = j+1
|
|
l0 = j
|
|
570 l1 = l0-4
|
|
if(l1.le.0) go to 590
|
|
if(l1.le.nv11) go to 580
|
|
l0 = l1-nv11
|
|
go to 570
|
|
580 h1(l1) = h(i)
|
|
go to 600
|
|
590 h2(l0) = h2(l0) + h(i)
|
|
600 continue
|
|
c rotate the new row of (avv) into triangle.
|
|
if(nv11.le.0) go to 670
|
|
c rotations with the rows 1,2,...,nv11 of (avv).
|
|
do 660 j=1,nv11
|
|
piv = h1(1)
|
|
i2 = min0(nv11-j,4)
|
|
if(piv.eq.0.) go to 640
|
|
c calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,av1(j,1),co,si)
|
|
c apply that transformation to the columns of matrix g.
|
|
ic = j
|
|
do 610 i=1,nuu
|
|
call fprota(co,si,right(i),c(ic))
|
|
ic = ic+nv7
|
|
610 continue
|
|
c apply that transformation to the rows of (avv) with respect to av2.
|
|
do 620 i=1,4
|
|
call fprota(co,si,h2(i),av2(j,i))
|
|
620 continue
|
|
c apply that transformation to the rows of (avv) with respect to av1.
|
|
if(i2.eq.0) go to 670
|
|
do 630 i=1,i2
|
|
i1 = i+1
|
|
call fprota(co,si,h1(i1),av1(j,i1))
|
|
630 continue
|
|
640 do 650 i=1,i2
|
|
h1(i) = h1(i+1)
|
|
650 continue
|
|
h1(i2+1) = 0.
|
|
660 continue
|
|
c rotations with the rows nv11+1,...,nv7 of avv.
|
|
670 do 700 j=1,4
|
|
ij = nv11+j
|
|
if(ij.le.0) go to 700
|
|
piv = h2(j)
|
|
if(piv.eq.0.) go to 700
|
|
c calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,av2(ij,j),co,si)
|
|
c apply that transformation to the columns of matrix g.
|
|
ic = ij
|
|
do 680 i=1,nuu
|
|
call fprota(co,si,right(i),c(ic))
|
|
ic = ic+nv7
|
|
680 continue
|
|
if(j.eq.4) go to 700
|
|
c apply that transformation to the rows of (avv) with respect to av2.
|
|
j1 = j+1
|
|
do 690 i=j1,4
|
|
call fprota(co,si,h2(i),av2(ij,i))
|
|
690 continue
|
|
700 continue
|
|
c we update the sum of squared residuals.
|
|
do 705 i=1,nuu
|
|
sq = sq+right(i)**2
|
|
705 continue
|
|
go to 750
|
|
c rotation into triangle of the new row of (avv), in case the elements
|
|
c corresponding to the b-splines n(j;v),j=nv7+1,...,nv4 are all zero.
|
|
710 irot =nrold
|
|
do 740 i=1,5
|
|
irot = irot+1
|
|
piv = h(i)
|
|
if(piv.eq.0.) go to 740
|
|
c calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,av1(irot,1),co,si)
|
|
c apply that transformation to the columns of matrix g.
|
|
ic = irot
|
|
do 720 j=1,nuu
|
|
call fprota(co,si,right(j),c(ic))
|
|
ic = ic+nv7
|
|
720 continue
|
|
c apply that transformation to the rows of (avv).
|
|
if(i.eq.5) go to 740
|
|
i2 = 1
|
|
i3 = i+1
|
|
do 730 j=i3,5
|
|
i2 = i2+1
|
|
call fprota(co,si,h(j),av1(irot,i2))
|
|
730 continue
|
|
740 continue
|
|
c we update the sum of squared residuals.
|
|
do 745 i=1,nuu
|
|
sq = sq+right(i)**2
|
|
745 continue
|
|
750 if(nrold.eq.number) go to 770
|
|
760 nrold = nrold+1
|
|
go to 450
|
|
770 continue
|
|
c test whether the b-spline coefficients must be determined.
|
|
if(iback.ne.0) return
|
|
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.
|
|
k = 1
|
|
do 780 i=1,nuu
|
|
call fpbacp(av1,av2,c(k),nv7,4,c(k),5,nv)
|
|
k = k+nv7
|
|
780 continue
|
|
c second step: solve the system (cr) (ru)' = (c1).
|
|
k = 0
|
|
do 795 j=1,nv7
|
|
k = k+1
|
|
l = k
|
|
do 785 i=1,nuu
|
|
right(i) = c(l)
|
|
l = l+nv7
|
|
785 continue
|
|
call fpback(au,right,nuu,5,right,nu)
|
|
l = k
|
|
do 790 i=1,nuu
|
|
c(l) = right(i)
|
|
l = l+nv7
|
|
790 continue
|
|
795 continue
|
|
c calculate from the conditions (2)-(3)-(4), the remaining b-spline
|
|
c coefficients.
|
|
800 ncof = nu4*nv4
|
|
j = ncof
|
|
do 805 l=1,nv4
|
|
q(l) = dr01
|
|
q(j) = dr11
|
|
j = j-1
|
|
805 continue
|
|
i = nv4
|
|
j = 0
|
|
if(iop0.eq.0) go to 815
|
|
do 810 l=1,nv4
|
|
i = i+1
|
|
q(i) = c0(l)
|
|
810 continue
|
|
815 if(nuu.eq.0) go to 835
|
|
do 830 l=1,nuu
|
|
ii = i
|
|
do 820 k=1,nv7
|
|
i = i+1
|
|
j = j+1
|
|
q(i) = c(j)
|
|
820 continue
|
|
do 825 k=1,3
|
|
ii = ii+1
|
|
i = i+1
|
|
q(i) = q(ii)
|
|
825 continue
|
|
830 continue
|
|
835 if(iop1.eq.0) go to 845
|
|
do 840 l=1,nv4
|
|
i = i+1
|
|
q(i) = c1(l)
|
|
840 continue
|
|
845 do 850 i=1,ncof
|
|
c(i) = q(i)
|
|
850 continue
|
|
c calculate the quantities
|
|
c res(i,j) = (r(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)
|
|
fp = 0.
|
|
do 890 i=1,nu
|
|
fpu(i) = 0.
|
|
890 continue
|
|
do 900 i=1,nv
|
|
fpv(i) = 0.
|
|
900 continue
|
|
ir = 0
|
|
nroldu = 0
|
|
c main loop for the different grid points.
|
|
do 950 i1=1,mu
|
|
numu = nru(i1)
|
|
numu1 = numu+1
|
|
nroldv = 0
|
|
do 940 i2=1,mv
|
|
numv = nrv(i2)
|
|
numv1 = numv+1
|
|
ir = ir+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.
|
|
k1 = numu*nv4+numv
|
|
do 920 l1=1,4
|
|
k2 = k1
|
|
fac = spu(i1,l1)
|
|
do 910 l2=1,4
|
|
k2 = k2+1
|
|
term = term+fac*spv(i2,l2)*c(k2)
|
|
910 continue
|
|
k1 = k1+nv4
|
|
920 continue
|
|
c calculate the squared residual at the current grid point.
|
|
term = (r(ir)-term)**2
|
|
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 930
|
|
fpv(numv1) = fpv(numv1)-fac
|
|
fpv(numv) = fpv(numv)+fac
|
|
930 nroldv = numv
|
|
if(numu.eq.nroldu) go to 940
|
|
fpu(numu1) = fpu(numu1)-fac
|
|
fpu(numu) = fpu(numu)+fac
|
|
940 continue
|
|
nroldu = numu
|
|
950 continue
|
|
return
|
|
end
|