!   D03EEF Example Program Text
!   Mark 28.3 Release. NAG Copyright 2022.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: bndy, fexact, pdef
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: nin = 5, nout = 6
    Contains
      Subroutine pdef(x,y,alpha,beta,gamma,delta,epslon,phi,psi)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: alpha, beta, delta, epslon, gamma, &
                                          phi, psi
        Real (Kind=nag_wp), Intent (In) :: x, y
!       .. Intrinsic Procedures ..
        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)

!       .. 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 Procedures ..
        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)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: fexact
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x, y
!       .. Intrinsic Procedures ..
        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 d03eefe_mod, Only: bndy, fexact, nin, nout, pdef, zero
      Use nag_library, Only: d03edf, d03eef, nag_wp
!     .. 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 Procedures ..
      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