
      SUBROUTINE INIT1 (NP, H, F, G, U, FR)
      IMPLICIT NONE
      INTEGER NP
      DOUBLE PRECISION F, G, H, FR (NP, NP), U (NP, NP)
      EXTERNAL F, G
C
C     COMPUTES INITIAL VALUES ON THE FINEST GRID
C     (THIS COULD BE PARALLELISED WITH FORALL IF `F' AND `G' ARE PURE).
C
        INTEGER I, J, N
        DOUBLE PRECISION H2, X, Y

        H2 = H * H
        N  = NP - 1
        DO J = 2, N
          Y = (J - 1) * H
          U (1, J) = G (0.0D0, Y)
          DO I = 2, N
            X = (I - 1) * H
            U  (I, J) = 0.0D0
            FR (I, J) = H2 * F (X, Y)
          ENDDO
          U (NP, J) = G (1.0D0, Y)
        ENDDO
        DO I = 1, NP
          X = (I - 1) * H
          U (I, 1)  = G (X, 0.0D0)
          U (I, NP) = G (X, 1.0D0)
        ENDDO
      END
C
C.......................................................................
C
      SUBROUTINE INT2A (NPC, NPF, UC, UF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION UC (NPC, NPC), UF (NPF, NPF)
C
C     COARSE-TO-FINE TRANSFER (BILINEAR INTERPOLATION) AND CORRECTION.
C     (THIS MODIFIES ODD SITES ONLY OF UF)
C
        INTEGER NC

        NC = NPC - 1
        UF (2 : 2 * NC : 2, 3 : 2 * NC - 1 : 2) =
     1      UF (2 : 2 * NC : 2, 3 : 2 * NC - 1 : 2) +
     2          0.5D0 * (UC (1 : NC, 2 : NC) + UC (2 : NC + 1, 2 : NC))
        UF (3 : 2 * NC - 1 : 2, 2 : 2 * NC : 2) =
     1      UF (3 : 2 * NC - 1 : 2, 2 : 2 * NC : 2) +
     2          0.5D0 * (UC (2 : NC, 1 : NC) + UC (2 : NC, 2 : NC + 1))
      END
C
C.......................................................................
C
      SUBROUTINE INT4 (NPC, NPF, UC, UF, FRF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION FRF (NPF, NPF), UC (NPC, NPC), UF (NPF, NPF)
C
C     4-TH ORDER FMG-INTERPOLATION
C     (THIS USES AUTOMATIC TEMPORARIES---THIS COULD BE AVOIDED
C     AT THE COST OF REORDERING SOME ARITHMETIC).
C
        INTEGER NC
        DOUBLE PRECISION HA (NPC, NPC), HN (NPC, NPC)
C
C       TRANSFER OF COARSE GRID VALUES TO THE FINE GRID
C
        NC = NPC - 1
        UF (3 : 2 * NC - 1 : 2, 3 : 2 * NC - 1 : 2) =
     1                                    UC (2 : NC, 2 : NC)
C
C       COMPUTATION OF THE REMAINING EVEN POINTS BY USING THE
C       ROTATED 5-POINT STAR
C
        HN ( : , 2 : NPC) =
     1      UC ( : , 1 : NPC - 1) + UC ( : , 2 : NPC)
        HA (2 : NPC, 2 : NPC) = HN (1 : NPC - 1, 2 : NPC)
        UF (2 : 2 * NPC - 2, 2 : 2 * NPC - 2) =
     1      0.25D0 * (2 * FRF (2 : 2 * NPC - 2, 2 : 2 * NPC - 2) +
     2           HA (2 : NPC, 2 : NPC) + HN (2 : NPC, 2 : NPC))
C
C       COMPUTATION OF THE ODD POINTS BY ONE HALF (ODD) RELAXATION STEP
C       **** WHAT?? -VE ITMAX DOESN'T DO THAT !! ****
C
        CALL RELAX (-1, .FALSE., NPF, UF, FRF)
      END
C
C.......................................................................
C
      SUBROUTINE PUTZB (NPC, UC)
      IMPLICIT NONE
      INTEGER NPC
      DOUBLE PRECISION UC (NPC, NPC)
C
C     PUTS ZERO TO BOUNDARY OF COARSER GRIDS
C
        INTEGER NC

        NC = NPC - 1
        UC (1,   2 : NC)  = 0.0D0
        UC (NPC, 2 : NC)  = 0.0D0
        UC (1 : NPC, 1)   = 0.0D0
        UC (1 : NPC, NPC) = 0.0D0
      END
C
C.......................................................................
C
      SUBROUTINE RELAX (ITMAX, IZ, NPF, UF, FF)
      IMPLICIT NONE
      INTEGER ITMAX, NPF
      LOGICAL IZ
      DOUBLE PRECISION FF (NPF, NPF), UF (NPF, NPF)
C
C     RED BLACK RELAXATION ON GRID NPF
C
        INTEGER IT, NXF, NYF, NXPF, NYPF, NXF2, NYF2

        NXPF = NPF
        NYPF = NPF
        NXF  = NXPF - 1
        NYF  = NYPF - 1
        NXF2 = NXF - 1
        NYF2 = NYF - 1
        DO IT = 1, ITMAX
          IF (IZ .AND. (IT .LE. 1)) THEN
C
C           ZERO STARTING VALUES
C           RELAXATION OF RED POINTS
C
            UF (2 : NXF : 2, 2 : NYF : 2) =
     1          FF (2 : NXF : 2, 2 : NYF : 2) * 0.25
            UF (3 : NXF : 2, 3 : NYF : 2) =
     1          FF (3 : NXF : 2, 3 : NYF : 2) * 0.25
          ELSE
C
C           NON-ZERO STARTING VALUES
C           RELAXATION OF RED POINTS
C
            UF (2 : NXF2 : 2, 2 : NYF2 : 2) =
     1          0.25 * (FF (2 : NXF2 : 2, 2 : NYF2 : 2) +
     2                 (UF (1 : NXF2 - 1 : 2, 2 : NYF2 : 2) +
     3                 (UF (2 : NXF2 : 2, 1 : NYF2 - 1 : 2) +
     4                 (UF (3 : NXF2 + 1 : 2, 2 : NYF2 : 2) +
     5                  UF (2 : NXF2 : 2, 3 : NYF2 + 1 : 2)))))

            UF (3 : NXF2 + 1 : 2, 3 : NYF2 + 1 : 2) =
     1          0.25 * (FF (3 : NXF2 + 1 : 2, 3 : NYF2 + 1 : 2) +
     2                 (UF (4 : NXF2 + 2 : 2, 3 : NYF2 + 1 : 2) +
     3                 (UF (3 : NXF2 + 1 : 2, 4 : NYF2 + 2 : 2) +
     4                 (UF (3 : NXF2 + 1 : 2, 2 : NYF2 : 2) +
     5                  UF (2 : NXF2 : 2, 3 : NYF2 + 1 : 2)))))

            UF (2 : NXF : 2, NYF) = 0.25 * (FF (2 : NXF : 2, NYF) +
     1          (UF (1 : NXF - 1 : 2, NYF) +
     2          (UF (3 : NXF + 1 : 2, NYF) +
     3          (UF (2 : NXF : 2, NYF - 1) +
     4           UF (2 : NXF : 2, NYF + 1)))))

            UF (NXF, 2 : NYF : 2) = 0.25 * (FF (NXF, 2 : NYF : 2) +
     1          (UF (NXF, 1 : NYF - 1 : 2) +
     2          (UF (NXF, 3 : NYF + 1 : 2) +
     3          (UF (NXF - 1, 2 : NYF : 2) +
     4           UF (NXF + 1, 2 : NYF : 2)))))
          ENDIF
C
C         RELAXATION OF BLACK POINTS
C
          UF (3 : NXF2 : 2, 2 : NYF2 : 2) =
     1        0.25 * (FF (3 : NXF2 : 2, 2 : NYF2 : 2) +
     2               (UF (3 : NXF2 : 2, 1 : NYF2 - 1 : 2) +
     3               (UF (4 : NXF2 + 1 : 2, 2 : NYF2 : 2) +
     4               (UF (2 : NXF2 - 1 : 2, 2 : NYF2 : 2) +
     5                UF (3 : NXF2 : 2, 3 : NYF2 + 1 : 2)))))

          UF (2 : NXF2 - 1 : 2, 3 : NYF2 + 1 : 2) =
     1        0.25 * (FF (2 : NXF2 - 1 : 2, 3 : NYF2 + 1 : 2) +
     2               (UF (1 : NXF2 - 2 : 2, 3 : NYF2 + 1 : 2) +
     3               (UF (2 : NXF2 - 1 : 2, 4 : NYF2 + 2 : 2) +
     4               (UF (2 : NXF2 - 1 : 2, 2 : NYF2 : 2) +
     5                UF (3 : NXF2 : 2, 3 : NYF2 + 1 : 2)))))

          UF (3 : NXF : 2, NYF) = 0.25 * (FF (3 : NXF : 2, NYF) +
     1        (UF (2 : NXF - 1 : 2, NYF) +
     2        (UF (4 : NXF + 1 : 2, NYF) +
     3        (UF (3 : NXF : 2, NYF - 1) +
     4         UF (3 : NXF : 2, NYF + 1)))))

          UF (NXF, 3 : NYF : 2) = 0.25 * (FF (NXF, 3 : NYF : 2) +
     1        (UF (NXF, 2 : NYF - 1 : 2) +
     2        (UF (NXF, 4 : NYF + 1 : 2) +
     3        (UF (NXF - 1, 3 : NYF : 2) +
     4         UF (NXF + 1, 3 : NYF : 2)))))
        ENDDO
C
      END
C
C.......................................................................
C
      SUBROUTINE RESTR (NPC, NPF, FC, UF, FF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION FC (NPC, NPC), FF (NPF, NPF), UF (NPF, NPF)
C
C     COMPUTATION OF THE DEFECT AND FINE-TO-COARSE TRANSFER
C     (HALF-INJECTION)
C
        INTEGER NC

        NC = NPC - 1
        FC (2 : NC, 2 : NC) = 
     1      2 * (FF (3 : 2 * NC - 1 : 2, 3 : 2 * NC - 1 : 2) -
     2      4.0D0 * UF (3 : 2 * NC - 1 : 2, 3 : 2 * NC - 1 : 2) +
     3              UF (3 : 2 * NC - 1 : 2, 2 : 2 * NC - 2 : 2) +
     4              UF (2 : 2 * NC - 2 : 2, 3 : 2 * NC - 1 : 2) +
     5              UF (4 : 2 * NC : 2, 3 : 2 * NC - 1 : 2) +
     6              UF (3 : 2 * NC - 1 : 2, 4 : 2 * NC : 2))
      END
C
C.......................................................................
C
      SUBROUTINE TRANS (NPC, NPF, UC, FC, UF, FF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION FC (NPC, NPC), FF (NPF, NPF),
     1                 UC (NPC, NPC), UF (NPF, NPF)
C
C     TRANSFER OF F AND U FROM GRID NPF TO NPC
C
        INTEGER NC

        NC = NPC - 1
        UC (2 : NC, 2 : NC) = 
     1      UF (3 : 2 * NC - 1 : 2, 3 : 2 * NC - 1 : 2)
        FC (2 : NC, 2 : NC) = 
     1      4.0D0 * FF (3 : 2 * NC - 1 : 2, 3 : 2 * NC - 1 : 2)
        UC (1,   2 : NC)  = UF (1,   3 : 2 * NC - 1 : 2)
        UC (NPC, 2 : NC)  = UF (NPF, 3 : 2 * NC - 1 : 2)
        UC (1 : NPC, 1)   = UF (1 : NPF : 2, 1)
        UC (1 : NPC, NPC) = UF (1 : NPF : 2, NPC)
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION DEFMX (UC, FC, NP)
      IMPLICIT NONE
      INTEGER NP
      DOUBLE PRECISION UC (NP, NP), FC (NP, NP)
C
C     COMPUTES THE MAXIMUM NORM OF THE DEFECT
C     (PARALLEL VERSION DOESN'T WORK FOR SOME REASON---WITH NAG COMPILER,
C     ANYWAY---BUT THIS ROUTINE ISN'T TIMED, SO DON'T WORRY).
C
        INTEGER N, I, J
        DOUBLE PRECISION D, HSQR

        N = NP - 1
        HSQR  = DBLE (FLOAT (N * N))
C        DEFMX = MAXVAL (DABS
C     1      ((FC (2 : N, 2 : N) - 4.0D0 * UC (2 : N, 2 : N) +
C     2        UC (1 : N - 1, 2 : N) + UC (2 : N, 1 : N - 1) +
C     3        UC (2 : N, 3 : N + 1) + UC (2 : N + 1, 2 : N)) * HSQR))
        DEFMX = 0.0D0
        DO J = 2, N
          DO I = 2, N
            D = (FC (I, J) - 4.0D0 * UC (I, J) +
     1           UC (I - 1, J) + UC (I, J - 1) +
     2           UC (I, J + 1) + UC (I + 1, J)) * HSQR
            DEFMX = DMAX1 (DEFMX, DABS (D))
          ENDDO
        ENDDO
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION DIFMX (UC, NP, SOL)
      IMPLICIT NONE
      INTEGER NP
      DOUBLE PRECISION UC (NP, NP), SOL
      EXTERNAL SOL
C
C     COMPUTES THE MAXIMUM NORM OF THE DIFFERENCE BETWEEN
C     SOL (=SOLUTION OF THE BVP) AND THE VALUES IN UC.
C     (THIS COULD BE PARALLELISED WITH FORALL IF `SOL' WAS PURE)
C
        INTEGER I, J
        DOUBLE PRECISION H, X, Y

        H     = 1.0D0 / (NP - 1)
        DIFMX = 0.0D0
        DO J = 1, NP
          Y = (J - 1) * H
          DO I = 1, NP
            X = (I - 1) * H
            DIFMX = DMAX1 (DIFMX, DABS (SOL (X, Y) - UC (I, J)))
          ENDDO
        ENDDO
      END

