! D02TGF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. Parameters .. Real (Kind=nag_wp), Parameter :: 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 :: n = 2, nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: x0, x1 Integer :: k1 ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: b(:,:) Integer :: l(n) = (/1,2/) Integer :: 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 ! Evaulate 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) rhs = 3.0_nag_wp Return End Subroutine bdyc End Module d02tgfe_mod Program d02tgfe ! D02TGF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02tgf, e02akf, nag_wp Use d02tgfe_mod, Only: b, bdyc, coeff, coeff_tol, k1, l, m, n, nin, & nout, x0, x1 ! .. 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