      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, IDC, IDF, IDR, 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
          IDF = ID (LEV - 1)
          IDC = ID (LEV)
          IDR = IDF + INITF
          CALL FMGDN (NPK (LEV - 1), W (IDF), W (IDR),
     1                NPK (LEV), W (IDC))
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
        INTEGER NP, L
        INCLUDE 'stack.inc'

C----> BEGIN STACK ACCESS
        NP = 2**M + 1
        ID  (1) = 1
        NPK (1) = NP
        ID  (2) = 1 + NPK (1)**2
        DO L = 2, M
          NPK (L) = (NPK (L - 1) + 1) / 2
          ID (L + 1) = ID (L) + NPK (L)**2
        ENDDO
C
C       CHECK OF ERROR CONDITIONS (INSUFFICIENT MEMORY)
C
        INITF = ID (M + 1)
        IF (2 * INITF .GT. IDIM) THEN
          IER  = 1
          JDIM = 2 * INITF
        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 IDC, IDF, IRC, IRF, 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
        CALL INIT1 (NPK (1), H, F, G, W (1), W (INITF + 1))
C----> END STACK ACCESS
C
C       TRANSFER OF F AND U TO COARSER GRIDS
C
        DO L = 2, LIN
C----> BEGIN STACK ACCESS
          IDC = ID (L)
          IDF = ID (L - 1)
          IRC = INITF + IDC
          IRF = INITF + IDF
          CALL TRANS (NPK (L), NPK (L - 1), W (IDC), W (IRC),
     1                W (IDF), W (IRF))
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
          IDC = ID (L)
          CALL PUTZB (NPK (L), W (IDC))
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
        INTEGER IDC, IDF, IRC, IRF, 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
            IDF = ID (K)
            IDC = ID (K + 1)
            IRF = INITF + IDF
            IRC = INITF + IDC
            CALL MGIUP (NPK (K + 1), W (IRC),
     1                  NPK (K), W (IDF), W (IRF), NY1, IZ)
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
          IDF = ID (M)
          IRF = INITF + IDF
          CALL RELAX(1, (M .GT. LEV .AND. ICGAM(M) .EQ. 0), NPK (M),
     1               W (IDF), W (IRF))
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
            IDF = ID (K)
            IDC = ID (K + 1)
            IRF = INITF + IDF
            IRC = INITF + IDC
            CALL MGIDN (NPK (K), W (IDF), W (IRF),
     1                  NPK (K + 1), W (IDC), NY2)
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
        GETDEF = DEFMX (W, W (INITF + 1), NPK (1))
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
        GETDIF = DIFMX (W, NPK (1), SOL)
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)
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)
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)
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

