    PROGRAM nag_pde_ell_mg_ex03

      ! Example Program Text for nag_pde_ell_mg
      ! NAG fl90, Release 4. NAG Copyright 2000.

      ! .. Use Statements ..
      USE nag_examples_io, ONLY : nag_std_out
      USE nag_pde_ell_mg, ONLY : nag_pde_ell_mg_sol
      USE nag_math_constants, ONLY : nag_pi
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC EPSILON, KIND, MAX, MIN, PRECISION, REAL, SIN, SQRT, SUM
      ! .. Parameters ..
      INTEGER, PARAMETER :: level_low = 3
      INTEGER, PARAMETER :: wp = KIND(1.0D0)
      INTEGER, PARAMETER :: level_up = MIN(8,PRECISION(0.0_wp))
      INTEGER, PARAMETER :: maxsize = (2**level_up+1)**2
      REAL (wp), PARAMETER :: one = 1.0_wp
      REAL (wp), PARAMETER :: two = 2.0_wp
      REAL (wp), PARAMETER :: zero = 0.0_wp
      REAL (wp), PARAMETER :: mone = -one
      ! .. Local Scalars ..
      INTEGER :: i, ix, iy, j, k1, k2, level, num_iter, nx, nxy, ny
      REAL (wp) :: acc, hx, hx2, hy, hy2, mu, pi, resid_norm, rms_err, tpi2, &
       xj
      ! .. Local Arrays ..
      REAL (wp) :: a(maxsize,7), rhs(maxsize), sol(maxsize), u(maxsize)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) 'Example Program Results for nag_pde_ell_mg_ex03'

      pi = nag_pi(zero)
      tpi2 = two*pi*pi

      ! Do-loop on the level

      DO level = level_low, level_up
        nx = 2**level + 1
        ny = nx
        nxy = nx*ny
        acc = MAX(1.0E-6_wp,SQRT(nx*EPSILON(zero)))

        ! Set-up operator, right-hand side and
        ! initial guess

        hx = one/REAL(nx-1,kind=wp)
        hx2 = hx*hx
        hy = one/REAL(ny-1,kind=wp)
        hy2 = hy*hy
        mu = two/hx2 + two/hy2

        a(1:nxy,1) = mone/hy2
        a(1:nxy,2) = zero
        a(1:nxy,3) = mone/hx2
        a(1:nxy,4) = mu
        a(1:nxy,5) = mone/hx2
        a(1:nxy,6) = zero
        a(1:nxy,7) = mone/hy2

        u(1:nxy) = zero

        ! Exact solution "sol"
        DO j = 1, ny
          xj = REAL(j-1,kind=wp)*hy
          k2 = j*nx
          k1 = k2 + 1 - nx
          sol(k1:k2) = SIN(pi*xj)*SIN(pi*hx*(/(REAL(i-1,kind=wp),i=1,nx)/))
        END DO

        rhs(1:nxy) = tpi2*sol(1:nxy)

        ! Correction for the boundary conditions
        ! Horizontal boundaries
        DO i = 1, nx
          ! Boundary condition Y = 0
          ix = i
          u(ix) = sol(ix)
          rhs(ix) = mu*sol(ix)
          a(ix,1:7:2) = zero
          ! Boundary condition Y = 1
          ix = i + (ny-1)*nx
          u(ix) = sol(ix)
          rhs(ix) = mu*sol(ix)
          a(ix,1:7:2) = zero
        END DO
        ! Vertical boundaries
        DO j = 1, ny
          ! Boundary condition X = 0
          iy = (j-1)*nx + 1
          u(iy) = sol(iy)
          rhs(iy) = mu*sol(iy)
          a(iy,1:7:2) = zero
          ! Boundary condition X = 1
          iy = j*nx
          u(iy) = sol(iy)
          rhs(iy) = mu*sol(iy)
          a(iy,1:7:2) = zero
        END DO
        ! Solve the equation
        CALL nag_pde_ell_mg_sol(nx,ny,a(1:nxy,1:7),rhs(1:nxy),u(1:nxy), &
         resid_norm=resid_norm,acc=acc,num_iter=num_iter)
        WRITE (nag_std_out,*) ' '
        WRITE (nag_std_out,999) ' Level ', level
        WRITE (nag_std_out,999) ' Order of the matrix ', nxy
        WRITE (nag_std_out,998) ' Accuracy ', acc
        WRITE (nag_std_out,999) ' Number of iteration ', num_iter
        WRITE (nag_std_out,998) ' Residual norm ', resid_norm

        rms_err = SUM((u(:nxy)-sol(:nxy))**2)
        rms_err = SQRT(rms_err/REAL(nxy,kind=wp))

        WRITE (nag_std_out,998) ' RMS Error ', rms_err
      END DO

999   FORMAT (1X,A,I7)
998   FORMAT (1X,A,1P,E10.4)

    END PROGRAM nag_pde_ell_mg_ex03
