      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              SUBSET FORTRAN 90 REWRITE:
C                 D.B. CARPENTER (1993, UNIV. OF SOUTHAMPTON)
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      IMPLICIT NONE
      DOUBLE PRECISION F, G, GETDIF, GETDEF, DIF, DEF
      INTEGER JDIM, IER, M, N, NP, NY1, NY2, NCYCLE, NFMG, IGAM
      INTEGER NOP, LEV
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
      INTEGER NW1
      PARAMETER(NW1=11)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DOUBLE PRECISION DEF1, DEF2, CONVR
      DOUBLE PRECISION SEC2, RMFLOP
      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 M = 3, MMAX
C
        NY1 = 2
        NY2 = 1
        NFMG = 0
        IGAM = 1
C
        CALL MG00D (M, NY1, NY2, NFMG, NCYCLE, IGAM, F, G, JDIM, IER)
C
        IF (IER .NE. 0) THEN
          WRITE (NW1,9000) IER, JDIM
          STOP
        ENDIF
        N = 2**M
        NP = N + 1
        DIF = GETDIF (M, G)
        DEF = GETDEF (M)
        CONVR=(DEF/DEF1)**(0.1)
        WRITE (NW1,9100) NP, NP**2, DIF, DEF, CONVR
C
        NOP = 1
        DO 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
        ENDDO
        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)
      ENDDO
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 F(X,Y)
      IMPLICIT NONE
      DOUBLE PRECISION X, Y
        F = 10.0D0 * SIN (3.0D0 * X + Y)
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION G (X, Y)
      IMPLICIT NONE
      DOUBLE PRECISION X, Y
        G = SIN (3.0D0 * X + Y)
      END
C
C***********************************************************************
C
      SUBROUTINE MG00D (M, NY1, NY2, NFMG, NCYCLE, IGAM, F, G,
     1                  JDIM, IER)
      IMPLICIT NONE
      INTEGER M, NY1, NY2, NFMG, NCYCLE, IGAM, JDIM, IER
      DOUBLE PRECISION F, G
      EXTERNAL F, G
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        JDIM    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 LIN, IT, LEV
        DOUBLE PRECISION  SEC1, SEC2, DEF1, DEF2
        COMMON /TIME/ SEC2
        COMMON /DEFECT/ DEF1,DEF2
        DOUBLE PRECISION GETDEF
        INCLUDE 'stack.inc'

        IER = 0
C
C       CHECK OF ERROR CONDITIONS (WRONG PARAMETERS)
C
        IF ((M .LE. 0) .OR. (M .GE. 110) .OR. (NY1 .LE. 0) .OR.
     1      (NY2 .LE. 0) .OR. (NCYCLE .LE. 0) .OR. (IGAM .LE. 0)) THEN
          IER = 2
          RETURN
        ENDIF
C
C       INITIALISE STACK...
C
        CALL INISTK (M, IER, JDIM)
        IF (IER .NE. 0) RETURN 
C
C       LIN = LEVEL FROM WHICH INITIAL GUESS PROPAGATED DOWN
C
        IF (NFMG .EQ. 0) THEN
          LIN = 1
        ELSE
          LIN = M
        ENDIF
C
C       SET UP ALL GRID VALUES NEEDED
C
        CALL INIT (M, LIN, F, G)
C
C       (FULL) MULTIGRID PROCEDURE...
C
C       INITIAL GUESS FROM COARSEST GRID
C
        DO LEV = LIN, 2, -1
          CALL MGI (LEV, M, NY1, NY2, IGAM)
C----> BEGIN STACK ACCESS
          GOTO (101, 102, 103, 104, 105, 106, 107) M - LEV + 1

 101      CONTINUE
            CALL FMGDN (NP2, U2, F2, NP1, U1)
          GOTO 150
 102      CONTINUE
            CALL FMGDN (NP3, U3, F3, NP2, U2)
          GOTO 150
 103      CONTINUE
            CALL FMGDN (NP4, U4, F4, NP3, U3)
          GOTO 150
 104      CONTINUE
            CALL FMGDN (NP5, U5, F5, NP4, U4)
          GOTO 150
 105      CONTINUE
            CALL FMGDN (NP6, U6, F6, NP5, U5)
          GOTO 150
 106      CONTINUE
            CALL FMGDN (NP7, U7, F7, NP6, U6)
          GOTO 150
 107      CONTINUE
            CALL FMGDN (NP8, U8, F8, NP7, U7)
          GOTO 150

 150      CONTINUE
C----> END STACK ACCESS
        ENDDO
C
C       MULTIGRID CYCLES
C
        DEF1 = GETDEF (M)
C----> TIMING
        CALL TIMER(SEC1)
        DO IT = 1, NCYCLE
          CALL MGI (1, M, NY1, NY2, IGAM)
        ENDDO
        CALL TIMER(SEC2)
        SEC2 = SEC2 - SEC1
C----> TIMING

      END
C
C.......................................................................
C
      SUBROUTINE INISTK (M, IER, JDIM)
      IMPLICIT NONE
      INTEGER M, IER, JDIM
C
C     INITIALISE STACK
C
C        INTEGER NP, L
C        INCLUDE 'stack.inc'

C----> BEGIN STACK ACCESS
        IF (M .GT. 8) THEN
          IER  = 1
          JDIM = 8
        ENDIF
C----> END STACK ACCESS
      END
C
C.......................................................................
C
      SUBROUTINE INIT (M, LIN, F, G)
      IMPLICIT NONE
      INTEGER M, LIN
      DOUBLE PRECISION F, G
      EXTERNAL F, G
C
C     COMPUTES INITIAL VALUES ON ALL GRIDS.
C     LIN = LEVEL FROM WHICH INITIAL GUESS PROPAGATED DOWN
C
        INTEGER L
        DOUBLE PRECISION H
        INCLUDE 'stack.inc'

        H = 1.0D0 / DBLE (FLOAT (2**M))
C
C       INITIALISATION ON FINEST LEVEL
C
C----> BEGIN STACK ACCESS
        GOTO (101, 102, 103, 104, 105, 106, 107, 108) M

 101    CONTINUE
          CALL INIT1 (NP1, H, F, G, U1, F1)
        GOTO 150
 102    CONTINUE
          CALL INIT1 (NP2, H, F, G, U2, F2)
        GOTO 150
 103    CONTINUE
          CALL INIT1 (NP3, H, F, G, U3, F3)
        GOTO 150
 104    CONTINUE
          CALL INIT1 (NP4, H, F, G, U4, F4)
        GOTO 150
 105    CONTINUE
          CALL INIT1 (NP5, H, F, G, U5, F5)
        GOTO 150
 106    CONTINUE
          CALL INIT1 (NP6, H, F, G, U6, F6)
        GOTO 150
 107    CONTINUE
          CALL INIT1 (NP7, H, F, G, U7, F7)
        GOTO 150
 108    CONTINUE
          CALL INIT1 (NP8, H, F, G, U8, F8)
        GOTO 150

 150    CONTINUE
C----> END STACK ACCESS
C
C       TRANSFER OF F AND U TO COARSER GRIDS
C
        DO L = 2, LIN
C----> BEGIN STACK ACCESS
          GOTO (201, 202, 203, 204, 205, 206, 207) M - L + 1

 201      CONTINUE
            CALL TRANS (NP1, NP2, U1, F1, U2, F2)
          GOTO 250
 202      CONTINUE
            CALL TRANS (NP2, NP3, U2, F2, U3, F3)
          GOTO 250
 203      CONTINUE
            CALL TRANS (NP3, NP4, U3, F3, U4, F4)
          GOTO 250
 204      CONTINUE
            CALL TRANS (NP4, NP5, U4, F4, U5, F5)
          GOTO 250
 205      CONTINUE
            CALL TRANS (NP5, NP6, U5, F5, U6, F6)
          GOTO 250
 206      CONTINUE
            CALL TRANS (NP6, NP7, U6, F6, U7, F7)
          GOTO 250
 207      CONTINUE
            CALL TRANS (NP7, NP8, U7, F7, U8, F8)
          GOTO 250

 250      CONTINUE
C----> END STACK ACCESS
        ENDDO
C
C       PUT ZERO TO BOUNDARY VALUES FOR MULTIGRID CORRECTIONS
C
        DO L = LIN + 1, M
C----> BEGIN STACK ACCESS
          GOTO (301, 302, 303, 304, 305, 306, 307, 308) M - L + 1

 301      CONTINUE
            CALL PUTZB (NP1, U1)
          GOTO 350
 302      CONTINUE
            CALL PUTZB (NP2, U2)
          GOTO 350
 303      CONTINUE
            CALL PUTZB (NP3, U3)
          GOTO 350
 304      CONTINUE
            CALL PUTZB (NP4, U4)
          GOTO 350
 305      CONTINUE
            CALL PUTZB (NP5, U5)
          GOTO 350
 306      CONTINUE
            CALL PUTZB (NP6, U6)
          GOTO 350
 307      CONTINUE
            CALL PUTZB (NP7, U7)
          GOTO 350
 308      CONTINUE
            CALL PUTZB (NP8, U8)
          GOTO 350

 350      CONTINUE
C----> END STACK ACCESS
        ENDDO
      END
C
C.......................................................................
C
      SUBROUTINE MGI(LEV, M, NY1, NY2, IGAM)
      IMPLICIT NONE
      INTEGER LEV, M, NY1, NY2, IGAM
C
C     ONE MULTIGRID ITERATION STEP
C
C        INTEGER IDC, IDF, IRC, IRF, K, ICGAM (10), IGAML (10)
        INTEGER K, ICGAM (10), IGAML (10)
        LOGICAL IZ
        INCLUDE 'stack.inc'

        DO K = LEV, M
          IF ((K .EQ. LEV) .OR. (K .EQ. M)) THEN
            IGAML (K) = 1
          ELSE
            IGAML (K) = IGAM
          ENDIF
          ICGAM(K) = 0
        ENDDO
        K = LEV

        DO WHILE (ICGAM (K) .LT. IGAML (K))
C
C         ``UP''
C
          DO WHILE (K .LT. M)
            IZ = K .GT. LEV .AND. ICGAM(K) .EQ. 0

C----> BEGIN STACK ACCESS
            GOTO (101, 102, 103, 104, 105, 106, 107) M - K

 101        CONTINUE
              CALL MGIUP (NP1, F1, NP2, U2, F2, NY1, IZ)
            GOTO 150
 102        CONTINUE
              CALL MGIUP (NP2, F2, NP3, U3, F3, NY1, IZ)
            GOTO 150
 103        CONTINUE
              CALL MGIUP (NP3, F3, NP4, U4, F4, NY1, IZ)
            GOTO 150
 104        CONTINUE
              CALL MGIUP (NP4, F4, NP5, U5, F5, NY1, IZ)
            GOTO 150
 105        CONTINUE
              CALL MGIUP (NP5, F5, NP6, U6, F6, NY1, IZ)
            GOTO 150
 106        CONTINUE
              CALL MGIUP (NP6, F6, NP7, U7, F7, NY1, IZ)
            GOTO 150
 107        CONTINUE
              CALL MGIUP (NP7, F7, NP8, U8, F8, NY1, IZ)
            GOTO 150

 150        CONTINUE
C----> END STACK ACCESS

            ICGAM(K) = ICGAM(K) + 1
            K = K + 1
          ENDDO
C
C         EXACT SOLUTION ON COARSEST GRID (K = M)
C
C----> BEGIN STACK ACCESS
          CALL RELAX(1, (M .GT. LEV .AND. ICGAM(M) .EQ. 0), NP1, U1, F1)
C----> END STACK ACCESS
          ICGAM (M) = ICGAM (M) + 1
C
C         ``DOWN''
C
          DO WHILE ((K .GT. LEV) .AND. (ICGAM (K) .EQ. IGAML (K)))
            ICGAM (K) = 0
            K = K - 1

C----> BEGIN STACK ACCESS
            GOTO (201, 202, 203, 204, 205, 206, 207) M - K

 201        CONTINUE
              CALL MGIDN (NP2, U2, F2, NP1, U1, NY2)
            GOTO 250
 202        CONTINUE
              CALL MGIDN (NP3, U3, F3, NP2, U2, NY2)
            GOTO 250
 203        CONTINUE
              CALL MGIDN (NP4, U4, F4, NP3, U3, NY2)
            GOTO 250
 204        CONTINUE
              CALL MGIDN (NP5, U5, F5, NP4, U4, NY2)
            GOTO 250
 205        CONTINUE
              CALL MGIDN (NP6, U6, F6, NP5, U5, NY2)
            GOTO 250
 206        CONTINUE
              CALL MGIDN (NP7, U7, F7, NP6, U6, NY2)
            GOTO 250
 207        CONTINUE
              CALL MGIDN (NP8, U8, F8, NP7, U7, NY2)
            GOTO 250

 250        CONTINUE
C----> END STACK ACCESS

          ENDDO
        ENDDO
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION GETDEF (M)
      IMPLICIT NONE
      INTEGER M
C
C     EXTRACT MAXIMUM NORM OF THE DEFECT FROM STACK
C
        DOUBLE PRECISION DEFMX
        INCLUDE 'stack.inc'

C----> BEGIN STACK ACCESS
        GOTO (101, 102, 103, 104, 105, 106, 107, 108) M

 101    CONTINUE
          GETDEF = DEFMX (U1, F1, NP1)
        GOTO 150
 102    CONTINUE
          GETDEF = DEFMX (U2, F2, NP2)
        GOTO 150
 103    CONTINUE
          GETDEF = DEFMX (U3, F3, NP3)
        GOTO 150
 104    CONTINUE
          GETDEF = DEFMX (U4, F4, NP4)
        GOTO 150
 105    CONTINUE
          GETDEF = DEFMX (U5, F5, NP5)
        GOTO 150
 106    CONTINUE
          GETDEF = DEFMX (U6, F6, NP6)
        GOTO 150
 107    CONTINUE
          GETDEF = DEFMX (U7, F7, NP7)
        GOTO 150
 108    CONTINUE
          GETDEF = DEFMX (U8, F8, NP8)
        GOTO 150

 150    CONTINUE
C----> END STACK ACCESS
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION GETDIF (M, SOL)
      IMPLICIT NONE
      INTEGER M
      DOUBLE PRECISION SOL
      EXTERNAL SOL
C
C     EXTRACT MAXIMUM NORM OF THE DIFFERENCE FROM STACK
C
        DOUBLE PRECISION DIFMX
        INCLUDE 'stack.inc'

C----> BEGIN STACK ACCESS
        GOTO (101, 102, 103, 104, 105, 106, 107, 108) M

 101    CONTINUE
          GETDIF = DIFMX (U1, NP1, SOL)
        GOTO 150
 102    CONTINUE
          GETDIF = DIFMX (U2, NP2, SOL)
        GOTO 150
 103    CONTINUE
          GETDIF = DIFMX (U3, NP3, SOL)
        GOTO 150
 104    CONTINUE
          GETDIF = DIFMX (U4, NP4, SOL)
        GOTO 150
 105    CONTINUE
          GETDIF = DIFMX (U5, NP5, SOL)
        GOTO 150
 106    CONTINUE
          GETDIF = DIFMX (U6, NP6, SOL)
        GOTO 150
 107    CONTINUE
          GETDIF = DIFMX (U7, NP7, SOL)
        GOTO 150
 108    CONTINUE
          GETDIF = DIFMX (U8, NP8, SOL)
        GOTO 150

 150    CONTINUE
C----> END STACK ACCESS
      END
C
C.......................................................................
C
      SUBROUTINE FMGDN (NPF, UF, FRF, NPC, UC)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION UF (NPF, NPF), FRF (NPF, NPF), UC (NPC, NPC)
CHPF$ ALIGN UF (I, J) WITH * FRF (I, J)
CHPF$ DISTRIBUTE FRF * (BLOCK, BLOCK)
CHPF$ DISTRIBUTE UC * (BLOCK, BLOCK)
C
C     ONE ``DOWN'' STEP IN FULL MULTIGRID
C
        CALL INT4  (NPC, NPF, UC, UF, FRF)
        CALL PUTZB (NPC, UC)
      END
C
C.......................................................................
C
      SUBROUTINE MGIUP (NPC, FRC, NPF, UF, FRF, NY1, IZ)
      IMPLICIT NONE
      INTEGER NY1, NPC, NPF
      LOGICAL IZ
      DOUBLE PRECISION UF (NPF, NPF), FRF (NPF, NPF), FRC (NPC, NPC)
CHPF$ ALIGN UF (I, J) WITH * FRF (I, J)
CHPF$ DISTRIBUTE FRF * (BLOCK, BLOCK)
CHPF$ DISTRIBUTE FRC * (BLOCK, BLOCK)
C
C     ONE ``UP'' STEP IN MULTIGRID ITERATION.
C
C       RELAXATIONS BEFORE COARSE-GRID CORRECTION...
C
        CALL RELAX (NY1, IZ, NPF, UF, FRF)
C
C       RESIDUAL TRANSFER TO NEXT COARSER GRID...
C
        CALL RESTR (NPC, NPF, FRC, UF, FRF)

      END
C
C.......................................................................
C
      SUBROUTINE MGIDN (NPF, UF, FRF, NPC, UC, NY2)
      IMPLICIT NONE
      INTEGER NY2, NPC, NPF
      DOUBLE PRECISION UF (NPF, NPF), FRF (NPF, NPF), UC (NPC, NPC)
CHPF$ ALIGN UF (I, J) WITH * FRF (I, J)
CHPF$ DISTRIBUTE FRF * (BLOCK, BLOCK)
CHPF$ DISTRIBUTE UC * (BLOCK, BLOCK)
C
C     ONE ``DOWN'' STEP IN MULTIGRID ITERATION.
C
C       LINEAR INTERPOLATION TO NEXT FINER GRID
C
        CALL INT2A (NPC, NPF, UC, UF)
C
C       RELAXATION AFTER CGC
C
        CALL RELAX (NY2, .FALSE., NPF, UF, FRF)
      END

