! E02DCF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE e02dcfe_mod ! E02DCF 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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,8F9.4) END SUBROUTINE cprint END MODULE e02dcfe_mod PROGRAM e02dcfe ! E02DCF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : e02dcf, e02dff, nag_wp USE e02dcfe_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, mx, & my, npx, npy, nx, nxest, ny, nyest CHARACTER (1) :: start ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: c(:), f(:), fg(:), lamda(:), & mu(:), px(:), py(:), wrk(:), & x(:), y(:) INTEGER, ALLOCATABLE :: iwrk(:) ! .. Intrinsic Functions .. INTRINSIC max, min, real ! .. Executable Statements .. WRITE (nout,*) 'E02DCF Example Program Results' ! Skip heading in data file READ (nin,*) ! Input the number of X, Y co-ordinates MX, MY. READ (nin,*) mx, my nxest = mx + 4 nyest = my + 4 lwrk = 4*(mx+my) + 11*(nxest+nyest) + nxest*my + max(my,nxest) + 54 liwrk = 3 + mx + my + nxest + nyest ALLOCATE (x(mx),y(my),f(mx*my),lamda(nxest),mu(nyest),wrk(lwrk), & iwrk(liwrk),c((nxest-4)*(nyest-4))) READ (nin,*) x(1:mx) READ (nin,*) y(1:my) ! Input the MX*MY function values F at the grid points. READ (nin,*) f(1:mx*my) start = 'C' READ (nin,*) s ! Determine the spline approximation. ifail = 0 CALL e02dcf(start,mx,x,my,y,f,s,nxest,nyest,nx,lamda,ny,mu,c,fp,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,*) 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,99997) j, lamda(j), j, mu(j) ELSE IF (j<=nx-3) THEN WRITE (nout,99997) j, lamda(j) ELSE IF (j<=ny-3) THEN WRITE (nout,99996) j, mu(j) END IF END DO CALL cprint(c,ny,nx,nout) WRITE (nout,*) WRITE (nout,99998) 'Sum of squared residuals FP =', fp IF (fp==0.0E+0_nag_wp) THEN WRITE (nout,*) '(The spline is an interpolating spline)' ELSE 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 delta = (yhi-ylo)/real(npy-1,kind=nag_wp) 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,99995) ' X', (px(i),i=1,npx) WRITE (nout,*) ' Y' DO i = npy, 1, -1 WRITE (nout,99994) 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,1P,E13.4) 99997 FORMAT (1X,I16,F12.4,I11,F12.4) 99996 FORMAT (1X,I39,F12.4) 99995 FORMAT (1X,A,7F8.2) 99994 FORMAT (1X,F8.2,3X,7F8.2) END PROGRAM e02dcfe