! E02DDF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module e02ddfe_mod ! E02DDF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 Contains Subroutine cprint(c,ny,nx,nout) ! .. Scalar Arguments .. Integer, Intent (In) :: nout, nx, ny ! .. Array Arguments .. Real (Kind=nag_wp), Intent (In) :: c(ny-4,nx-4) ! .. Local Scalars .. Integer :: i ! .. Executable Statements .. Write (nout,*) Write (nout,*) 'The B-spline coefficients:' Write (nout,*) Do i = 1, ny - 4 Write (nout,99999) c(i,1:(nx-4)) End Do Return 99999 Format (1X,7F9.2) End Subroutine cprint End Module e02ddfe_mod Program e02ddfe ! E02DDF Example Main Program ! .. Use Statements .. Use nag_library, Only: e02ddf, e02dff, nag_wp Use e02ddfe_mod, Only: cprint, nin, nout ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: delta, fp, s, xhi, xlo, yhi, ylo Integer :: i, ifail, j, liwrk, lwrk, m, & npx, npy, nx, nxest, ny, nyest, & rank, u, v, ww Character (1) :: start ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: c(:), f(:), fg(:), lamda(:), & mu(:), px(:), py(:), w(:), & wrk(:), x(:), y(:) Integer, Allocatable :: iwrk(:) ! .. Intrinsic Procedures .. Intrinsic :: max, min, real ! .. Executable Statements .. Write (nout,*) 'E02DDF Example Program Results' ! Skip heading in data file Read (nin,*) Read (nin,*) m nxest = m nyest = nxest liwrk = m + 2*(nxest-7)*(nyest-7) u = nxest - 4 v = nyest - 4 ww = max(u,v) lwrk = (7*u*v+25*ww)*(ww+1) + 2*(u+v+4*m) + 23*ww + 56 Allocate (x(m),y(m),f(m),w(m),lamda(nxest),mu(nyest),c((nxest-4)*(nyest- & 4)),iwrk(liwrk),wrk(lwrk)) ! Input the data-points and the weights. Do i = 1, m Read (nin,*) x(i), y(i), f(i), w(i) End Do start = 'C' Read (nin,*) s ! Determine the spline approximation. ifail = 0 Call e02ddf(start,m,x,y,f,w,s,nxest,nyest,nx,lamda,ny,mu,c,fp,rank,wrk, & lwrk,iwrk,liwrk,ifail) Deallocate (wrk,iwrk) Write (nout,*) Write (nout,99999) 'Calling with smoothing factor S =', s, ': NX =', nx, & ', NY =', ny, ',' Write (nout,99998) 'rank deficiency =', (nx-4)*(ny-4) - rank Write (nout,*) Write (nout,*) ' I Knot LAMDA(I) J Knot MU(J)' Write (nout,*) Do j = 4, max(nx,ny) - 3 If (j<=nx-3 .And. j<=ny-3) Then Write (nout,99996) j, lamda(j), j, mu(j) Else If (j<=nx-3) Then Write (nout,99996) j, lamda(j) Else If (j<=ny-3) Then Write (nout,99995) j, mu(j) End If End Do Call cprint(c,ny,nx,nout) Write (nout,*) Write (nout,99997) ' Sum of squared residuals FP =', fp If (nx==8 .And. ny==8) Then Write (nout,*) & ' ( The spline is the least-squares bi-cubic polynomial )' End If ! Evaluate the spline on a rectangular grid at NPX*NPY points ! over the domain (XLO to XHI) x (YLO to YHI). Read (nin,*) npx, xlo, xhi Read (nin,*) npy, ylo, yhi lwrk = min(4*npx+nx,4*npy+ny) If (4*npx+nx>4*npy+ny) Then liwrk = npy + ny - 4 Else liwrk = npx + nx - 4 End If Allocate (px(npx),py(npy),fg(npx*npy),wrk(lwrk),iwrk(liwrk)) delta = (xhi-xlo)/real(npx-1,kind=nag_wp) Do i = 1, npx px(i) = min(xlo+real(i-1,kind=nag_wp)*delta,xhi) End Do Do i = 1, npy py(i) = min(ylo+real(i-1,kind=nag_wp)*delta,yhi) End Do ifail = 0 Call e02dff(npx,npy,nx,ny,px,py,lamda,mu,c,fg,wrk,lwrk,iwrk,liwrk,ifail) Write (nout,*) Write (nout,*) 'Values of computed spline:' Write (nout,*) Write (nout,99994) ' X', (px(i),i=1,npx) Write (nout,*) ' Y' Do i = npy, 1, -1 Write (nout,99993) py(i), (fg(npy*(j-1)+i),j=1,npx) End Do 99999 Format (1X,A,1P,E13.4,A,I5,A,I5,A) 99998 Format (1X,A,I5) 99997 Format (1X,A,1P,E13.4,A) 99996 Format (1X,I16,F12.4,I11,F12.4) 99995 Format (1X,I39,F12.4) 99994 Format (1X,A,7F8.2) 99993 Format (1X,F8.2,3X,7F8.2) End Program e02ddfe