!   D02UAF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d02uafe_mod

!      D02UAF 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       :: a = 0.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: b = 1.5_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: zero = 0.0_nag_wp
       INTEGER, PARAMETER                  :: m = 1, nin = 5, nout = 6
       LOGICAL, PARAMETER                  :: reqerr = .FALSE.
    CONTAINS
       FUNCTION exact(x)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: exact
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..	
          exact = exp(-x-one) + one
          RETURN
       END FUNCTION exact
       FUNCTION deriv(x)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: deriv
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..	
          deriv = -exp(-x-one)
          RETURN
       END FUNCTION deriv
       SUBROUTINE bndary(m,y,bmat,bvec)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: m
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: bmat(m,m+1), bvec(m), y(m)
!         .. Executable Statements ..	
!         Boundary condition on left side of domain
          y(1) = a
!         Set up Dirichlet condition using exact solution at x=a.
          bmat(1:m,1:m+1) = zero
          bmat(1,1) = one
          bvec(1) = exact(a)

          RETURN
       END SUBROUTINE bndary
       SUBROUTINE pdedef(m,f)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: m
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(m+1)
!         .. Executable Statements ..	
          f(1:m+1) = zero
          f(1) = one
          f(2) = one

          RETURN
       END SUBROUTINE pdedef

    END MODULE d02uafe_mod
    PROGRAM d02uafe

!      D02UAF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02uaf, d02ubf, d02ucf, d02uef, x02ajf
       USE d02uafe_mod, ONLY : a, b, bndary, deriv, exact, m, nag_wp, nin,     &
                               nout, one, pdedef, reqerr, zero
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: resid, teneps, uerr, uxerr
       INTEGER                             :: i, ifail, n
!      .. Local Arrays ..
       REAL (KIND=nag_wp)                  :: bmat(m,m+1), bvec(m), f(m+1), y(m)
       REAL (KIND=nag_wp), ALLOCATABLE     :: c(:), f0(:), u(:), uc(:,:),      &
                                              ux(:), x(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              abs, int, max
!      .. Executable Statements ..	

       WRITE (nout,*) ' D02UAF Example Program Results '
       WRITE (nout,*)

       READ (nin,*)
       READ (nin,*) n

       ALLOCATE (u(n+1),ux(n+1),f0(n+1),c(n+1),uc(n+1,m+1),x(n+1))

!      Setup problem boundary conditions and definition
       CALL bndary(m,y,bmat,bvec)
       CALL pdedef(m,f)

!      Set up solution grid.
       ifail = 0
       CALL d02ucf(n,a,b,x,ifail)

!      Set up problem right hand sides for grid and transform.
       f0(1:n+1) = one
       ifail = 0
       CALL d02uaf(n,f0,c,ifail)

!      Solve in coefficient space.
       ifail = 0
       CALL d02uef(n,a,b,m,c,bmat,y,bvec,f,uc,resid,ifail)
!      Transform solution and derivative back to real space.
       ifail = 0
       CALL d02ubf(n,a,b,0,uc(1,1),u,ifail)
       ifail = 0
       CALL d02ubf(n,a,b,1,uc(1,2),ux,ifail)

!      Print solution.
       WRITE (nout,*) ' Numerical solution U and derivative Ux'
       WRITE (nout,*)
       WRITE (nout,99999)
       DO i = 1, n + 1
          WRITE (nout,99998) x(i), u(i), ux(i)
       END DO
       IF (reqerr) THEN
          uerr = zero
          uxerr = zero
          DO i = 1, n + 1
             uerr = max(uerr,abs(u(i)-exact(x(i))))
             uxerr = max(uxerr,abs(ux(i)-deriv(x(i))))
          END DO
          teneps = 10.0_nag_wp*x02ajf()
          WRITE (nout,99997) 10*(int(uerr/teneps)+1)
          WRITE (nout,99996) 10*(int(uxerr/teneps)+1)
       END IF

99999  FORMAT (1X,T8,'x',T18,'U',T28,'Ux')
99998  FORMAT (1X,3F10.4)
99997  FORMAT (//1X,'U  is within a multiple ',I8,' of machine precision.')
99996  FORMAT (1X,'Ux is within a multiple ',I8,' of machine precision.')

    END PROGRAM d02uafe