*     D03EDF Example Program Text
*     Mark 14 Revised. NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      DOUBLE PRECISION ALPHA
      PARAMETER        (ALPHA=1.7D0)
      INTEGER          LEVELS, NGX, NGY, LDA
      PARAMETER        (LEVELS=3,NGX=2**LEVELS+1,NGY=NGX,LDA=4*(NGX+1)
     +                 *(NGY+1)/3)
*     .. Local Scalars ..
      DOUBLE PRECISION ACC, HX, HY
      INTEGER          I, IFAIL, IOUT, IX, IY, J, K, MAXIT, NUMIT,
     +                 OUTCHN
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,7), RHS(LDA), U(LDA), UB(NGX*NGY), US(LDA)
*     .. External Subroutines ..
      EXTERNAL         D03EDF, X04ABF
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D03EDF Example Program Results'
      OUTCHN = NOUT
      WRITE (NOUT,*)
      ACC = 1.0D-4
      CALL X04ABF(1,OUTCHN)
      MAXIT = 15
*     ** Set IOUT.GE.2 to obtain intermediate output **
      IOUT = 0
      HX = 1.0D0/DBLE(NGX+1)
      HY = 1.0D0/DBLE(NGY+1)
*     Set up operator, right-hand side and initial guess for
*     step-lengths HX and HY
      DO 40 J = 1, NGY
         DO 20 I = 1, NGX
            K = (J-1)*NGX + I
            A(K,1) = 1.0D0 - 0.5D0*ALPHA
            A(K,2) = 0.5D0*ALPHA
            A(K,3) = 1.0D0 - 0.5D0*ALPHA
            A(K,4) = -4.0D0 + ALPHA
            A(K,5) = 1.0D0 - 0.5D0*ALPHA
            A(K,6) = 0.5D0*ALPHA
            A(K,7) = 1.0D0 - 0.5D0*ALPHA
            RHS(K) = -4.0D0*HX*HY
            UB(K) = 0.0D0
   20    CONTINUE
   40 CONTINUE
*     Correction for the boundary conditions
*     Horizontal boundaries --
      DO 60 I = 2, NGX - 1
*        Boundary condition on Y=0 -- U=0
         IX = I
         A(IX,1) = 0.0D0
         A(IX,2) = 0.0D0
*        Boundary condition on Y=1 -- U=0
         IX = I + (NGY-1)*NGX
         A(IX,6) = 0.0D0
         A(IX,7) = 0.0D0
   60 CONTINUE
*     Vertical boundaries --
      DO 80 J = 2, NGY - 1
*        Boundary condition on X=0 -- U=0
         IY = (J-1)*NGX + 1
         A(IY,3) = 0.0D0
         A(IY,6) = 0.0D0
*        Boundary condition on X=1 -- U=1
         IY = J*NGX
         RHS(IY) = RHS(IY) - A(IY,5) - A(IY,2)
         A(IY,2) = 0.0D0
         A(IY,5) = 0.0D0
   80 CONTINUE
*     Now the four corners --
*     Bottom left corner
      K = 1
      A(K,1) = 0.0D0
      A(K,2) = 0.0D0
      A(K,3) = 0.0D0
      A(K,6) = 0.0D0
*     Top left corner
      K = 1 + (NGY-1)*NGX
      A(K,3) = 0.0D0
      A(K,6) = 0.0D0
      A(K,7) = 0.0D0
*     Bottom right corner
*     Use average value at discontinuity ( = 0.5 )
      K = NGX
      RHS(K) = RHS(K) - A(K,2)*0.5D0 - A(K,5)
      A(K,1) = 0.0D0
      A(K,2) = 0.0D0
      A(K,5) = 0.0D0
*     Top right corner
      K = NGX*NGY
      RHS(K) = RHS(K) - A(K,2) - A(K,5)
      A(K,2) = 0.0D0
      A(K,5) = 0.0D0
      A(K,6) = 0.0D0
      A(K,7) = 0.0D0
*     Solve the equations
      IFAIL = 1
*
      CALL D03EDF(NGX,NGY,LDA,A,RHS,UB,MAXIT,ACC,US,U,IOUT,NUMIT,IFAIL)
*
      IF (IFAIL.EQ.0) THEN
         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 100 J = 1, NGY
            WRITE (NOUT,99995) J, (U(I+(J-1)*NGX),I=1,NGX)
  100    CONTINUE
      ELSE
         WRITE (NOUT,99994) IFAIL
      END IF
*
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,:)
99994 FORMAT (1X,' ** D03EDF returned with IFAIL = ',I5)
      END