```    Program d03edfe

!     D03EDF Example Program Text

!     Mark 27.0 Release. NAG Copyright 2019.

!     .. Use Statements ..
Use nag_library, Only: d03edf, nag_wp, x04abf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Real (Kind=nag_wp), Parameter    :: four = 4.0_nag_wp
Real (Kind=nag_wp), Parameter    :: half = 0.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               :: iset = 1, nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: acc, alpha, hx, hy
Integer                          :: i, ifail, iout, ix, iy1, iy2, j, k,  &
lda, levels, maxit, ngx, ngy, numit, &
outchn
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: a(:,:), rhs(:), u(:), ub(:), us(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: real
!     .. Executable Statements ..
Write (nout,*) 'D03EDF Example Program Results'
!     Skip heading in data file
ngx = 2**levels + 1
ngy = ngx
lda = 4*(ngx+1)*(ngy+1)/3
Allocate (a(lda,7),rhs(lda),u(lda),ub(ngx*ngy),us(lda))

outchn = nout
Write (nout,*)
Call x04abf(iset,outchn)
!     ** Set iout > 2 to obtain intermediate output **
iout = 0
hx = one/real(ngx+1,kind=nag_wp)
hy = one/real(ngy+1,kind=nag_wp)

!     Set up operator, right-hand side and initial guess for
!     step-lengths HX and HY
a(1:ngx*ngy,1) = one - half*alpha
a(1:ngx*ngy,2) = half*alpha
a(1:ngx*ngy,3) = one - half*alpha
a(1:ngx*ngy,4) = -four + alpha
a(1:ngx*ngy,5) = one - half*alpha
a(1:ngx*ngy,6) = half*alpha
a(1:ngx*ngy,7) = one - half*alpha
rhs(1:ngx*ngy) = -four*hx*hy
ub(1:ngx*ngy) = zero

!     Correction for the boundary conditions
!         Horizontal boundaries --
!            Boundary condition on Y=0 -- U=0
a(2:ngx-1,1:2) = zero
!     Boundary condition on Y=1 -- U=0
ix = (ngy-1)*ngx
a(ix+1:ix+ngx-1,6:7) = zero

!     Vertical boundaries --
iy1 = 1
iy2 = ngx
Do j = 2, ngy - 1
!       Boundary condition on X=0 -- U=0
iy1 = iy1 + ngx
a(iy1,3) = zero
a(iy1,6) = zero
!       Boundary condition on X=1 -- U=1
iy2 = iy2 + ngx
rhs(iy2) = rhs(iy2) - a(iy2,5) - a(iy2,2)
a(iy2,2) = zero
a(iy2,5) = zero
End Do

!     Now the four corners --
!     Bottom left corner
a(1,1:3) = zero
a(1,6) = zero
!     Top left corner
a(ix+1,3) = zero
a(ix+1,6:7) = zero
!     Bottom right corner
!     Use average value at discontinuity ( = 0.5 )
k = ngx
rhs(k) = rhs(k) - a(k,2)*half - a(k,5)
a(k,1:2) = zero
a(k,5) = zero
!     Top right corner
k = ngx*ngy
rhs(k) = rhs(k) - a(k,2) - a(k,5)
a(k,2) = zero
a(k,5:7) = zero

!     Solve the equations
!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03edf(ngx,ngy,lda,a,rhs,ub,maxit,acc,us,u,iout,numit,ifail)

Write (nout,99999) ngx, ngy, acc, maxit
Write (nout,*)
Write (nout,99998) 'Residual norm =', us(1)
Write (nout,99997) 'Number of iterations =', numit
Write (nout,*)
Write (nout,*) 'Solution'
Write (nout,*)
Write (nout,99996) '  I/J', (i,i=1,ngx)
Do j = 1, ngy
Write (nout,99995) j, (u(i+(j-1)*ngx),i=1,ngx)
End Do

99999 Format (1X,'NGX = ',I3,'  NGY = ',I3,'  ACC =',1P,E10.2,'  MAXIT',' = ', &
I3)
99998 Format (1X,A,1P,E12.2)
99997 Format (1X,A,I5)
99996 Format (1X,A,10I7,:)
99995 Format (1X,I3,2X,10F7.3,:)
End Program d03edfe
```