      PROGRAM PDE2
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     MAIN PROGRAM FOR MG BENCHMARK - SIMD/VECTOR VERSION
C
C     SUPPLIED BY   MAX LEMKE
C                   SUPRENUM GMBH
C                   HOHE STRASSE 73
C                   D-5300 BONN 1
C
C     VERSION FEBRUARY 15, 1989
C
C     REFERENCES:
C               - K. STUEBEN, U. TROTTENBERG: MULTIGRID METHODS:
C                 FUNDAMENTAL ALGORITHMS, MODEL PROBLEM ANALYSIS
C                 AND APPLICATIONS, DECEMBER 84, GMD-STUDIEN NR. 96,
C                 GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG MBH,
C                 BONN
C
C               - M. LEMKE: EXPERIMENTS WITH A VECTORIZED MULTIGRID
C                 POISSON SOLVER ON THE CDC CYBER 205, CRAY X-MP AND
C                 FUJITSU VP 200;
C                 NOVEMBER 85, ARBEITSPAPIERE DER GMD, NR. 179,
C                 GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG MBH,
C                 BONN
C
C               - M. LEMKE: ERFAHRUNGEN MIT MEHRGITTERVERFAHREN FUER
C                 HELMHOLTZ-AEHNLICHE PROBLEME AUF VEKTORRECHNERN
C                 UND MEHRPROZESSOR-VEKTORRECHNERN
C                 DECEMBER 87, ARBEITSPAPIERE DER GMD, NR. 278,
C                 GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG MBH,
C                 BONN
C
C     AUTHORS: FIRST VERSION:
C                 K. WITSCH, H. FOERSTER (1984, GMD)
C              VECTORIZED (AND MULTITASKING) VERSION:
C                 M. LEMKE (1984 - 1987, GMD)
C              GENESIS BENCHMARK VERSION:
C                 M. LEMKE (1989, SUPRENUM GMBH)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     THE TIMER (DOUBLE PRECISION IN SECONDS) IS CALLED IN MG00D.
C     THE MACHINE DEPENDENT TIMING ROUTINE HAS TO BE INCLUDED.
C
C     ATTENTION: ONLY THE MG KERNEL IS VECTORIZED. SOME OTHER ROUTINES
C                ARE PURELY SEQUENTIAL. THEY NEED NOT BE VECTORIZED
C                BECAUSE THEY ARE NOT INCLUDED IN TIMING.
C
C                THE FMG OPTION IS NOT AVAILABLE IN THIS VERSION
C
C     FOR NUMERICAL VERIFICATION:
C        THE AVERAGE CONVERGENCE RATE SHOLD BE SOMEWHERE IN BETWEEN
C        0.02 AND 0.07.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      DOUBLE PRECISION F, G, DIFMX, DEFMX, DIF, DEF
      INTEGER IDIM, IER, M, N, NP, NY1, NY2, NCYCLE, NFMG, IGAM, INITF
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     DEFINITION OF FINEST GRIDSIZE: (2**MMAX+1)*(2**MMAX+1)
C
C     RESULTS ARE OBTAINED FOR THE FINEST GRID RANGING FROM 4 TO MMAX.
C     MMAX SHOULD BE CHOSEN AS LARGE AS POSSIBLE, DEPENDING ON MEMORY
C     CAPACITIES OF THE MACHINE.
C
      INCLUDE 'pde2.inc'
C
C---> the output unit number
      PARAMETER(NW1=11)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      PARAMETER (IDIM=(2**MMAX+1)*(2**MMAX+1)*3)
      DOUBLE PRECISION W(IDIM), DEF1,DEF2
      COMMON/DEMOR/ W
      DOUBLE PRECISION SEC2
      COMMON /TIME/ SEC2
      COMMON /DEFECT/ DEF1,DEF2
      EXTERNAL F, G
      SEC2=1.
      OPEN(NW1,FILE='result')
      REWIND(NW1)
      CALL HEADER(NW1)
      WRITE(NW1,91000)
91000 FORMAT(' ************* MULTIGRID BENCHMARK ******************* '/
     *     )
C
C---> LOOP OVER THE FINEST GRIDS IN ORDER TO GET RESULTS FOR DIFFERENT
C     PROBLEM SIZES
C
      NCYCLE = NITER

      DO 11 M=3,MMAX
C
        NY1 = 2
        NY2 = 1
        NFMG = 0
        IGAM = 1
C
        CALL MG00D(M,NY1,NY2,NFMG,NCYCLE,IGAM,F,G,W,IDIM,INITF,IER)
C
        IF (IER.EQ.0) GOTO 10
        WRITE (NW1,9000) IER, IDIM
        STOP
 10     N = 2**M
        NP = N+1
        DIF = DIFMX(W,NP,G)
        DEF = DEFMX(W,W(INITF+1),NP)
        CONVR=(DEF/DEF1)**(0.1)
        WRITE (NW1,9100)NP,NP**2, DIF, DEF, CONVR
C
        NOP = 1
        DO 123 LEV = 2,M
            NOP = NOP + (2**LEV-1)*(2**LEV-1)*(2.75+9.)
     *                + (2**(LEV-1)-1)*(2**(LEV-1)-1)*7
     *                + (2**(LEV-1)-1)*(2**(LEV-1)-1)*6
123     CONTINUE
            NOP = NOP + (2**M-1)*(2**M-1)*1.75
C
        RMFLOP = (NOP/1000000.)/SEC2*NCYCLE
        WRITE(NW1,900) NCYCLE, SEC2, SEC2/NCYCLE, NOP, RMFLOP
        WRITE(NW1,91100)
C
        NCYCLE = MAX(NCYCLE/4,10)
C
11    CONTINUE
C
      STOP
91100 FORMAT(' +++++++++++++++++++++++++++++++++++++++++++++++++++++ '/
     *      )
 9000 FORMAT (' *** ERROR, IER=', I3,' IDIM=', I6, ' ***')
 9100 FORMAT (/' NUMBER OF POINTS PER DIRECTION = ', I12 /
     *         ' GRIDSIZE OF SQUARE GRID =        ', I12/
     *         ' MAXIMUM NORM OF THE ERROR  =     ', D12.4 /
     *         ' MAXIMUM NORM OF THE DEFECT =     ', D12.4/
     *         ' AVERAGE CONV. RATE             = ', D12.4/)
900     FORMAT(/' TIME FOR ', I4, ' CYCLES (SEC.):', E15.8/
     *           ' TIME FOR     1 CYCLE (SEC.):', E15.8/
     *           ' NUMBER OF OPERATIONS:       ',I15/
     *           ' PERFORMANCE IN MFLOP/S:     ',F15.5/)
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION DEFMX(UC, FC, NP)
C
C       COMPUTES THE MAXIMUM NORM OF THE DEFECT
C
      INTEGER N, NP, I, J
      DOUBLE PRECISION UC(NP,NP), FC(NP,NP), D, HSQR
C
      N = NP - 1
      HSQR = DBLE(FLOAT(N*N))
      DEFMX = 0.0D0
      DO 20 J=2,N
        DO 10 I=2,N
          D = (FC(I,J)-4.0D0*UC(I,J)+UC(I-1,J)+UC(I,J-1)+UC(I,J+1)
     *         +UC(I+1,J))*HSQR
          DEFMX = DMAX1(DEFMX,DABS(D))
   10   CONTINUE
   20 CONTINUE
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION DIFMX(UC, NP, SOL)
C
C       COMPUTES THE MAXIMUM NORM OF THE DIFFERENCE BETWEEN
C       SOL (=SOLUTION OF THE BVP) AND THE VALUES IN UC
C
      INTEGER NP, I, J
      DOUBLE PRECISION UC(NP,NP), SOL, DFLOAT, H, X, Y
      DFLOAT(K) = DBLE(FLOAT(K))
C
      H = 1.0D0/DFLOAT(NP-1)
      DIFMX = 0.0D0
      DO 20 J=1,NP
        Y = DFLOAT(J-1)*H
        DO 10 I=1,NP
          X = DFLOAT(I-1)*H
          DIFMX = DMAX1(DIFMX,DABS(SOL(X,Y)-UC(I,J)))
   10   CONTINUE
   20 CONTINUE
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION F(X,Y)
      DOUBLE PRECISION X, Y
      F = 10.0D0*SIN(3.0D0*X+Y)
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION G(X,Y)
      DOUBLE PRECISION X, Y
      G = SIN(3.0D0*X+Y)
      END
C
C.......................................................................
C
      SUBROUTINE INIT1(NP, H, F, G, U, FR)
C
C        COMPUTES INITIAL VALUES ON THE FINEST GRID
C
      INTEGER I, J, N, NP
      DOUBLE PRECISION FR(NP,NP), U(NP,NP), F, G, DFLOAT, H, H2, X, Y
      DFLOAT(K) = DBLE(FLOAT(K))
C
      H2 = H*H
      N = NP - 1
      DO 20 J=2,N
        Y = DFLOAT(J-1)*H
        U(1,J) = G(0.0D0,Y)
        DO 10 I=2,N
          X = DFLOAT(I-1)*H
          U(I,J) = 0.0D0
          FR(I,J) = H2*F(X,Y)
   10   CONTINUE
        U(NP,J) = G(1.0D0,Y)
   20 CONTINUE
      DO 30 I=1,NP
        X = DFLOAT(I-1)*H
        U(I,1) = G(X,0.0D0)
        U(I,NP) = G(X,1.0D0)
   30 CONTINUE
      END
C
C.......................................................................
C
      SUBROUTINE INIT2(M, LIN, NPK, ID, U, FR, IDIM)
C
C        COMPUTES INITIAL VALUES ON THE COARSER GRIDS
C
      INTEGER IDIM, LIN, M, IDC, IDF, L, LH, ID(125), NPK(65)
      DOUBLE PRECISION FR(IDIM), U(IDIM)
C
      IF (LIN.EQ.1) GO TO 20
C
C        TRANSFER OF F AND U TO COARSER GRIDS
C
      DO 10 L=2,LIN
        IDC = ID(L)
        IDF = ID(L-1)
        CALL TRANS(NPK(L), NPK(L-1), U(IDC), FR(IDC), U(IDF), FR(IDF))
   10 CONTINUE
C
C        PUT ZERO TO BOUNDARY VALUES FOR MULTIGRID CORRECTIONS
C
   20 IF (LIN.EQ.M) GO TO 40
      LH = LIN + 1
      DO 30 L=LH,M
        IDC = ID(L)
        CALL PUTZB(NPK(L), U(IDC))
   30 CONTINUE
   40 RETURN
      END
C
C.......................................................................
C
      SUBROUTINE INT2A(NPC, NPF, UC, UF)
C
C       COARSE-TO-FINE TRANSFER (BILINEAR INTERPOLATION) AND CORRECTION
C
      INTEGER NC, NPC, NPF, I, IF, J, JF
      DOUBLE PRECISION UC(NPC,NPC), UF(NPF,NPF)
C
      NC = NPC - 1
      DO 20 J=2,NC
        JF = J + J - 1
        DO 10 I=1,NC
          IF = I + I
          UF(IF,JF) = UF(IF,JF) + 0.5D0*(UC(I,J)+UC(I+1,J))
   10   CONTINUE
   20 CONTINUE
      DO 40 J=1,NC
        JF = J + J
        DO 30 I=1,NC
          IF = I + I - 1
          UF(IF,JF) = UF(IF,JF) + 0.5D0*(UC(I,J)+UC(I,J+1))
   30   CONTINUE
   40 CONTINUE
      END
C
C.......................................................................
C
      SUBROUTINE INT4(NPC, NPF, UC, UF, FRF)
C
C       4-TH ORDER FMG-INTERPOLATION
C
      INTEGER NC, NPC, NPF, I, IF, J, JF, JM
      DOUBLE PRECISION FRF(NPF,NPF), UC(NPC,NPC), UF(NPF,NPF), HA, HN
C
C     TRANSFER OF COARSE GRID VALUES TO THE FINE GRID
C
      NC = NPC - 1
      JF = 3
      DO 20 J=2,NC
        IF = 3
        DO 10 I=2,NC
          UF(IF,JF) = UC(I,J)
          IF = IF + 2
   10   CONTINUE
        JF = JF + 2
   20 CONTINUE
C
C     COMPUTATION OF THE REMAINING EVEN POINTS BY USING THE
C     ROTATED 5-POINT STAR
C
      JF = 2
      JM = 1
      DO 40 J=2,NPC
        IF = 2
        HA = UC(1,JM) + UC(1,J)
        DO 30 I=2,NPC
          HN = UC(I,JM) + UC(I,J)
          UF(IF,JF) = 0.25D0*(FRF(IF,JF)+FRF(IF,JF)+HA+HN)
          HA = HN
          IF = IF + 2
   30   CONTINUE
        JM = JM + 1
        JF = JF + 2
   40 CONTINUE
C
C     COMPUTATION OF THE ODD POINTS BY ONE HALF (ODD) RELAXATION STEP
C
      CALL RELAX(-1, 1, NPF, UF, FRF, FRF)
      END
C
C***********************************************************************
C
      SUBROUTINE MG00D(M, NY1, NY2, NFMG, NCYCLE, IGAM, F, G, W, IDIM,
     *                 INITF, IER)
C
C     MULTIGRID MODULE FOR THE FAST SOLUTION OF POISSON'S EQUATION
C     WITH DIRICHLET BOUNDARY CONDITIONS ON THE UNIT SQUARE
C
C          DIFFERENTIAL EQUATION :   -DELTA U(X,Y) = F(X,Y)
C          BOUNDARY CONDITION :             U(X,Y) = G(X,Y)
C
C     CYCLE STRUCTURE:
C
C                COARSENING: H, 2H, 4H, ...
C                DIFFERENCE OPERATORS: USUAL 5-POINT STARS ON ALL GRIDS
C                RELAXATION: RED-BLACK
C                FINE-TO-COARSE: HALF INJECTION
C                COARSE-TO-FINE: BILINEAR INTERPOLATION
C
C     FULL MULTIGRID MODE:
C
C                FULL MULTIGRID INTERPOLATION USES GRID EQUATION
C                (4-TH ORDER)
C
C     INPUT:
C
C        M       NUMBER OF GRIDS (0 < M < 11). FOR GIVEN M:
C                NUMBER OF P0INTS ON THE FINEST GRID (INCLUDING BOUNDARY
C                POINTS) = 2**M + 1 IN BOTH DIRECTIONS
C        NY1     NUMBER OF RELAXATIONS BEFORE COARSE-GRID CORRECTION
C                (NY1 > 0)
C        NY2     NUMBER OF RELAXATIONS AFTER COARSE-GRID CORRECTION
C                (NY2 > 0)
C        NCYCLE  NUMBER OF MULTIGRID ITERATIONS (NCYCLE > 0)
C        NFMG    .EQ. 0 : NCYCLE MULTIGRID ITERATIONS ARE PERFORMED
C                .NE. 0 : FMG-VERSION IS PERFORMED PLUS NCYCLE-1
C                         ADDITIONAL MG ITERATIONS AFTERWARDS
C        IGAM    TYPE OF CYCLING (IGAM > 0).
C                E.G.: IGAM=1 FOR V-CYCLES, IGAM=2 FOR W-CYCLES
C        W       DOUBLE PRECISION WORK ARRAY OF DIMENSION IDIM
C        IDIM    DIMENSION OF W. APPROXIMATELY: IDIM > 2.8*4**M
C
C     EXTERNALS:
C
C        F       DOUBLE PRECISION FUNCTION F(X,Y), RIGHT HAND SIDE OF
C                THE DIFFERENTIAL EQUATION
C        G       DOUBLE PRECISION FUNCTION G(X,Y), BOUNDARY VALUES
C
C     REMARK:    GRID #1 AND #M ARE THE FINEST AND COARSEST GRID USED,
C                RESPECTIVELY
C
C     OUTPUT:
C
C        IER     ERROR INDICATOR
C                = 0  NO ERRORS
C                = 1  INSUFFICIENT MEMORY , I.E. IDIM TOO SMALL.
C                     IN THIS CASE IDIM IS USED AS OUTPUT PARAMETER
C                     TO SHOW THE MINIMAL DIMENSION
C                = 2  M, NY1, NY2, NCYCLE OR IGAM WRONG
C        IDIM    ONLY IN CASE IER=1 : MINIMAL LENGTH OF W
C        INITF   CF. DESCRIPTION OF W
C        W       W CONTAINS THE DISCRETE APPROXIMATION TO THE GIVEN
C                BOUNDARY VALUE PROBLEM ON THE FINEST GRID. THE GRID
C                VALUES ARE STORED ROWWISE FROM LEFT TO RIGHT AND
C                FROM BOTTOM TO TOP.
C                I.E. THE GRID VALUE CORRESPONDING TO THE GRID POINT
C                (XI,YJ) :
C
C                     XI = (I-1)*H,  YJ = (J-1)*H   (H=1/N, N=2**M)
C
C                IS STORED AT
C
C                        W((J-1)*(N+1)+I)     (0 < I,J < N+2)
C
C                THE CORRESPONDING VALUES OF THE RIGHT HAND SIDE ARE
C                STORED IN THE SAME MANNER AT
C
C                        W((J-1)*(N+1)+I+INITF).
C
C                THE REMAINING STORAGE CONTAINS COARSE GRID VALUES.
C
C     REMARK:    THIS IS A SPECIALIZED PROGRAM FROM THE MG00 PROGRAM
C                COLLECTION. SEE: FOERSTER, WITSCH IN
C                'ON EFFICIENT MULTIGRID SOFTWARE FOR ELLIPTIC
C                PROBLEMS ON RECTANGULAR DOMAINS', MATH. COMPUT.
C                SIMULATION XXIII, PP. 293-298, 1981
C
C
      INTEGER IDC, IDF, IDR, IT, ITM, INITF, L, L1, LEV, LIN, NP,
     * NY2, IGAM, IDIM, IER, NCYCLE, NFMG, M, ID(125), NPK(65)
      DOUBLE PRECISION  SEC1, SEC2, DEF1, DEF2
      COMMON /TIME/ SEC2
      COMMON /DEFECT/ DEF1,DEF2
      DOUBLE PRECISION W(   1), HK(50), DFLOAT
      DOUBLE PRECISION F,G, DEFMX
      EXTERNAL F, G, DEFMX
      DFLOAT(L1)=DBLE(FLOAT(L1))
      IER = 0
C
C        CHECK OF ERROR CONDITIONS (WRONG PARAMETERS)
C
      IF (M.LE.0 .OR. M.GE.110 .OR. NY1.LE.0 .OR. NY2.LE.0 .OR.
     *    NCYCLE.LE.0 .OR. IGAM.LE.0 ) IER=2
      IF (IER.GT.0) RETURN
C
C        DETERMINATION OF COARSER GRIDS
C
      NP = 2**M + 1
      HK(1) = 1.0D0/DFLOAT(NP-1)
      ID(1) = 1
      ID(2) = NP*NP + 1
      NPK(1) = NP
CVD$L NOVECTOR
      DO 10 L=2,M
        NPK(L) = (NPK(L-1)+1)/2
        HK(L) = HK(L-1) + HK(L-1)
        ID(L+1) = ID(L) + NPK(L)**2
   10 CONTINUE
C
C        CHECK OF ERROR CONDITIONS (INSUFFICIENT MEMORY)
C
      INITF = ID(M+1)
      IF (2*INITF.GT.IDIM) THEN
        IER = 1
        IDIM = 2*INITF
        RETURN 
      ENDIF
      LIN = 1
      IF (NFMG.NE.0) LIN = M
C
C        SET UP ALL GRID VALUES NEEDED
C
      CALL INIT1(NP, HK(1), F, G, W(1), W(INITF+1))
      CALL INIT2(M, LIN, NPK, ID, W(1), W(INITF+1), INITF)
C
C        (FULL) MULTIGRID PROCEDURE
C
      DO 40 L=1,LIN
        LEV = LIN - L + 1
        ITM = 1
        IF (LEV.EQ.1) ITM = NCYCLE
        DEF1= DEFMX(W,W(INITF+1),NP)
C----> TIMING
        CALL TIMER(SEC1)
        DO 30 IT=1,ITM
        CALL MGI(LEV, M, NY1, NY2, IGAM, NPK, ID, W, W(INITF+1), INITF)
   30   CONTINUE
        CALL TIMER(SEC2)
        SEC2 = SEC2 - SEC1
C----> TIMING
        IF (L.LT.LIN) THEN
          IDF = ID(LEV-1)
          IDC = ID(LEV)
          IDR = IDF+INITF
          CALL INT4(NPK(LEV), NPK(LEV-1), W(IDC), W(IDF), W(IDR))
          CALL PUTZB(NPK(LEV), W(IDC))
        ENDIF
   40 CONTINUE

      END
C
C.......................................................................
C
      SUBROUTINE MGI(LEV, M, NY1, NY2, IGAM, NPK, ID, U, FR, IDIM)
C
C        ONE MULTIGRID ITERATION STEP (ON ACTUAL FINEST GRID LEV)
C
      INTEGER IDIM, LEV, M, NY1, NY2, ID(125), NPK(65), IDC, IDF, IZ,
     *        K, IGAM, ICGAM(10)
      DOUBLE PRECISION FR(IDIM), U(IDIM)
C
      IZ = 1
      DO 10 K=LEV,M
        ICGAM(K) = 0
   10 CONTINUE
      K = LEV
      IF (K.EQ.M) GO TO 30
   20   IDF = ID(K)
C
C        RELAXATIONS BEFORE COARSE-GRID CORRECTION
C
        IZ = 1
        IF (K.GT.LEV .AND. ICGAM(K).EQ.0) IZ = 0
        CALL RELAX(NY1, IZ, NPK(K), U(IDF), FR(IDF), FR(1))
        ICGAM(K) = ICGAM(K) + 1
C
C        RESIDUAL TRANSFER TO NEXT COARSER GRID
C
        IDC = ID(K+1)
        CALL RESTR(NPK(K+1), NPK(K), FR(IDC), U(IDF), FR(IDF))
        K = K + 1
        IF (K.LT.M) GOTO 20
      IZ = 0
   30 IDC = ID(M)
C
C        EXACT SOLUTION ON COARSEST GRID
C
      CALL RELAX(1,IZ,NPK(M),U(IDC),FR(IDC),FR(1))
      IF (K.EQ.LEV) GOTO 50
   40   K = K - 1
        IDF = ID(K)
        IDC = ID(K+1)
C
C        LINEAR INTERPOLATION TO NEXT FINER GRID
C
        CALL INT2A(NPK(K+1), NPK(K), U(IDC), U(IDF))
C
C        RELAXATION AFTER CGC
C
        CALL RELAX(NY2, 1, NPK(K), U(IDF), FR(IDF), FR(1))
      IF (K.EQ.LEV) GOTO 50
      IF (ICGAM(K).LT.IGAM) GOTO 20
      ICGAM(K) = 0
      GOTO 40
   50 RETURN
      END
C
C.......................................................................
C
      SUBROUTINE PUTZB(NPC, UC)
C
C       PUTS ZERO TO BOUNDARY OF COARSER GRIDS
C
      INTEGER NC, NPC, I, J
      DOUBLE PRECISION UC(NPC,NPC)
C
      NC = NPC - 1
      DO 10 J=2,NC
        UC(1,J) = 0.0D0
        UC(NPC,J) = 0.0D0
   10 CONTINUE
      DO 20 I=1,NPC
        UC(I,1) = 0.0D0
        UC(I,NPC) = 0.0D0
   20 CONTINUE
      END
C
C.......................................................................
C
      SUBROUTINE RELAX(ITMAX, IZ, NPF,UF,FF, W)
C
C     RED BLACK RELAXATION ON GRID NPF
C
      DOUBLE PRECISION FF(NPF,NPF), UF(NPF,NPF), W(1)
      INTEGER I, IS, IT, J, NXF, NYF
C     ..
      NXPF = NPF
      NYPF = NPF
      NXF = NXPF - 1
      NYF = NYPF - 1
      NXF2 = NXF - 1
      NYF2 = NYF - 1
      DO 200 IT=1,ITMAX
        IS = 2
        IF (IZ.NE.0 .OR. IT.GT.1) GO TO 30
C
C       ZERO STARTING VALUES
C       RELAXATION OF RED POINTS
C
        DO 20 J=2,NYF,2
          DO 10 I=2,NXF,2
            UF(I,J) = FF(I,J)*0.25
   10     CONTINUE
   20   CONTINUE
        DO 21 J=3,NYF,2
          DO 11 I=3,NXF,2
            UF(I,J) = FF(I,J)*0.25
   11     CONTINUE
   21   CONTINUE
        GO TO 70
   30   CONTINUE
C
C       NON-ZERO STARTING VALUES
C       RELAXATION OF RED POINTS
C
CDIR$ IVDEP
        DO 60 J=2,NYF2,2
          DO 50 I=IS,NXF2,2
            UF(I  ,J  ) = 0.25*(FF(I,J)+(UF(I-1,J  )+(UF(I  ,J-1)+
     +                         (UF(I+1,J  )+ UF(I  ,J+1)))))
            UF(I+1,J+1) = 0.25*(FF(I+1,J+1)+(UF(I+2,J+1)+(UF(I+1,J+2)+
     +                         (UF(I+1,J  )+ UF(I  ,J+1)))))
   50     CONTINUE
   60   CONTINUE
CDIR$ IVDEP
        DO 150 I=2,NXF,2
            UF(I  ,NYF) = 0.25*(FF(I  ,NYF)+(UF(I-1,NYF)+(UF(I+1,NYF)+
     +                         (UF(I  ,NYF-1)+ UF(I  ,NYF+1)))))
  150   CONTINUE
CDIR$ IVDEP
        DO 160 J=2,NYF,2
            UF(NXF,J  ) = 0.25*(FF(NXF,J  )+(UF(NXF,J-1)+(UF(NXF,J+1)+
     +                         (UF(NXF-1,J  )+ UF(NXF+1,J  )))))
  160   CONTINUE
C
   70   IS = 3
C
C       RELAXATION OF BLACK POINTS
C
CDIR$ IVDEP
        DO 100 J=2,NYF2,2
          DO 90 I=IS,NXF2,2
            UF(I  ,J  ) = 0.25*(FF(I,J)+(UF(I  ,J-1)+(UF(I+1,J  )+
     +                         (UF(I-1,J  )+ UF(I  ,J+1)))))
            UF(I-1,J+1) = 0.25*(FF(I-1,J+1)+(UF(I-2,J+1)+(UF(I-1,J+2)+
     +                         (UF(I-1,J  )+ UF(I  ,J+1)))))
   90     CONTINUE
  100   CONTINUE
CDIR$ IVDEP
        DO 180 I=3,NXF,2
            UF(I  ,NYF) = 0.25*(FF(I  ,NYF)+(UF(I-1,NYF)+(UF(I+1,NYF)+
     +                         (UF(I  ,NYF-1)+ UF(I  ,NYF+1)))))
  180   CONTINUE
CDIR$ IVDEP
        DO 190 J=3,NYF,2
            UF(NXF,J  ) = 0.25*(FF(NXF,J  )+(UF(NXF,J-1)+(UF(NXF,J+1)+
     +                         (UF(NXF-1,J  )+ UF(NXF+1,J  )))))
  190   CONTINUE
  110 CONTINUE
  200 CONTINUE
C
      END
C
C.......................................................................
C
      SUBROUTINE RESTR(NPC, NPF, FC, UF, FF)
C
C       COMPUTATION OF THE DEFECT AND FINE-TO-COARSE TRANSFER
C       (HALF-INJECTION)
C
      INTEGER NC, NPC, NPF, I, IF, J, JF
      DOUBLE PRECISION FC(NPC,NPC), FF(NPF,NPF), UF(NPF,NPF)
C
      NC = NPC - 1
      DO 20 J=2,NC
        JF = J + J - 1
        DO 10 I=2,NC
          IF = I + I - 1
          FC(I,J) = 2 * ( FF(IF,JF) -
     *        4.0D0*UF(IF,JF) + UF(IF,JF-1) + UF(IF-1,JF)
     *                        + UF(IF+1,JF) + UF(IF,JF+1))
   10   CONTINUE
   20 CONTINUE
      END
C
C.......................................................................
C
      SUBROUTINE TRANS(NPC, NPF, UC, FC, UF, FF)
C
C       TRANSFER OF F AND U FROM GRID NPF TO NPC
C
      INTEGER NC, NPC, NPF, I, IF, J, JF
      DOUBLE PRECISION FC(NPC,NPC), FF(NPF,NPF), UC(NPC,NPC),
     *                 UF(NPF,NPF)
C
      NC = NPC - 1
C
      DO 20 J=2,NC
        JF = J + J - 1
        UC(1,J) = UF(1,JF)
        DO 10 I=2,NC
          IF = I + I - 1
          UC(I,J) = UF(IF,JF)
          FC(I,J) = 4.0D0*FF(IF,JF)
   10   CONTINUE
        UC(NPC,J) = UF(NPF,JF)
   20 CONTINUE
      DO 30 I=1,NPC
        IF = I + I - 1
        UC(I,1) = UF(IF,1)
        UC(I,NPC) = UF(IF,NPF)
   30 CONTINUE
      END
