! D02TGF Example Program Text
! Mark 29.2 Release. NAG Copyright 2023.
Module d02tgfe_mod
! D02TGF 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 :: bdyc, coeff
! .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: coeff_tol = 1.0E-5_nag_wp
Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp
Integer, Parameter, Public :: n = 2, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: x0, x1
Integer, Public, Save :: k1
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable, Public, Save :: b(:,:)
Integer, Public, Save :: l(n) = (/1,2/)
Integer, Public, Save :: m(n) = (/1,2/)
Contains
Subroutine coeff(x,i,a,ia,ia1,rhs)
! .. Use Statements ..
Use nag_library, Only: e02akf
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: rhs
Real (Kind=nag_wp), Intent (In) :: x
Integer, Intent (In) :: i, ia, ia1
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: a(ia,ia1)
! .. Local Scalars ..
Real (Kind=nag_wp) :: z1, z2
Integer :: ifail
! .. Executable Statements ..
If (i<=1) Then
! Evaluate z1, z2 at x using previous coeffs b.
ifail = 0
Call e02akf(k1,x0,x1,b(1,1),1,k1,x,z1,ifail)
Call e02akf(k1,x0,x1,b(1,2),1,k1,x,z2,ifail)
a(1,1) = z2*z2 - one
a(1,2) = two
a(2,1) = two*z1*z2 + one
rhs = two*z1*z2*z2
Else
a(1,2) = -one
a(2,3) = two
End If
Return
End Subroutine coeff
Subroutine bdyc(x,i,j,a,ia,ia1,rhs)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: rhs
Real (Kind=nag_wp), Intent (Out) :: x
Integer, Intent (In) :: i, ia, ia1, j
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: a(ia,ia1)
! .. Executable Statements ..
x = -one
a(i,j) = one
If (i==2 .And. j==1) Then
rhs = 3.0_nag_wp
End If
Return
End Subroutine bdyc
End Module d02tgfe_mod
Program d02tgfe
! D02TGF Example Main Program
! .. Use Statements ..
Use d02tgfe_mod, Only: b, bdyc, coeff, coeff_tol, k1, l, m, n, nin, &
nout, x0, x1
Use nag_library, Only: d02tgf, e02akf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: emax, x, xinc
Integer :: i, ia1, ifail, iter, j, k, kp, ldc, &
liw, lsum, lw, mimax
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: c(:,:), w(:), y(:)
Integer, Allocatable :: iw(:)
! .. Intrinsic Procedures ..
Intrinsic :: abs, max, real, sum
! .. Executable Statements ..
Write (nout,*) 'D02TGF Example Program Results'
! Skip heading in data file
Read (nin,*)
Read (nin,*) mimax, kp
Read (nin,*) x0, x1
lsum = sum(l(1:n))
k1 = mimax + 1
ldc = k1
liw = n*k1
lw = 2*(n*kp+lsum)*(n*k1+1) + 7*n*k1
Allocate (b(k1,n),c(ldc,n),w(lw),y(n),iw(liw))
! Initialize coefficients b(:,:) such that z1 = 0 and z2 = 3.
b(1:k1,1:n) = 0.0_nag_wp
b(1,2) = 6.0_nag_wp
! Iterate until coefficients of linearized systems converge.
iter = 0
emax = 1.0_nag_wp
iters: Do While (emax>=coeff_tol)
iter = iter + 1
Write (nout,*)
Write (nout,99999) ' Iteration', iter, ' Chebyshev coefficients are'
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02tgf(n,m,l,x0,x1,k1,kp,c,ldc,coeff,bdyc,w,lw,iw,liw,ifail)
Write (nout,99998)(c(1:k1,j),j=1,n)
emax = 0.0_nag_wp
Do j = 1, n
Do i = 1, k1
emax = max(emax,abs(c(i,j)-b(i,j)))
End Do
End Do
b(1:k1,1:n) = c(1:k1,1:n)
End Do iters
Deallocate (b)
! Print solution on uniform mesh.
k = 9
ia1 = 1
Write (nout,*)
Write (nout,99999) 'Solution evaluated at', k, ' equally spaced points'
Write (nout,*)
Write (nout,99997) ' X ', (j,j=1,n)
xinc = (x1-x0)/real(k-1,kind=nag_wp)
x = x0
Do i = 1, k
Do j = 1, n
ifail = 0
Call e02akf(k1,x0,x1,c(1,j),ia1,k1,x,y(j),ifail)
End Do
Write (nout,99996) x, (y(j),j=1,n)
x = x + xinc
End Do
99999 Format (1X,A,I3,A)
99998 Format (1X,9F8.4)
99997 Format (1X,A,2(' Y(',I1,')'))
99996 Format (1X,3F10.4)
End Program d02tgfe