      PROGRAM PDE1
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     MAIN PROGRAM FOR PDE1 BENCHMARK - SIMD/VECTOR VERSION
C
C     SUPPLIED BY   PALLAS GmbH
C
C     AUTHOR:  GENESIS BENCHMARK VERSION:
C		  M. LEMKE (1989, PALLAS GMBH)
C
C     FOR NUMERICAL VERIFICATION:
C	 THE AVERAGE CONVERGENCE RATE SHOULD BE SOMEWHERE BETWEEN
C	 0.7 AND 0.98
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     THE TIMER (DOUBLE PRECISION IN SECONDS) IS CALLED IN MAIN PROGRAM.
C     THE MACHINE DEPENDENT TIMING ROUTINE HAS TO BE INCLUDED.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      INCLUDE 'pde1.inc'
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      PARAMETER (NPMAX=2**(MMAX+6)+1)
      PARAMETER (NPMIN=2**(MMAX-3)+1)
      PARAMETER (IDIM=NPMAX*NPMIN**2)
      PARAMETER(NW1=11)

      DOUBLE PRECISION U(IDIM),F(IDIM)
      COMMON U,F
      ITER = MAXIT

      OPEN(NW1,FILE='result')
      CALL HEADER(NW1)
      WRITE(NW1,91000)
91000 FORMAT(' ***** 3D RED BLACK RELAXATION BENCHMARK ************* '/
     *    )
C
C---> LOOP OVER THE FINEST GRIDS IN ORDER TO GET RESULTS FOR DIFFERENT
C     PROBLEM SIZES
      DO 1 N=3,MMAX
        NP=2**N+1
        CALL TESTA(U,F,NP,NP,NP,ITER,NW1)
        CALL TESTB(U,F,NP,NP,NP,ITER,NW1)
        CALL TESTC(U,F,NP,NP,NP,ITER,NW1)
        ITER = MAX(ITER/8,5)
C
1     CONTINUE
      NX = NP-1
      NY = NP-1
      NZ = NP-1
      DO 2 N=1,3
         NX = NX * 4
         NY = NY / 2
         NZ = NZ / 2
         IF (NY.GE.3.AND.NZ.GE.3) THEN
            CALL TESTA(U,F,NX+1,NY+1,NZ+1,ITER,NW1)
            CALL TESTB(U,F,NX+1,NY+1,NZ+1,ITER,NW1)
            CALL TESTC(U,F,NX+1,NY+1,NZ+1,ITER,NW1)
         ENDIF
2     CONTINUE
C
      CLOSE(NW1)

      STOP
      END
C
C     ONE TESTCASE  STANDARD IMPLEMENTATION
C
      SUBROUTINE TESTA(U,F,NX,NY,NZ,ITER,NW1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NX,NY,NZ), F(NX,NY,NZ)
C
        H=1./(NX-1)
        CALL INIT(U,F,NX,NY,NZ)
        RES1=RES(U,F,NX,NY,NZ,H)
C----------------> TIMER IN SECONDS
        CALL TIMER(SEC1)
        DO 2 NREL=1,ITER
          CALL RELAXA(U,F,NX,NY,NZ,H)
2       CONTINUE
        CALL TIMER(SEC2)
C----------------> TIMER IN SECONDS (IN SEC2)
        RES2=RES(U,F,NX,NY,NZ,H)
        CONVR=(RES2/RES1)**(0.1)
        WRITE(NW1,91100)
        WRITE(NW1,9100)NX,NY,NZ,NX*NY*NZ, ITER, CONVR
C
        NOP = (NX-2)*(NY-2)*(NZ-2)*8
C
        RMFLOP = (NOP/1000000.)/(SEC2-SEC1)*ITER
        WRITE(NW1,900) (SEC2-SEC1)/ITER,NOP, RMFLOP
        WRITE(NW1,91100)
91100 FORMAT(' +++++++++   STANDARD IMPLEMENTATION  ++++++++++++++++ '/
     *     )
 9100 FORMAT (/' NUMBER OF POINTS PER DIRECTION = ', 3I12 /
     *        ' GRIDSIZE OF CUBIC  GRID =        ', I12/
     *        ' AVERAGE CONV. RATE (', I6, 'ITER.)  = ', D12.4/)
900    FORMAT(/' TIME FOR 1 RELAX (SEC.):', E15.8/
     *           ' NUMBER OF OPERATIONS:   ',I15/
     *           ' PERFORMANCE IN MFLOP/S: ',F15.5/)
      END
      DOUBLE PRECISION FUNCTION RES(U,F,NX,NY,NZ,H)
C
C	COMPUTES THE MAXIMUM NORM OF THE DEFECT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION U(NX,NY,NZ), F(NX,NY,NZ)
C
      RES = 0.
      FAC=1/(H*H)
      DO 30 K=2,NZ-1
        DO 20 J=2,NY-1
          DO 10 I=2,NX-1
            D = F(I,J,K)+FAC*(-6.0*U(I,J,K)+U(I-1,J,K)+U(I+1,J,K)
     *                                     +U(I,J-1,K)+U(I,J+1,K)
     *                                     +U(I,J,K-1)+U(I,J,K+1))
               RES = MAX(ABS(D),RES)
   10     CONTINUE
   20   CONTINUE
   30 CONTINUE
      END
C
C.......................................................................
C
      SUBROUTINE INIT(U,F,NX,NY,NZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION U(NX,NY,NZ),F(NX,NY,NZ)
      DO 10 K = 1, NZ-1
        DO 20 J = 1, NY-1
          DO 30 I = 1, NX-1
            U(I,J,K) = 0.
30        CONTINUE
20      CONTINUE
10    CONTINUE
C
      DO 11 K = 2, NZ-1
        DO 21 J = 2, NY-1
          DO 31 I = 2, NX-1
            U(I,J,K) = 1000.
            F(I,J,K) = 10.
31        CONTINUE
21      CONTINUE
11    CONTINUE
C
      END
      SUBROUTINE TESTB(U,F,NX,NY,NZ,ITER,NW1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NX,NY,NZ), F(NX,NY,NZ)
C
        H=1./(NX-1)
        CALL INIT(U,F,NX,NY,NZ)
        RES1=RES(U,F,NX,NY,NZ,H)
C----------------> TIMER IN SECONDS
        CALL TIMER(SEC1)
        DO 2 NREL=1,ITER
          CALL RELAXB(U,F,NX,NY,NZ,H)
2       CONTINUE
        CALL TIMER(SEC2)
C----------------> TIMER IN SECONDS (IN SEC2)
        RES2=RES(U,F,NX,NY,NZ,H)
        CONVR=(RES2/RES1)**(0.1)
        WRITE(NW1,91100)
        WRITE(NW1,9100)NX,NY,NZ,NX*NY*NZ, ITER, CONVR
C
        NOP = (NX-2)*(NY-2)*(NZ-2)*8
C
        RMFLOP = (NOP/1000000.)/(SEC2-SEC1)*ITER
        WRITE(NW1,900) (SEC2-SEC1)/ITER,NOP, RMFLOP
        WRITE(NW1,91100)
91100 FORMAT(' +++++++++   UNROLLING OF DEPTH 2 IN EACH DIMENSION  + '/
     *     )
 9100 FORMAT (/' NUMBER OF POINTS PER DIRECTION = ', 3I12 /
     *        ' GRIDSIZE OF CUBIC  GRID =        ', I12/
     *        ' AVERAGE CONV. RATE (', I6, 'ITER.)  = ', D12.4/)
900     FORMAT(/' TIME FOR 1 RELAX (SEC.):', E15.8/
     *         ' NUMBER OF OPERATIONS:   ',I15/
     *         ' PERFORMANCE IN MFLOP/S: ',F15.5/)
      END
C
C     ONE TESTCASE  SINGLE LOOP IN TWO DIMENSION
C
      SUBROUTINE TESTC(U,F,NX,NY,NZ,ITER,NW1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NX,NY,NZ), F(NX,NY,NZ)
C
        H=1./(NX-1)
        CALL INIT(U,F,NX,NY,NZ)
        RES1=RES(U,F,NX,NY,NZ,H)
C----------------> TIMER IN SECONDS
        CALL TIMER(SEC1)
        DO 2 NREL=1,ITER
          CALL RELAXC(U,F,NX,NY,NZ,H)
2       CONTINUE
        CALL TIMER(SEC2)
C----------------> TIMER IN SECONDS (IN SEC2)
        RES2=RES(U,F,NX,NY,NZ,H)
        CONVR=(RES2/RES1)**(0.1)
        WRITE(NW1,91100)
        WRITE(NW1,9100)NX,NY,NZ,NX*NY*NZ, ITER, CONVR
C
        NOP = (NX-2)*(NY-2)*(NZ-2)*8
C
        RMFLOP = (NOP/1000000.)/(SEC2-SEC1)*ITER
        WRITE(NW1,900) (SEC2-SEC1)/ITER,NOP, RMFLOP
        WRITE(NW1,91100)
91100 FORMAT(' +++++++++   SINGLE LOOP IN 2 DIMENSIONS  +++++++++++ '/
     *     )
 9100 FORMAT (/' NUMBER OF POINTS PER DIRECTION = ', 3I12 /
     *        ' GRIDSIZE OF CUBIC  GRID =        ', I12/
     *        ' AVERAGE CONV. RATE (', I6, 'ITER.)  = ', D12.4/)
900     FORMAT(/' TIME FOR 1 RELAX (SEC.):', E15.8/
     *        ' NUMBER OF OPERATIONS:   ',I15/
     *          ' PERFORMANCE IN MFLOP/S: ',F15.5/)
      END
C
C
C
      SUBROUTINE RELAXA(U,F,NX,NY,NZ,H)
C
C     3D RED BLACK RELAXATION    STANDARD IMLEMENTATION
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION U(NX,NY,NZ),F(NX,NY,NZ)
      H2 = H*H
      FACTOR = 1.0D0/6.0D0
C
C     RELAXATION OF THE RED POINTS
C
      IA = 0
      DO 10 K = 2, NZ-1
        DO 20 J = 2, NY-1
          DO 30 I = 2 + IA, NX-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
30        CONTINUE
        IA = MOD (IA+1,2)
20      CONTINUE
10    CONTINUE
C
C     RELAXATION OF THE BLACK POINTS
C
      IA = 1
      DO 11 K = 2, NZ-1
        DO 21 J = 2, NY-1
          DO 31 I = 2 + IA, NX-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
31        CONTINUE
        IA = MOD (IA+1,2)
21      CONTINUE
11    CONTINUE
      END
C
C.......................................................................
C
      SUBROUTINE RELAXB(U,F,NX,NY,NZ,H)
C
C     3D RED BLACK RELAXATION    STANDARD IMLEMENTATION
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION U(NX,NY,NZ),F(NX,NY,NZ)
      INTEGER COLOR
      H2 = H*H
      FACTOR = 1.0D0/6.0D0
      DO 1 COLOR = 0,1
      NZ1 = NZ
      IF (MOD(NZ,2).EQ.1) THEN
        K = NZ-1
        DO 120 J = 2, NY-1
          IINC = MOD( K+J+2+COLOR,2)
          DO 130 I = 2 + IINC, NX-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
130       CONTINUE
120     CONTINUE
        NZ1 = NZ1 - 1
      ENDIF
      NY1 = NY
      IF (MOD(NY,2).EQ.1) THEN
        J = NY-1
        DO 210 K = 2, NZ1-1
          IINC = MOD( K+J+2+COLOR,2)
          DO 230 I = 2 + IINC, NX-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
230       CONTINUE
210     CONTINUE
        NY1 = NY1 - 1
      ENDIF
      NX1 = NX
      IF (MOD(NX,2).EQ.1) THEN
        I = NX-1
        DO 310 K = 2, NZ1-1
          JINC = MOD( K+2+I+COLOR,2)
          DO 320 J = 2 + JINC, NY1-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
320       CONTINUE
310     CONTINUE
        NX1 = NX1 - 1
      ENDIF
      IF (COLOR.EQ.0) THEN
      DO 10 K = 2, NZ1-1,2
        DO 20 J = 2, NY1-1,2
          DO 30 I = 2 , NX1-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
            U(I+1,J+1,K) = FACTOR*(H2*F(I+1,J+1,K) +
     *                    (U(I,J+1,K) + U(I+2,J+1,K) +
     *                     U(I+1,J,K) + U(I+1,J+2,K) +
     *                     U(I+1,J+1,K-1) + U(I+1,J+1,K+1)))
            U(I,J+1,K+1) = FACTOR*(H2*F(I,J+1,K+1) +
     *                    (U(I-1,J+1,K+1) + U(I+1,J+1,K+1) +
     *                     U(I,J,K+1) + U(I,J+2,K+1) +
     *                     U(I,J+1,K) + U(I,J+1,K+2)))
            U(I+1,J,K+1) = FACTOR*(H2*F(I+1,J,K+1) +
     *                    (U(I,J,K+1) + U(I+2,J,K+1) +
     *                     U(I+1,J-1,K+1) + U(I+1,J+1,K+1) +
     *                     U(I+1,J,K) + U(I+1,J,K+2)))
30        CONTINUE
20      CONTINUE
10    CONTINUE
      ELSE
      DO 1010 K = 2, NZ1-1,2
        DO 1020 J = 2, NY1-1,2
          DO 1030 I = 2 , NX1-1, 2
            U(I+1,J,K) = FACTOR*(H2*F(I+1,J,K) +
     *                    (U(I,J,K) + U(I+2,J,K) +
     *                     U(I+1,J-1,K) + U(I+1,J+1,K) +
     *                     U(I+1,J,K-1) + U(I+1,J,K+1)))
            U(I,J+1,K) = FACTOR*(H2*F(I,J+1,K) +
     *                    (U(I-1,J+1,K) + U(I+1,J+1,K) +
     *                     U(I,J,K) + U(I,J+2,K) +
     *                     U(I,J+1,K-1) + U(I,J+1,K+1)))
            U(I,J,K+1) = FACTOR*(H2*F(I,J,K+1) +
     *                    (U(I-1,J,K+1) + U(I+1,J,K+1) +
     *                     U(I,J-1,K+1) + U(I,J+1,K+1) +
     *                     U(I,J,K) + U(I,J,K+2)))
            U(I+1,J+1,K+1) = FACTOR*(H2*F(I+1,J+1,K+1) +
     *                    (U(I,J+1,K+1) + U(I+2,J+1,K+1) +
     *                     U(I+1,J,K+1) + U(I+1,J+2,K+1) +
     *                     U(I+1,J+1,K) + U(I+1,J+1,K+2)))
1030      CONTINUE
1020    CONTINUE
1010  CONTINUE
      ENDIF
 1    CONTINUE
      END
C
C.......................................................................
C
      SUBROUTINE RELAXC(U,F,NX,NY,NZ,H)
C
C     3D RED BLACK RELAXATION    SINGLE LOOP IN 2 DIMENSIONS
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION U(NX*NY,NZ),F(NX*NY,NZ)
      INTEGER COLOR
      INTEGER COLOR1
      IK(I,J) = I + (J-1) * NX
      H2 = H*H
      FACTOR = 1.0D0/6.0D0
      IF (MOD(NX,2).NE.1) THEN
         WRITE(6,*) 'WRONG NUMBER OF GRIDPOINTS FOR TESTC'
      ENDIF
      DO 1 COLOR=0,1
      DO 10 K = 2, NZ-1
        COLOR1 = MOD(K+COLOR,2)
        I = 1
        DO 21 J = 2 + 1-COLOR1, NY-1
           F(I+(J-1)*NX,K) = U( I+(J-1)*NX,K)
 21     CONTINUE
        I = NX
        DO 22 J = 2 + 1-COLOR1, NY-1
          F(I+(J-1)*NX,K) = U( I+(J-1)*NX,K)
 22     CONTINUE
          DO 30 I = 2 + COLOR1+NX , NX-1+(NY-2)*NX, 2
            U(I,K) = FACTOR*(H2*F(I,K) +
     *                    (U(I-1,K) + U(I+1,K) +
     *                     U(I-NX,K) + U(I+NX,K) +
     *                     U(I,K-1) + U(I,K+1)))
30        CONTINUE
20      CONTINUE
        I = 1
        DO 23 J = 2 + 1-COLOR1, NY-1
           U(I+(J-1)*NX,K) = F( I+(J-1)*NX,K)
 23     CONTINUE
        I = NX
        DO 24 J = 2 + 1-COLOR1, NY-1
           U(I+(J-1)*NX,K) = F(I+(J-1)*NX,K)
 24     CONTINUE
10    CONTINUE
1     CONTINUE 	
      END
