! D02TZF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
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
! .. Accessibility Statements ..
Private
Public :: ffun, fjac, gafun, gajac, gbfun, &
gbjac, guess
! .. Parameters ..
Integer, Parameter, Public :: mmax = 2, neq = 1, nin = 5, &
nlbc = 1, nout = 6, nrbc = 1
! .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: alpha, beta, eps
! .. Local Arrays ..
Integer, Public, Save :: m(1) = (/2/)
Contains
Subroutine ffun(x,y,neq,m,f,iuser,ruser)
! .. 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 (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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,iuser,ruser)
! .. Use Statements ..
Use nag_library, Only: x02ajf
! .. 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:*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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 Procedures ..
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,iuser,ruser)
yp(i,j) = y(i,j) - ptrb
Call ffun(x,yp,neq,m,f2,iuser,ruser)
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,iuser,ruser)
! .. Scalar Arguments ..
Integer, Intent (In) :: neq, nlbc
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: ga(nlbc)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
Integer, Intent (In) :: m(neq)
! .. Executable Statements ..
ga(1) = ya(1,0) - alpha
Return
End Subroutine gafun
Subroutine gbfun(yb,neq,m,nrbc,gb,iuser,ruser)
! .. Scalar Arguments ..
Integer, Intent (In) :: neq, nrbc
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: gb(nrbc)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
Integer, Intent (In) :: m(neq)
! .. Executable Statements ..
gb(1) = yb(1,0) - beta
Return
End Subroutine gbfun
Subroutine gajac(ya,neq,m,nlbc,dgady,iuser,ruser)
! .. 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:*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
Integer, Intent (In) :: m(neq)
! .. Executable Statements ..
dgady(1,1,0) = one
Return
End Subroutine gajac
Subroutine gbjac(yb,neq,m,nrbc,dgbdy,iuser,ruser)
! .. 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:*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
Integer, Intent (In) :: m(neq)
! .. Executable Statements ..
dgbdy(1,1,0) = one
Return
End Subroutine gbjac
Subroutine guess(x,neq,m,y,dym,iuser,ruser)
! .. 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) :: ruser(*), y(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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 d02tzfe_mod, Only: alpha, beta, eps, ffun, fjac, gafun, gajac, &
gbfun, gbjac, guess, m, mmax, neq, nin, nlbc, &
nout, nrbc
Use nag_library, Only: d02tlf, d02tvf, d02txf, d02tyf, d02tzf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: ermx
Integer :: i, iermx, ifail, ijermx, j, licomm, &
lrcomm, mxmesh, ncol, nmesh
Logical :: failed
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: mesh(:), rcomm(:), tol(:), y(:,:)
Real (Kind=nag_wp) :: ruser(1)
Integer, Allocatable :: icomm(:), ipmesh(:)
Integer :: iuser(2)
! .. Executable Statements ..
Write (nout,*) 'D02TZF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) ncol, nmesh, mxmesh
Allocate (mesh(mxmesh),tol(neq),y(neq,0:mmax-1),ipmesh(mxmesh))
Read (nin,*) alpha, beta, eps
Read (nin,*) mesh(1:nmesh)
Read (nin,*) ipmesh(1:nmesh)
Read (nin,*) tol(1:neq)
! Workspace query to get size of rcomm and icomm
ifail = 0
Call d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,ruser,0, &
iuser,2,ifail)
lrcomm = iuser(1)
licomm = iuser(2)
Allocate (rcomm(lrcomm),icomm(licomm))
! Initialize
ifail = 0
Call d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rcomm, &
lrcomm,icomm,licomm,ifail)
eps = 0.1_nag_wp*eps
contn: Do j = 1, 2
Write (nout,99997) tol(1), eps
! Solve
ifail = -1
Call d02tlf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rcomm,icomm,iuser, &
ruser,ifail)
failed = ifail /= 0
! Extract mesh.
ifail = -1
Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rcomm,icomm, &
ifail)
! Print mesh statistics.
Write (nout,99996) nmesh, ermx, iermx, ijermx
If (failed) Then
Exit contn
End If
! 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,rcomm,icomm,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,rcomm,icomm,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