Program d03edfe

!     D03EDF Example Program Text

!     Mark 29.1 Release. NAG Copyright 2023.

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