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