C*********************************************************************
C
C This benchmark test program is measuring a cpu performance
C of floating point operation by a Poisson equation solver.
CC
C If you have any question, please ask me via email.
C written by Ryutaro HIMENO, November 26, 2001.
C Version 3.0
C ----------------------------------------------
C Ryutaro Himeno, Dr. of Eng.
C Head of Computer Information Division,
C RIKEN (The Institute of Pysical and Chemical Research)
C Email : himeno@postman.riken.go.jp
C ---------------------------------------------------------------
C You can adjust the size of this benchmark code to fit your target
C computer. In that case, please chose following sets of
C (mimax,mjmax,mkmax):
C small : 65,33,33
C small : 129,65,65
C midium: 257,129,129
C large : 513,257,257
C ext.large: 1025,513,513
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 eq. which appears in an incompressible Navier-Stokes solver.
C A point-Jacobi method is employed in this solver as this method can
C be easyly vectrized and be parallelized.
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 -------------------
C -------------------
C "use portlib" statement on the next line is for Visual fortran
C to use UNIX libraries. Please remove it if your system is UNIX.
C -------------------
      IMPLICIT REAL*4(a-h,o-z)
C
      PARAMETER (mimax=513,mjmax=257,mkmax=257)
C      PARAMETER (mimax=257,mjmax=129,mkmax=129)
C      PARAMETER (mimax=129,mjmax=65,mkmax=65)
C      PARAMETER (mimax=65,mjmax=33,mkmax=33)
C
C     ttarget specifys the measuring period in sec
        PARAMETER (ttarget=60.0)
CC Arrey
      common /pres/ p(mimax,mjmax,mkmax)
      common /mtrx/ a(mimax,mjmax,mkmax,4),
     +     b(mimax,mjmax,mkmax,3),c(mimax,mjmax,mkmax,3)
      common /bound/ bnd(mimax,mjmax,mkmax)
      common /work/ wrk1(mimax,mjmax,mkmax),wrk2(mimax,mjmax,mkmax)
CC Other constants
      common /others/ imax,jmax,kmax,omega
C
      real*8 t1,t2
!HPF$ DISTRIBUTE (*,*,BLOCK):: p,bnd,wrk1,wrk2
!HPF$ DISTRIBUTE (*,*,BLOCK,*) :: a,b,c
C
      omega=0.8
      imax=mimax-1
      jmax=mjmax-1
      kmax=mkmax-1
CC Initializing matrixes
      call initmt
      write(*,*) ' mimax=',mimax,' mjmax=',mjmax,' mkmax=',mkmax
      write(*,*) ' imax=',imax,' jmax=',jmax,' kmax=',kmax
CC Start measuring
C
      nn=1000
        write(*,*) ' Start rehearsal measurement process.'
        write(*,*) ' Measure the performance in 1000 times.'
C
      gosa=0.0
      call etime(t1)
C
C Jacobi iteration
      call jacobi(nn,gosa)
C
      call etime(t2)
      flop=real(kmax-2)*real(jmax-2)*real(imax-2)*34.0*real(nn)
      xmflops2=flop/(t2-t1)*1.0e-6
        write(*,*) '  MFLOPS:',xmflops2,'  time(s):',t2-t1,gosa
C
      stop
      END
C
C
C**************************************************************
      subroutine initmt
C**************************************************************
      IMPLICIT REAL*4(a-h,o-z)
C
      PARAMETER (mimax=513,mjmax=257,mkmax=257)
C      PARAMETER (mimax=257,mjmax=129,mkmax=129)
C      PARAMETER (mimax=129,mjmax=65,mkmax=65)
C      PARAMETER (mimax=65,mjmax=33,mkmax=33)
C


CC Arrey
      common /pres/ p(mimax,mjmax,mkmax)
      common /mtrx/ a(mimax,mjmax,mkmax,4),
     +     b(mimax,mjmax,mkmax,3),c(mimax,mjmax,mkmax,3)
      common /bound/ bnd(mimax,mjmax,mkmax)
      common /work/ wrk1(mimax,mjmax,mkmax),wrk2(mimax,mjmax,mkmax)
CC other constants
      common /others/ imax,jmax,kmax,omega
!HPF$ DISTRIBUTE (*,*,BLOCK):: p,bnd,wrk1,wrk2
!HPF$ DISTRIBUTE (*,*,BLOCK,*) :: a,b,c
C
      do k=1,mkmax
         do j=1,mjmax
            do i=1,mimax
               a(i,j,k,1)=0.0
               a(i,j,k,2)=0.0
               a(i,j,k,3)=0.0
               a(i,j,k,4)=0.0
               b(i,j,k,1)=0.0
               b(i,j,k,2)=0.0
               b(i,j,k,3)=0.0
               c(i,j,k,1)=0.0
               c(i,j,k,2)=0.0
               c(i,j,k,3)=0.0
               p(i,j,k) =0.0
               wrk1(i,j,k)=0.0
               bnd(i,j,k)=0.0
            enddo
         enddo
      enddo
C
      do k=1,kmax
         do j=1,jmax
            do i=1,imax
               a(i,j,k,1)=1.0
               a(i,j,k,2)=1.0
               a(i,j,k,3)=1.0
               a(i,j,k,4)=1.0/6.0
               b(i,j,k,1)=0.0
               b(i,j,k,2)=0.0
               b(i,j,k,3)=0.0
               c(i,j,k,1)=1.0
               c(i,j,k,2)=1.0
               c(i,j,k,3)=1.0
               p(i,j,k) =float((k-1)*(k-1))/float((kmax-1)*(kmax-1))
               wrk1(i,j,k)=0.0
               bnd(i,j,k)=1.0
            enddo
         enddo
      enddo
C
      return
      end
C
C*************************************************************
      subroutine jacobi(nn,gosa)
C*************************************************************
      IMPLICIT REAL*4(a-h,o-z)
C
      PARAMETER (mimax=513,mjmax=257,mkmax=257)
C      PARAMETER (mimax=257,mjmax=129,mkmax=129)
C      PARAMETER (mimax=129,mjmax=65,mkmax=65)
C      PARAMETER (mimax=65,mjmax=33,mkmax=33)
C

CC Arrey
      common /pres/ p(mimax,mjmax,mkmax)
      common /mtrx/ a(mimax,mjmax,mkmax,4),
     +     b(mimax,mjmax,mkmax,3),c(mimax,mjmax,mkmax,3)
      common /bound/ bnd(mimax,mjmax,mkmax)
      common /work/ wrk1(mimax,mjmax,mkmax),wrk2(mimax,mjmax,mkmax)
CC other constants
      common /others/ imax,jmax,kmax,omega
!HPF$ DISTRIBUTE (*,*,BLOCK):: p,bnd,wrk1,wrk2
!HPF$ DISTRIBUTE (*,*,BLOCK,*) :: a,b,c
C
C
      DO loop=1,nn
!HPF$ INDEPENDENT,REDUCTION(gosa)
         DO K=2,kmax-1
            DO J=2,jmax-1
               DO 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)
     2                 +b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K)
     3                 -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)
     5                 -p(I,J+1,K-1)+p(I,J-1,K-1))
     6                 +b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1)
     7                 -p(I+1,J,K-1)+p(I-1,J,K-1))
     8                 +c(I,J,K,1)*p(I-1,J,K)+c(I,J,K,2)*p(I,J-1,K)
     9                 +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
               enddo
            enddo
         enddo
C
         DO K=2,kmax-1
            DO J=2,jmax-1
               DO I=2,imax-1
                  p(I,J,K)=wrk2(I,J,K)
               enddo
            enddo
         enddo
C
      enddo
CC End of iteration
      return
      end
