! D02TYF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
Module d02tyfe_mod
! D02TYF 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, nmesh_out = 11, &
nout = 6, nrbc = 1
! .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: a
! .. 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)
! .. Intrinsic Procedures ..
Intrinsic :: sqrt
! .. Executable Statements ..
If (y(1,0)<=0.0E0_nag_wp) Then
f(1) = 0.0_nag_wp
Else
f(1) = (y(1,0))**1.5_nag_wp/sqrt(x)
End If
Return
End Subroutine ffun
Subroutine fjac(x,y,neq,m,dfdy,iuser,ruser)
! .. 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)
! .. Intrinsic Procedures ..
Intrinsic :: sqrt
! .. Executable Statements ..
If (y(1,0)<=0.0E0_nag_wp) Then
dfdy(1,1,0) = 0.0_nag_wp
Else
dfdy(1,1,0) = 1.5_nag_wp*sqrt(y(1,0))/sqrt(x)
End If
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) - 1.0_nag_wp
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)
Return
End Subroutine gbfun
Subroutine gajac(ya,neq,m,nlbc,dgady,iuser,ruser)
! .. 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) = 1.0_nag_wp
Return
End Subroutine gajac
Subroutine gbjac(yb,neq,m,nrbc,dgbdy,iuser,ruser)
! .. 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) = 1.0_nag_wp
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) = 1.0_nag_wp - x/a
y(1,1) = -1.0_nag_wp/a
dym(1) = 0.0_nag_wp
Return
End Subroutine guess
End Module d02tyfe_mod
Program d02tyfe
! D02TYF Example Main Program
! .. Use Statements ..
Use nag_library, Only: d02tlf, d02tvf, d02tyf, d02tzf, nag_wp
Use d02tyfe_mod, Only: a, ffun, fjac, gafun, gajac, gbfun, gbjac, guess, &
m, mmax, neq, nin, nlbc, nmesh_out, nout, nrbc
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: ainc, ermx, x
Integer :: i, iermx, ifail, ijermx, 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)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*) 'D02TYF 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,*) a
Read (nin,*) tol(1:neq)
ainc = a/real(nmesh-1,kind=nag_wp)
mesh(1) = 0.0E0_nag_wp
Do i = 2, nmesh - 1
mesh(i) = mesh(i-1) + ainc
End Do
mesh(nmesh) = a
ipmesh(1) = 1
ipmesh(2:nmesh-1) = 2
ipmesh(nmesh) = 1
! 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)
Write (nout,99999) tol(1), a
! 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)
If (ifail/=1) Then
! Print mesh statistics
Write (nout,99998) nmesh, ermx, iermx, ijermx
Write (nout,99997)(i,ipmesh(i),mesh(i),i=1,nmesh)
End If
If (.Not. failed) Then
! Print solution on output mesh.
Write (nout,99996)
x = 0.0_nag_wp
ainc = a/real(nmesh_out-1,kind=nag_wp)
Do i = 1, nmesh_out
ifail = 0
Call d02tyf(x,y,neq,mmax,rcomm,icomm,ifail)
Write (nout,99995) x, y(1,0:1)
x = x + ainc
End Do
End If
99999 Format (//' Tolerance = ',E8.1,' A = ',F8.2)
99998 Format (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, &
' in interval ',I4,' for component ',I4/)
99997 Format (/' Mesh points:'/4(I4,'(',I1,')',E11.4))
99996 Format (/' Computed solution'/' x solution derivative')
99995 Format (' ',F8.2,2F11.5)
End Program d02tyfe