

C      THIS PROGRAM IS TO MEASURE A COMPUTER PERFORMANCE IN MFLOPS
C      BY USING A KERNEL WHICH APPEARS IN A LINEAR SOLVER OF PRESSURE 
C      POISSON INCLUDED IN AN INCOMPRESSIBLE NAVIER-STOKES SOLVER. 
C      A POINT-JACOBI METHOD IS EMPLOYED IN THIS SOLVER.
C      ------------------
C      CODED BY RYUTARO HIMENO, NISSAN RESEARCH CENTER, NISSAN MOTOR CO.
C      E-MAIL:himeno@nova.lab.nissan.co.jp
C      FAX: 0468-65-5699 (FROM OUTSIDE OF JAPAN:+81-468-65-5699)
C      TEL.: 0468-67-5151 (FROM OUTSIDE OF JAPAN:+81-468-67-5151)
C      DATE NOVEMBER 24, 1994
C      ------------------
C      FINITE-DIFFERENCE METHOD, CURVILINEAR COODINATE SYSTEM
C      VECTORIZABLE AND PARALLELIZABLE ON EACH GRID POINT
C      NO. OF GRID POINTS : IMAX X JMAX X KMAX INCLUDING BOUNDARIES
C      ------------------
C      A,B,C:COEFFICIENT MATRIX, WRK1: SOURCE TERM OF POISSON EQUATION
C      WRK2 : WORKING AREA, OMEGA : RELAXATION PARAMETER
C      BND:CONTROL VARIABLE FOR BOUNDARIES AND OBJECTS ( = 0 OR 1)
C      P: PRESSURE 
C       -----------------
      PARAMETER (MIMAX=49,MJMAX=49,MKMAX=49)  
      DIMENSION P(MIMAX,MJMAX,MKMAX)
      DIMENSION A(MIMAX,MJMAX,MKMAX,4),
     1               B(MIMAX,MJMAX,MKMAX,3),C(MIMAX,MJMAX,MKMAX,3)
      DIMENSION BND(MIMAX,MJMAX,MKMAX),WRK1(MIMAX,MJMAX,MKMAX)
      DIMENSION WRK2(MIMAX,MJMAX,MKMAX)
      DIMENSION TIME1(2)
      DATA OMEGA /0.8/
C
      IMAX=MIMAX-1
      JMAX=MJMAX-1
      KMAX=MKMAX-1
      GOSA=0.0
C     INITIALIZING MATRIXES
      DO 10 K=1,MKMAX
        DO 10 J=1,MJMAX
          DO 10 I=1,MIMAX
            A(I,J,K,1)=FLOAT(I)/FLOAT(IMAX)
            A(I,J,K,2)=FLOAT(J)/FLOAT(JMAX)
            A(I,J,K,3)=FLOAT(K)/FLOAT(KMAX)
            A(I,J,K,4)=1.0/(A(I,J,K,1)+A(I,J,K,2)+A(I,J,K,3))
            B(I,J,K,1)=FLOAT(I)/FLOAT(IMAX)
            B(I,J,K,2)=FLOAT(J)/FLOAT(JMAX)
            B(I,J,K,3)=FLOAT(K)/FLOAT(KMAX)
            C(I,J,K,1)=FLOAT(I)/FLOAT(IMAX)
            C(I,J,K,2)=FLOAT(J)/FLOAT(JMAX)
            C(I,J,K,3)=FLOAT(K)/FLOAT(KMAX)
            P(I,J,K)=FLOAT(J*I)/FLOAT(JMAX*IMAX)
                        WRK1(I,J,K)=FLOAT(K*J*I)/FLOAT(KMAX*JMAX*IMAX)
            BND(I,J,K)=1.0
 10   CONTINUE
      WRITE(*,6003)MIMAX,MJMAX,MKMAX
 6003 FORMAT('  MIMAX=',I3,4X,'MJMAX=',I3,4X,'MKMAX=',I3)
      WRITE(*,*) ' IMAX=',IMAX,' JMAX=',JMAX,' KMAX=',KMAX
C
C     START  MEASURING (THREE TIMES)
        GOSA=0.0

      DO 900 KK=1,3
        EE0=DTIME(TIME0)
CF      CALL TIME(IEE0)  
CF      CALL CLOCKM(IE0)
C
        DO 100 K=2,KMAX-1
        DO 100 J=2,JMAX-1
        DO 100 I=2,IMAX-1
          S0=A(I,J,K,1)*P(I+1,J,K)+A(I,J,K,2)*P(I,J+1,K)
     1      +A(I,J,K,3)*P(I,J,K+1)
     3      +B(I,J,K,1)*(P(I+1,J+1,K)-P(I+1,J-1,K)
     *                  -P(I-1,J+1,K)+P(I-1,J-1,K))
     4      +B(I,J,K,2)*(P(I,J+1,K+1)-P(I,J-1,K+1)
     *                  -P(I,J+1,K-1)+P(I,J-1,K-1))
     5      +B(I,J,K,3)*(P(I+1,J,K+1)-P(I-1,J,K+1)
     *                  -P(I+1,J,K-1)+P(I-1,J,K-1))
     6      +C(I,J,K,1)*P(I-1,J,K)+C(I,J,K,2)*P(I,J-1,K)
     *      +C(I,J,K,3)*P(I,J,K-1)+WRK1(I,J,K)
          SS=(S0*A(I,J,K,4)-P(I,J,K))*BND(I,J,K)
          GOSA=GOSA+SS*SS
          WRK2(I,J,K)=P(I,J,K)+OMEGA *SS
  100   CONTINUE
        DO 200 K=2,KMAX-1
        DO 200 J=2,JMAX-1
        DO 200 I=2,IMAX-1
          P(I,J,K)=WRK2(I,J,K)
  200   CONTINUE
C
C       END OF ONE ITERATION
C
C**** CASE SUN CONVEX COMPILER ****
C
        EE1=DTIME(TIME1)
C
C**** CASE FUJITSU FRT COMPILER *****
C
CF      CALL TIME(IEE1)
CF      CALL CLOCKM(IE1)
CF      EE1 = 1.0E-3*FLOAT(IEE1-IEE0)
CF	FTIME1(1) = 1.0E-3*FLOAT(IE1-IE0)
C
        NFLOP=(KMAX-2)*(JMAX-2)*(IMAX-2)*34
        IF(EE1.NE.0.0) XMFLOPS=NFLOP/EE1*1.0E-6
        IF(TIME1(1).NE.0.0) XMFLOPS2=NFLOP/TIME1(1)*1.0E-6
        WRITE(*,6001) TIME1(1), EE1, GOSA
        WRITE(*,6004) XMFLOPS2, XMFLOPS
 6001 FORMAT(' CPUTIME=',F18.15,3X,'ELAPSED=',F18.15,3X,'GOSA=',E13.4)
 6004 FORMAT(' MFLOPS(CPUTIME)=',F10.5,7X,'MFLOPS(ELAPSTIME)=',F10.5)
        WRITE(*,*)' '
  900 CONTINUE
      STOP
      END
