! D02TZF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02tzfe_mod ! D02TZF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: mmax = 2, neq = 1, nin = 5, & nlbc = 1, nout = 6, nrbc = 1 ! .. Local Scalars .. REAL (KIND=nag_wp) :: alpha, beta, eps ! .. Local Arrays .. INTEGER :: m(1) = (/ 2/) CONTAINS SUBROUTINE ffun(x,y,neq,m,f) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: f(neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. f(1) = (y(1,0)-y(1,0)*y(1,1))/eps RETURN END SUBROUTINE ffun SUBROUTINE fjac(x,y,neq,m,dfdy) ! .. Use Statements .. USE nag_library, ONLY : x02ajf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: dfdy(neq,neq,0:*) REAL (KIND=nag_wp), INTENT (IN) :: y(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Local Scalars .. REAL (KIND=nag_wp) :: epsh, fac, ptrb INTEGER :: i, j, k ! .. Local Arrays .. REAL (KIND=nag_wp) :: f1(1), f2(1), yp(1,0:3) ! .. Intrinsic Functions .. INTRINSIC abs, max, sqrt ! .. Executable Statements .. epsh = 100.0_nag_wp*x02ajf() fac = sqrt(x02ajf()) DO i = 1, neq DO j = 0, m(i) - 1 yp(i,j) = y(i,j) END DO END DO DO i = 1, neq DO j = 0, m(i) - 1 ptrb = max(epsh,fac*abs(y(i,j))) yp(i,j) = y(i,j) + ptrb CALL ffun(x,yp,neq,m,f1) yp(i,j) = y(i,j) - ptrb CALL ffun(x,yp,neq,m,f2) DO k = 1, neq dfdy(k,i,j) = 0.5_nag_wp*(f1(k)-f2(k))/ptrb END DO yp(i,j) = y(i,j) END DO END DO RETURN END SUBROUTINE fjac SUBROUTINE gafun(ya,neq,m,nlbc,ga) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nlbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: ga(nlbc) REAL (KIND=nag_wp), INTENT (IN) :: ya(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. ga(1) = ya(1,0) - alpha RETURN END SUBROUTINE gafun SUBROUTINE gbfun(yb,neq,m,nrbc,gb) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nrbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: gb(nrbc) REAL (KIND=nag_wp), INTENT (IN) :: yb(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. gb(1) = yb(1,0) - beta RETURN END SUBROUTINE gbfun SUBROUTINE gajac(ya,neq,m,nlbc,dgady) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nlbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: dgady(nlbc,neq,0:*) REAL (KIND=nag_wp), INTENT (IN) :: ya(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. dgady(1,1,0) = one RETURN END SUBROUTINE gajac SUBROUTINE gbjac(yb,neq,m,nrbc,dgbdy) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nrbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: dgbdy(nrbc,neq,0:*) REAL (KIND=nag_wp), INTENT (IN) :: yb(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. dgbdy(1,1,0) = one RETURN END SUBROUTINE gbjac SUBROUTINE guess(x,neq,m,y,dym) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: dym(neq) REAL (KIND=nag_wp), INTENT (INOUT) :: y(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. y(1,0) = alpha + (beta-alpha)*x y(1,1) = (beta-alpha) dym(1) = 0.0_nag_wp RETURN END SUBROUTINE guess END MODULE d02tzfe_mod PROGRAM d02tzfe ! D02TZF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02tkf, d02tvf, d02txf, d02tyf, d02tzf USE d02tzfe_mod, ONLY : alpha, beta, eps, ffun, fjac, gafun, gajac, & gbfun, gbjac, guess, m, mmax, nag_wp, neq, nin, & nlbc, nout, nrbc ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: ermx INTEGER :: i, iermx, ifail, ijermx, j, & liwork, lrwork, mxmesh, ncol, & nmesh LOGICAL :: failed ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: mesh(:), rwork(:), tol(:), y(:,:) INTEGER, ALLOCATABLE :: ipmesh(:), iwork(:) ! .. Executable Statements .. WRITE (nout,*) 'D02TZF Example Program Results' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) ncol, nmesh, mxmesh liwork = mxmesh*(11*neq+6) lrwork = mxmesh*(109*neq**2+78*neq+7) ALLOCATE (mesh(mxmesh),tol(neq),rwork(lrwork),y(neq,0:mmax-1), & ipmesh(mxmesh),iwork(liwork)) READ (nin,*) alpha, beta, eps READ (nin,*) mesh(1:nmesh) READ (nin,*) ipmesh(1:nmesh) READ (nin,*) tol(1:neq) ! Initialize ifail = 0 CALL d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rwork, & lrwork,iwork,liwork,ifail) eps = 0.1_nag_wp*eps CONTN: DO j = 1, 2 WRITE (nout,99997) tol(1), eps ! Solve ifail = -1 CALL d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rwork,iwork, & ifail) failed = ifail /= 0 ! Extract mesh. ifail = -1 CALL d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rwork,iwork, & ifail) ! Print mesh statistics. WRITE (nout,99996) nmesh, ermx, iermx, ijermx IF (failed) EXIT CONTN ! Print solution at every second point on final mesh. WRITE (nout,99999) DO i = 1, nmesh, 2 ifail = -1 CALL d02tyf(mesh(i),y,neq,mmax,rwork,iwork,ifail) WRITE (nout,99998) mesh(i), y(1,0), y(1,1) END DO IF (j==1) THEN ! Halve final mesh for new initial mesh and set up for continuation. nmesh = (nmesh+1)/2 ifail = 0 CALL d02txf(mxmesh,nmesh,mesh,ipmesh,rwork,iwork,ifail) ! Reduce continuation parameter. eps = 0.1_nag_wp*eps END IF END DO CONTN 99999 FORMAT (/' Solution and derivative at every second point:'/ & ' x u u''') 99998 FORMAT (' ',F8.4,2F11.5) 99997 FORMAT (//' Tolerance = ',E8.1,' EPS = ',E10.3) 99996 FORMAT (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, & ' in interval ',I4,' for component ',I4) END PROGRAM d02tzfe