PROGRAM d03edfe

!      D03EDF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. 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 Functions ..
       INTRINSIC                          real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03EDF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) levels
       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,*)
       READ (nin,*) alpha, acc
       READ (nin,*) maxit
       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