1!*********************************************************************
2!
3! This benchmark test program is measuring a cpu performance
4! of floating point operation by a Poisson equation solver.
5!!
6! If you have any question, please ask me via email.
7! written by Ryutaro HIMENO, November 26, 2001.
8! Version 3.0
9! ----------------------------------------------
10! Ryutaro Himeno, Dr. of Eng.
11! Head of Computer Information Division,
12! RIKEN (The Institute of Pysical and Chemical Research)
13! Email : himeno@postman.riken.go.jp
14! ---------------------------------------------------------------
15! You can adjust the size of this benchmark code to fit your target
16! computer. In that case, please chose following sets of
17! (mimax,mjmax,mkmax):
18! small : 65,33,33
19! small : 129,65,65
20! midium: 257,129,129
21! large : 513,257,257
22! ext.large: 1025,513,513
23! This program is to measure a computer performance in MFLOPS
24! by using a kernel which appears in a linear solver of pressure
25! Poisson eq. which appears in an incompressible Navier-Stokes solver.
26! A point-Jacobi method is employed in this solver as this method can
27! be easyly vectrized and be parallelized.
28! ------------------
29! Finite-difference method, curvilinear coodinate system
30! Vectorizable and parallelizable on each grid point
31! No. of grid points : imax x jmax x kmax including boundaries
32! ------------------
33! A,B,C:coefficient matrix, wrk1: source term of Poisson equation
34! wrk2 : working area, OMEGA : relaxation parameter
35! BND:control variable for boundaries and objects ( = 0 or 1)
36! P: pressure
37! -------------------
38! -------------------
39! "use portlib" statement on the next line is for Visual fortran
40! to use UNIX libraries. Please remove it if your system is UNIX.
41! -------------------
42! use portlib
43!
44module pres
45  real(4),dimension(:,:,:),allocatable :: p
46end module pres
47!
48module mtrx
49  real(4),dimension(:,:,:,:),allocatable :: a,b,c
50end module mtrx
51!
52module bound
53  real(4),dimension(:,:,:),allocatable :: bnd
54end module bound
55!
56module work
57  real(4),dimension(:,:,:),allocatable :: wrk1,wrk2
58end module work
59!
60module others
61  integer :: mimax,mjmax,mkmax
62  integer :: imax,jmax,kmax
63  real(4),parameter :: omega=0.8
64end module others
65!
66!
67program HimenoBMTxp_F90
68  use others
69!
70  implicit none
71!
72  integer :: nn
73  integer :: ic,icr,icm
74  real(4) :: flop,xmflops2,score,gosa
75  real(8) :: cpu0,cpu1,cpu,dt
76  ! ttarget specifys the measuring period in sec
77  real(4),parameter :: ttarget=60.0
78  real(8),external :: second
79!
80  call readparam
81!
82!! Initializing matrixes
83  call initmem
84  call initmt
85  print *,' mimax=',mimax,' mjmax=',mjmax,' mkmax=',mkmax
86  print *,' imax=',imax,' jmax=',jmax,' kmax=',kmax
87!
88  call system_clock(ic,icr,icm)
89  dt= 1.0/real(icr,8)
90  print '(2x,a,e10.5)','Time measurement accuracy : ',dt
91!! Start measuring
92!
93  nn=3
94  print *,' Start rehearsal measurement process.'
95  print *,' Measure the performance in 3 times.'
96!
97  cpu0=second()
98!! Jacobi iteration
99  call jacobi(nn,gosa)
100!
101  cpu1=second()
102  cpu = (cpu1 - cpu0)*dt
103  flop=real(kmax-2,4)*real(jmax-2,4)*real(imax-2,4)*34.0
104  if(cpu /= 0.0)  xmflops2=flop/real(cpu,4)*1.0e-6*real(nn,4)
105  print *,'  MFLOPS:',xmflops2,'  time(s):',cpu,gosa
106!
107!! end the test loop
108  nn=int(ttarget/(cpu/3.0))
109  print *,'Now, start the actual measurement process.'
110  print *,'The loop will be excuted in',nn,' times.'
111  print *,'This will take about one minute.'
112  print *,'Wait for a while.'
113!
114!! Jacobi iteration
115  cpu0= second()
116!! Jacobi iteration
117  call jacobi(nn,gosa)
118!
119  cpu1= second()
120  cpu = (cpu1 - cpu0)*dt
121  if(cpu /= 0.0)  xmflops2=flop*1.0e-6/real(cpu,4)*real(nn,4)
122!
123!!!      xmflops2=nflop/cpu*1.0e-6*float(nn)
124!
125  print *,' Loop executed for ',nn,' times'
126  print *,' Gosa :',gosa
127  print *,' MFLOPS:',xmflops2, '  time(s):',cpu
128  score=xmflops2/82.84
129  print *,' Score based on Pentium III 600MHz :',score
130!
131  pause
132!   stop
133end program HimenoBMTxp_F90
134!
135!
136subroutine readparam
137!
138  use others
139!
140  implicit none
141!
142  character(10) :: size
143!
144  print *,'Select Grid-size:'
145!  print *,'For example:'
146  print *,'Grid-size= '
147  print *,'           XS (64x32x32)'
148  print *,'           S  (128x64x64)'
149  print *,'           M  (256x128x128)'
150  print *,'           L  (512x256x256)'
151  print *,'           XL (1024x512x512)'
152!  print *,' Grid-size = '
153  read(*,*) size
154  call grid_set(size)
155!
156  imax=mimax-1
157  jmax=mjmax-1
158  kmax=mkmax-1
159!
160  return
161end subroutine readparam
162!
163!
164subroutine grid_set(size)
165!
166  use others
167!
168  implicit none
169!
170  character(10),intent(in) :: size
171!
172  select case(size)
173  case("xs")
174     mimax=65
175     mjmax=33
176     mkmax=33
177  case("XS")
178     mimax=65
179     mjmax=33
180     mkmax=33
181  case("s")
182     mimax=129
183     mjmax=65
184     mkmax=65
185  case("S")
186     mimax=129
187     mjmax=65
188     mkmax=65
189  case("m")
190     mimax=257
191     mjmax=129
192     mkmax=129
193  case("M")
194     mimax=257
195     mjmax=129
196     mkmax=129
197  case("l")
198     mimax=513
199     mjmax=257
200     mkmax=257
201  case("L")
202     mimax=513
203     mjmax=257
204     mkmax=257
205  case("xl")
206     mimax=1025
207     mjmax=513
208     mkmax=513
209  case("XL")
210     mimax=1025
211     mjmax=513
212     mkmax=513
213  case default
214     print *,'Invalid input character !!'
215     stop
216  end select
217!
218  return
219end subroutine grid_set
220!
221!
222!**************************************************************
223subroutine initmt
224!**************************************************************
225  use pres
226  use mtrx
227  use bound
228  use work
229  use others
230!
231  implicit none
232!
233  integer :: i,j,k
234!
235  a=0.0
236  b=0.0
237  c=0.0
238  p=0.0
239  wrk1=0.0
240  wrk2=0.0
241  bnd=0.0
242!
243  a(1:imax,1:jmax,1:kmax,1:3)=1.0
244  a(1:imax,1:jmax,1:kmax,4)=1.0/6.0
245  c(1:imax,1:jmax,1:kmax,:)=1.0
246  bnd(1:imax,1:jmax,1:kmax)=1.0
247  do k=1,kmax
248     p(:,:,k)=real((k-1)*(k-1),4)/real((kmax-1)*(kmax-1),4)
249  enddo
250!
251  return
252end subroutine initmt
253!
254!*************************************************************
255subroutine initmem
256!*************************************************************
257  use pres
258  use mtrx
259  use bound
260  use work
261  use others
262!
263  implicit none
264!
265  allocate(p(mimax,mjmax,mkmax))
266  allocate(a(mimax,mjmax,mkmax,4),b(mimax,mjmax,mkmax,3), &
267       c(mimax,mjmax,mkmax,3))
268  allocate(bnd(mimax,mjmax,mkmax))
269  allocate(wrk1(mimax,mjmax,mkmax),wrk2(mimax,mjmax,mkmax))
270!
271  return
272end subroutine initmem
273!
274!*************************************************************
275subroutine jacobi(nn,gosa)
276!*************************************************************
277  use pres
278  use mtrx
279  use bound
280  use work
281  use others
282!
283  implicit none
284!
285  integer,intent(in) :: nn
286  real(4),intent(inout) :: gosa
287  integer :: i,j,k,loop
288  real(4) :: s0,ss
289!
290!
291  do loop=1,nn
292     gosa= 0.0
293     do k=2,kmax-1
294        do j=2,jmax-1
295           do i=2,imax-1
296              s0=a(I,J,K,1)*p(I+1,J,K) &
297                   +a(I,J,K,2)*p(I,J+1,K) &
298                   +a(I,J,K,3)*p(I,J,K+1) &
299                   +b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K) &
300                               -p(I-1,J+1,K)+p(I-1,J-1,K)) &
301                   +b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1) &
302                               -p(I,J+1,K-1)+p(I,J-1,K-1)) &
303                   +b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1) &
304                               -p(I+1,J,K-1)+p(I-1,J,K-1)) &
305                   +c(I,J,K,1)*p(I-1,J,K) &
306                   +c(I,J,K,2)*p(I,J-1,K) &
307                   +c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
308              ss=(s0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
309              GOSA=GOSA+SS*SS
310              wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
311           enddo
312        enddo
313     enddo
314!
315     p(2:imax-1,2:jmax-1,2:kmax-1)= &
316          wrk2(2:imax-1,2:jmax-1,2:kmax-1)
317!
318  enddo
319!! End of iteration
320  return
321end subroutine jacobi
322!
323!
324!
325function second() result(rtime)
326!
327  implicit none
328!
329  integer :: ic,ir,im
330  real(8) :: rtime
331!
332  call system_clock(ic,ir,im)
333!
334  rtime= real(ic,8)
335!
336end function second
337