! D02TGF Example Program Text ! Mark 23 Release. NAG Copyright 2011. 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 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 USE d02tgfe_mod, ONLY : b, bdyc, coeff, coeff_tol, k1, l, m, n, nag_wp, & 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 Functions .. 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