tugas3-openMP/tugas.f
2022-11-27 13:18:52 +07:00

249 lines
7.6 KiB
Fortran

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 -------------------
! use portlib
use omp_lib
IMPLICIT REAL*4(a-h,o-z)
real*8 t1,t2
C
C PARAMETER (mimax=513,mjmax=257,mkmax=257)
C PARAMETER (mimax=257,mjmax=129,mkmax=129)
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
dimension time0(2),time1(2)
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=10000
write(*,*) ' Start rehearsal measurement process.'
write(*,*) ' Measure the performance in 10000 times.'
C
! cpu0=dtime(time0)
t1 = omp_get_wtime()
C
C Jacobi iteration
call jacobi(nn,gosa)
C
! cpu1= dtime(time1)
t2 = omp_get_wtime()
! cpu = cpu1
cpu = t2-t1
flop=real(kmax-2)*real(jmax-2)*real(imax-2)*34.0*real(nn)
xmflops2=flop/cpu*1.0e-6
write(*,*) ' MFLOPS:',xmflops2,' time(s):',cpu,gosa
C
C end the test loop
! nn=ifix(ttarget/(cpu/3.0))
! write(*,*) 'Now, start the actual measurement process.'
! write(*,*) 'The loop will be excuted in',nn,' times.'
! write(*,*) 'This will take about one minute.'
! write(*,*) 'Wait for a while.'
C
C Jacobi iteration
! cpu0=dtime(time0)
! call jacobi(nn,gosa)
C
! cpu1= dtime(time1)
! cpu = cpu1
! flop=real(kmax-2)*real(jmax-2)*real(imax-2)*34.0*real(nn)
! xmflops2=flop*1.0e-6/cpu
C
CCC xmflops2=nflop/cpu*1.0e-6*float(nn)
C
! write(*,*) ' Loop executed for ',nn,' times'
! write(*,*) ' Gosa :',gosa
! write(*,*) ' MFLOPS:',xmflops2, ' time(s):',cpu
! score=xmflops2/82.84
! write(*,*) ' Score based on Pentium III 600MHz :',score
C
! pause
stop
END
C
C
C**************************************************************
subroutine initmt
C**************************************************************
IMPLICIT REAL*4(a-h,o-z)
C
C PARAMETER (mimax=513,mjmax=257,mkmax=257)
C PARAMETER (mimax=257,mjmax=129,mkmax=129)
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
C
!$OMP do
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
!$OMP end do
C
!$OMP do
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
!$OMP end do
C
return
end
C
C*************************************************************
subroutine jacobi(nn,gosa)
C*************************************************************
IMPLICIT REAL*4(a-h,o-z)
C
C PARAMETER (mimax=513,mjmax=257,mkmax=257)
C PARAMETER (mimax=257,mjmax=129,mkmax=129)
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
C
C
DO loop=1,nn
gosa=0.0
!$OMP do 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
!$OMP end do
C
!$OMP do
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
!$OMP end do
C
enddo
CC End of iteration
return
end