cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc cc cc mnpogr : pogrid test program cc cc cc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c ..local scalars.. real cv,del,ermax,er0,exz0,fp,pi,r,sum,sv,x,y,z0,one,ai,s integer i,ier,is,j,k,kwrk,lwrk,m,mu,mv,nc,nuest,nu,nvest,nv c ..local arrays.. integer ider(2),iopt(3),iwrk(100),iw(29) real u(9),v(20),z(180),c(300),tu(50),tv(50),f(180),wk(116), * exact(180),err(9),sp(9),wrk(1600) c ..function references.. real abs,atan2,cos,sin,tespog c .. c set constants one = 1 pi = atan2(0.,-one) c we set up the radius of the disc r = one c we set up the number of u (radius)-values of the grid. mu = 9 c we set up the u-coordinates of the grid. do 10 i=1,mu ai = i u(i) = ai*0.1 10 continue c we set up the number of v (angle)-values of the grid mv = 20 c we set up the v-coordinates of the grid. del = pi*0.1 do 20 j=1,mv ai = j-1 v(j) = ai*del-pi 20 continue c we fetch the data values at the grid points. m = mu*mv read(5,900) (z(i),i=1,m) c we fetch the data value at the origin. read(5,900) z0 c we print the data values at the grid points. we also compute and print c the exact value of the test function underlying the data. write(6,905) write(6,910) (i,i=1,mu) write(6,915) exz0 = tespog(0.,0.) er0 = abs(exz0-z0) ermax = er0 sum = er0 do 40 j=1,mv cv = cos(v(j)) sv = sin(v(j)) k = j do 30 i=1,mu x = u(i)*cv y = u(i)*sv exact(k) = tespog(x,y) err(i) = abs(exact(k)-z(k)) sum = sum+err(i) if(err(i).gt.ermax) ermax = err(i) k = k+mv 30 continue write(6,920) j,(z(k),k=j,m,mv) write(6,925) (exact(k),k=j,m,mv) 40 continue ai = m+1 sum = sum/ai write(6,930) z0,exz0 write(6,935) sum,ermax c we set up the dimension information nuest = 16 nvest = 27 kwrk = 100 lwrk = 1600 c main loop for the different spline approximations do 300 is=1,6 go to (110,120,130,140,150,160),is c we start computing a set of spline approximations with c only c0-continuity at the origin, 110 iopt(2) = 0 ider(2) = 0 c non-vanishing at the boundary of the disc, iopt(3) = 0 c with a data value at the origin. ider(1) = 0 c initialisation iopt(1) = 0 c a large value for s for computing the least-squares polynomial s = 5. go to 200 c iopt(1) = 1 from the second call on 120 s = 0.1 iopt(1) = 1 go to 200 c an interpolating spline 130 s = 0. go to 200 c a second set of approximations with c1-continuity at the origin 140 iopt(2) = 1 c vanishing at the boundary of the disc. iopt(3) = 1 c exact value at the origin. ider(1) = 1 z0 = exz0 c reinitialization iopt(1) = 0 s = 0.1 go to 200 c no data value at the origin 150 ider(1) = -1 c vanishing partial derivatives at the origin ider(2) = 1 c reinitialization iopt(1) = 0 go to 200 c finally we calculate the least-squares spline according to the current c set of knots 160 iopt(1) = -1 200 call pogrid(iopt,ider,mu,u,mv,v,z,z0,r,s,nuest,nvest, * nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier) c printing of the fitting results. if(iopt(1).ge.0) go to 210 write(6,940) go to 220 210 write(6,945) s 220 write(6,950) iopt(2) if(ider(2).eq.1) write(6,955) if(iopt(3).eq.1) write(6,960) write(6,965) fp,ier write(6,970) nu write(6,975) write(6,980) (tu(i),i=1,nu) write(6,985) nv write(6,975) write(6,980) (tv(i),i=1,nv) nc = (nu-4)*(nv-4) write(6,990) write(6,980) (c(i),i=1,nc) c evaluation of the spline approximation call bispev(tu,nu,tv,nv,c,3,3,u,mu,v,mv,f,wk,116,iw,29,ier) write(6,995) write(6,910) (i,i=1,mu,2) write(6,915) er0 = abs(exz0-c(1)) ermax = er0 sum = er0 do 240 j=1,mv k = j do 230 i=1,mu sp(i) = f(k) err(i) = abs(exact(k)-f(k)) sum = sum+err(i) if(err(i).gt.ermax) ermax = err(i) k = k+mv 230 continue if( (j/3)*3 .ne.j ) go to 240 write(6,920) j,(sp(i),i=1,mu,2) write(6,925) (err(i),i=1,mu,2) 240 continue sum = sum/ai write(6,1000) c(1),er0 write(6,935) sum,ermax 300 continue stop 900 format(10f8.3) 905 format(49h1data value (exact function value) at grid points) 910 format(8h u(i),i=,3x,9(i1,7x)) 915 format(8h v(j),j=) 920 format(1x,i5,9(2x,f6.3)) 925 format(7x,9(2h (,f5.3,1h))) 930 format(23h0data value at (0,0) = ,f7.3,5x,14hexact value = ,f7.3) 935 format(19h0mean abs. error = ,f9.3,5x,18hmax. abs. error = ,f9.3) 940 format(21h0least-squares spline) 945 format(25h0smoothing spline with s=,f7.2) 950 format(1x,35horder of continuity at the origin =,i3) 955 format(1x,43hvanishing partial derivatives at the origin) 960 format(1x,37hvanishing at the boundary of the disc) 965 format(1x,23hsum squared residuals =,e15.6,5x,11herror flag=,i3) 970 format(1x,42htotal number of knots in the u-direction =,i3) 975 format(1x,22hposition of the knots ) 980 format(5x,8f9.4) 985 format(1x,42htotal number of knots in the v-direction =,i3) 990 format(23h0b-spline coefficients ) 995 format(50h0spline value (approximation error) at grid points) 1000 format(25h0spline value at (0,0) = ,f7.3,5x,8herror = ,f7.3) end c real function tespog(x,y) c function program tespog calculates the value of the test function c underlying the data. c .. c ..scalar arguments.. real x,y,f c .. f = 1.-((3.*x-1.)**2+(3.*y-1.)**2)/(11.-6.*(x+y)) tespog = f-(1.-x**2-y**2)*(x+y)*54./121. return end