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

    MODULE d02ubfe_mod

!      D02UBF Example Program Module:
!             Parameter and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER       :: four = 4.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: three = 3.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: two = 2.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: zero = 0.0_nag_wp
       INTEGER, PARAMETER                  :: m = 3, nin = 5, nout = 6
       LOGICAL, PARAMETER                  :: reqerr = .FALSE.
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: a, b
    CONTAINS
       FUNCTION exact(x,q)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: exact
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
          INTEGER, INTENT (IN)                :: q
!         .. Local Scalars ..
          INTEGER                             :: lq
!         .. Intrinsic Functions ..
          INTRINSIC                              cos, mod, sin
!         .. Executable Statements ..
          lq = mod(q,4)
          SELECT CASE (lq)
          CASE (0)
             exact = cos(x)
          CASE (1)
             exact = -sin(x)
          CASE (2)
             exact = -cos(x)
          CASE (3)
             exact = sin(x)
          END SELECT
          RETURN
       END FUNCTION exact
       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:2) = a
          y(3) = b
!         Set up Dirichlet condition using exact solution at x=a.
          bmat(1:m,1:m+1) = zero
          bmat(1:3,1) = one
          bmat(2:3,2) = two
          bmat(2:3,3) = three
          bvec(1) = zero
          bvec(2) = two
          bvec(3) = -two

          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) = one
          f(2) = two
          f(3) = three
          f(4) = four

          RETURN
       END SUBROUTINE pdedef

    END MODULE d02ubfe_mod
    PROGRAM d02ubfe

!      D02UBF Example Main Program

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

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

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

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

!      Set up problem boundary conditions and definition
       pi = x01aaf(one)
       a = -pi/two
       b = pi/two
       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) = two*sin(x(1:n+1)) - two*cos(x(1:n+1))
       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 derivatives back to real space
       DO q = 0, m
          ifail = 0
          CALL d02ubf(n,a,b,q,uc(1,q+1),u(1:n+1,q+1),ifail)
       END DO

!      Print solution
       WRITE (nout,*) ' Numerical Solution U and its first three derivatives'
       WRITE (nout,*)
       WRITE (nout,99999)
       DO i = 1, n + 1
          WRITE (nout,99998) x(i), u(i,1:m+1)
       END DO

       IF (reqerr) THEN
          uerr(1:m+1) = zero
          DO i = 1, n + 1
             DO q = 0, m
                q1 = q + 1
                uerr(q1) = max(uerr(q1),abs(u(i,q1)-exact(x(i),q)))
             END DO
          END DO
          teneps = 10.0_nag_wp*x02ajf()
          WRITE (nout,'(//)')
          DO q = 0, m
             WRITE (nout,99997) q, 10*(int(uerr(q+1)/teneps)+1)
          END DO
       END IF

99999  FORMAT (1X,T8,'X',T18,'U',T28,'Ux',T37,'Uxx',T47,'Uxxx')
99998  FORMAT (1X,5F10.4)
99997  FORMAT (1X,'Error in the order ',I1,' derivative of U is < ',I8, &
          ' * machine precision.')

    END PROGRAM d02ubfe