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

    MODULE d03eefe_mod

!      D03EEF 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       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: zero = 0.0_nag_wp
       INTEGER, PARAMETER                  :: nin = 5, nout = 6
    CONTAINS
       SUBROUTINE pdef(x,y,alpha,beta,gamma,delta,epslon,phi,psi)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: alpha, beta, delta, epslon,   &
                                                 gamma, phi, psi
          REAL (KIND=nag_wp), INTENT (IN)     :: x, y
!         .. Intrinsic Functions ..
          INTRINSIC                              cos, sin
!         .. Executable Statements ..
          alpha = one
          beta = zero
          gamma = one
          delta = 50.0_nag_wp
          epslon = 50.0_nag_wp
          phi = zero

          psi = sin(x)*((-alpha-gamma+phi)*sin(y)+epslon*cos(y)) + &
             cos(x)*(delta*sin(y)+beta*cos(y))

          RETURN
       END SUBROUTINE pdef
       SUBROUTINE bndy(x,y,a,b,c,ibnd)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          INTEGER, PARAMETER                  :: bottom = 0, left = 3,         &
                                                 right = 1, top = 2
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: a, b, c
          REAL (KIND=nag_wp), INTENT (IN)     :: x, y
          INTEGER, INTENT (IN)                :: ibnd
!         .. Intrinsic Functions ..
          INTRINSIC                              sin
!         .. Executable Statements ..
          IF (ibnd==top .OR. ibnd==right) THEN

!            Solution prescribed

             a = one
             b = zero
             c = sin(x)*sin(y)
          ELSE IF (ibnd==bottom) THEN

!            Derivative prescribed

             a = zero
             b = one
             c = -sin(x)
          ELSE IF (ibnd==left) THEN

!            Derivative prescribed

             a = zero
             b = one
             c = -sin(y)
          END IF

          RETURN
       END SUBROUTINE bndy
       FUNCTION fexact(x,y)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: fexact
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x, y
!         .. Intrinsic Functions ..
          INTRINSIC                              sin
!         .. Executable Statements ..
          fexact = sin(x)*sin(y)
          RETURN
       END FUNCTION fexact
    END MODULE d03eefe_mod

    PROGRAM d03eefe

!      D03EEF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03edf, d03eef, nag_wp
       USE d03eefe_mod, ONLY : bndy, fexact, nin, nout, pdef, zero
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: acc, hx, hy, rmserr, xmax, xmin, &
                                              xx, ymax, ymin, yy
       INTEGER                             :: i, icase, ifail, iout, ix, j,    &
                                              lda, levels, maxit, ngx, ngxy,   &
                                              ngy, numit
       CHARACTER (7)                       :: scheme
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: a(:,:), rhs(:), u(:), ub(:),     &
                                              us(:), x(:), y(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real, sqrt
!      .. Executable Statements ..
       WRITE (nout,*) 'D03EEF Example Program Results'
       WRITE (nout,*)
       FLUSH (nout)

!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) levels
       ngx = 2**levels + 1
       ngy = ngx
       lda = 4*(ngx+1)*(ngy+1)/3
       ngxy = ngx*ngy
       ALLOCATE (a(lda,7),rhs(lda),u(lda),ub(ngxy),us(lda),x(ngxy),y(ngxy))
       READ (nin,*) xmin, xmax
       READ (nin,*) ymin, ymax
       hx = (xmax-xmin)/real(ngx-1,kind=nag_wp)
       DO i = 1, ngx
          xx = xmin + real(i-1,kind=nag_wp)*hx
          x(i:ngxy:ngx) = xx
       END DO
       hy = (ymax-ymin)/real(ngy-1,kind=nag_wp)
       DO j = 1, ngy
          yy = ymin + real(j-1,kind=nag_wp)*hy
          y((j-1)*ngx+1:j*ngx) = yy
       END DO
!      ** set iout > 2 to obtain intermediate output from D03EDF **
       iout = 0
       READ (nin,*) acc
       READ (nin,*) maxit


CASES: DO icase = 1, 2

          SELECT CASE (icase)
          CASE (1)
!            Central differences
             scheme = 'Central'
          CASE (2)
!            Upwind differences
             scheme = 'Upwind'
          END SELECT
!         Discretize the equations
!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft     
          ifail = -1
          CALL d03eef(xmin,xmax,ymin,ymax,pdef,bndy,ngx,ngy,lda,a,rhs,scheme, &
             ifail)

          IF (ifail<0) THEN
             WRITE (nout,99995) ifail
             EXIT CASES
          END IF

!         Set the initial guess to zero
          ub(1:ngxy) = zero

!         Solve the equations
          ifail = 0
          CALL d03edf(ngx,ngy,lda,a,rhs,ub,maxit,acc,us,u,iout,numit,ifail)

!         Print out the solution
          WRITE (nout,*)
          WRITE (nout,*) 'Exact solution above computed solution'
          WRITE (nout,*)
          WRITE (nout,99998) '  I/J', (i,i=1,ngx)
          rmserr = zero
          DO j = ngy, 1, -1
             ix = (j-1)*ngx
             WRITE (nout,*)
             WRITE (nout,99999) j, (fexact(x(ix+i),y(ix+i)),i=1,ngx)
             WRITE (nout,99999) j, u(ix+1:ix+ngx)
             DO i = 1, ngx
                rmserr = rmserr + (fexact(x(ix+i),y(ix+i))-u(ix+i))**2
             END DO
          END DO
          rmserr = sqrt(rmserr/real(ngxy,kind=nag_wp))
          WRITE (nout,*)
          WRITE (nout,99997) 'Number of Iterations = ', numit
          WRITE (nout,99996) 'RMS Error = ', rmserr
       END DO CASES

99999  FORMAT (1X,I3,2X,10F7.3:/(6X,10F7.3))
99998  FORMAT (1X,A,10I7:/(6X,10I7))
99997  FORMAT (1X,A,I3)
99996  FORMAT (1X,A,1P,E10.2)
99995  FORMAT (1X,' ** D03EEF returned with IFAIL = ',I5)
    END PROGRAM d03eefe