姫野ベンチマークをPythonで実行し、[numpy, numba]で高速化。
姫野ベンチは、3次元ポアソン方程式をヤコビ法で解く反復処理の速度を図るベンチマークです。
本家のホームページは→こちら Fortran90とCの公式版があります。
20数年前からある古いベンチマークですが、流体解析に必要なメモリ律速のベンチマークに向いています。
Pythonを始めるにあたり、どの程度の速度がでるか興味があって始めました。
ベンチマーク環境
- CPU: Intel(R) Core(TM) i5-3570K CPU @ 3.40GHz(Boost時3.8GHz) (4core)
- RAM: 32GB
- OS: CentOS Linux release 7.5.1804 (64bit)
- Fortran compiler: gfortran 4.8.5 20150623 (Red Hat 4.8.5-28) (GCC)
- C compiler: gcc 4.8.5 20150623 (Red Hat 4.8.5-28) (GCC)
- Python: Python 3.6.3 (default, May 13 2018, 22:01:23) (Anaconda 5.1.0)
かなり古いCPUです。1coreでの計測です。
オリジナル版の計測結果
使用したソースとcompile option
- gfortran -Ofast himenoBMTxp.f90 コードのダウンロード(Fortran90)
- gcc -Ofast himenoBMTxpa.c コードのダウンロードの（C, dynamic allocate version）
ベンチマーク結果
| 単精度 | XS | S | M | L | XL |
|---|---|---|---|---|---|
| Fortran[Mflops] | 4011 | 3845 | 2591 | 3040 | 2881 |
| C [Mflops] | 1808 | 1790 | 1250 | 791 | 407 |
| 配列サイズ(IxJxK) | 65x33x33 | 129x65x65 | 257x129x129 | 513x257x257 | 1025x513x513 |
| メモリ使用量 | 5MB | 31MB | 235MB | 1,854MB | 14,752MB |
| メモリ使用量は、Fortran版で/usr/timeの出力。C版は10%弱使用量が少ない。 | |||||
| Fortranは安定して速い。Cは配列サイズが大きいと速度低下が大きい。 |
Python版
| 単精度 | XS | S | M | L | XL |
|---|---|---|---|---|---|
| pure_python | 1.9 | 2.1 | 2.1 | 2.1 | 2.1 |
| numba order='C' | 1281 | 611 | 483 | 420 | 368 |
| numba order='F' | 2216 | 1955 | 1558 | 1161 | 870 |
| numpy order='C' | 1043 | 954 | 941 | 929 | 715 |
| numpy order='F' | 1186 | 1180 | 1099 | 1078 | 821 |
ソースは、後で示します。試したのは、
- pure_python: Fortran版をなるべく変えないようPythonに移植したもの。配列はnumpy.ndarrayを使用
- numba : 上記のpure_python版に@numba.jit(nopython=True)を追加し、エラーが出たところを必要最小限で修正したもの。(nopython=True)は、numbaが失敗した際にpythonにデグレードせずエラーを出させるため。
- numpy : pure_pnuython版のループ部をnumpyの添え字をずらしたviewを与えることで、一行にまとめたもの
※:numba版とnumpy版はndarray作成時のorderを'C'と'F'の両方実施
実行方法
python himenoBMT_pure.py [XS|S|M|L|XL]サイズを指定してください
憶えたてのf-stringsを使っているので、 python3.6以上です。
結果(単精度)
- pure_pythonは、目も当てられない。
- numbaはCとほぼ互角の速度。ソースコードの変更も少なく、forループをそのまま書けば良いので、見通しが良く、コードのメンテナンスも容易。Fortranに対しても半分程度の速度が出るので、予想以上の結果。ifortと比較したいところ。
- numpyでも結構行ける。配列サイズが大きい場合には、numbaよりも速い場合がある。
- orderは'C'よりも'F'の方が速い。これは元のコードのループ順に依存すると思われる。
結論
PythonでもCに匹敵する速度で実行できるコードが書けそう
numbaはすごい、でも@numba.jit(nopython=True)で動かないと戸惑う
以下、ソースコード
まずは、pure_python版、なるべくオリジナルのFortran版を素直に移植したつもりだが、
Classを作る必要があったのか…。Python始めて2カ月の初心者ということで大目に見て下さい。
himenoBMT_pure.py
'''
Created on 2018/06/10
@author: akio
'''
import sys
import time
import numpy as np
class pres:
pres = []
class mtrx:
a = []
b = []
c = []
class bound:
bnd = []
class work:
wrk1 = []
wrk2 = []
class others:
imax = -1
jmax = -1
kmax = -1
mimax = -1
mjmax = -1
mkmax = -1
omega = 0.8
def readparam(size):
grid_set(size)
others.imax = others.mimax-1
others.jmax = others.mjmax-1
others.kmax = others.mkmax-1
return
def grid_set(size):
if size == "xs":
mimax = 65
mjmax = 33
mkmax = 33
elif size == "s":
mimax = 129
mjmax = 65
mkmax = 65
elif size == "m":
mimax=257
mjmax=129
mkmax=129
elif size == "l":
mimax=513
mjmax=257
mkmax=257
elif size == "xl":
mimax=1025
mjmax=513
mkmax=513
else:
print("Invalid input character !! {}".format(size))
sys.exit("-1")
others.mimax = mimax
others.mjmax = mjmax
others.mkmax = mkmax
def initmt():
pres.p[:,:,:] = 0
mtrx.a[:,:,:,:] = 0
mtrx.b[:,:,:,:] = 0
mtrx.c[:,:,:] = 0
bound.bnd[:,:,:] = 0.0
work.wrk1[:,:,:] = 0.0
work.wrk2[:,:,:] = 0.0
mtrx.a[0:others.imax,0:others.jmax,0:others.kmax,0:3] = 1.0
mtrx.a[0:others.imax,0:others.jmax,0:others.kmax,3] = 1.0/6.0
mtrx.c[0:others.imax,0:others.jmax,0:others.kmax,:] = 1.0
bound.bnd[0:others.imax,0:others.jmax,0:others.kmax] = 1.0
for k in range(others.kmax):
pres.p[:,:,k] = ((k-1)*(k-1))/((others.kmax-1)*(others.kmax-1))
def initmem():
dim = (others.mimax, others.mjmax, others.mkmax)
dim4 = (others.mimax, others.mjmax, others.mkmax,4)
dim3 = (others.mimax, others.mjmax, others.mkmax,3)
pres.p = np.ndarray(dim, dtype=np.float32)
mtrx.a = np.ndarray(dim4, dtype=np.float32)
mtrx.b = np.ndarray(dim3, dtype=np.float32)
mtrx.c = np.ndarray(dim3, dtype=np.float32)
bound.bnd = np.ndarray(dim, dtype=np.float32)
work.wrk1 = np.ndarray(dim, dtype=np.float32)
work.wrk2 = np.ndarray(dim, dtype=np.float32)
def jacobi(nn):
#gosa = np.array(0.0, dtype=np.float32)
for loop in range(nn):
gosa = np.zeros((1), dtype=np.float32)
#ss = 0.0 #np.array(0.0, np.float32)
for k in range(1,others.kmax-1):
for j in range(1,others.jmax-1):
for i in range(1,others.imax-1):
s0 = mtrx.a[i,j,k,0]*pres.p[i+1,j,k] \
+mtrx.a[i,j,k,1]*pres.p[i,j+1,k] \
+mtrx.a[i,j,k,2]*pres.p[i,j,k+1] \
+mtrx.b[i,j,k,0]*(pres.p[i+1,j+1,k]-pres.p[i+1,j-1,k] \
-pres.p[i-1,j+1,k]+pres.p[i-1,j-1,k]) \
+mtrx.b[i,j,k,1]*(pres.p[i,j+1,k+1]-pres.p[i,j-1,k+1] \
-pres.p[i,j+1,k-1]+pres.p[i,j-1,k-1]) \
+mtrx.b[i,j,k,2]*(pres.p[i+1,j,k+1]-pres.p[i-1,j,k+1] \
-pres.p[i+1,j,k-1]+pres.p[i-1,j,k-1]) \
+mtrx.c[i,j,k,0]*pres.p[i-1,j,k] \
+mtrx.c[i,j,k,1]*pres.p[i,j-1,k] \
+mtrx.c[i,j,k,2]*pres.p[i,j,k-1]+work.wrk1[i,j,k]
ss=(s0*mtrx.a[i,j,k,3]-pres.p[i,j,k])*bound.bnd[i,j,k]
gosa[0] = gosa[0] + ss*ss
work.wrk2[i,j,k]=pres.p[i,j,k]+ others.omega*ss
if loop < 4: print('loop: {} gosa: {}'.format(loop,gosa[0]))
pres.p[1:others.imax-1,1:others.jmax-1,1:others.kmax-1] = \
work.wrk2[1:others.imax-1,1:others.jmax-1,1:others.kmax-1]
return float(gosa[0])
def second():
return time.time()
if __name__ == '__main__':
ttarget = 10.0 # 60
xmflops2 = 0.0
gosa = 0.0
if len(sys.argv) > 1:
size = sys.argv[1].lower()
else:
size = 's'
readparam(size)
initmem()
initmt()
print(' mimax={} mjmax={} mkmax={}'
.format(others.mimax,others.mjmax,others.mkmax))
print(' imax={} jmax={} kmax={}'
.format(others.imax,others.jmax,others.kmax))
dt= time.get_clock_info('time').resolution
print(' {:s} , {:12.5e}'.format('Time measurement accuracy : ',dt))
## Start measuring
nn=3
print(' Start rehearsal measurement process.')
print(' Measure the performance in 3 times.')
cpu0=second()
## Jacobi iteration
gosa = jacobi(nn)
cpu1=second()
cpu = (cpu1 - cpu0)
flop=(others.kmax-2)*(others.jmax-2)*(others.imax-2)*34.0
if (cpu != 0.0): xmflops2 = flop/cpu*1.0e-6*nn
print(' MFLOPS:{} time(s):{}'
.format(xmflops2,cpu))
#
## end the test loop
nn=int(ttarget/(cpu/3.0))+1
print(' Gosa : {}'.format(gosa))
print('Now, start the actual measurement process.')
print('The loop will be excuted in {} times.'.format(nn))
print('This will take about one minute.')
print('Wait for a while.')
#
## Jacobi iteration
cpu0= second()
## Jacobi iteration
gosa = jacobi(nn)
#
cpu1= second()
cpu = (cpu1 - cpu0)
if (cpu != 0.0): xmflops2 = flop*1.0e-6/cpu*nn
#
### xmflops2=nflop/cpu*1.0e-6*float(nn)
#
print('The loop will be excuted in {} times.'.format(nn))
print(' Gosa : {}'.format(gosa))
print(' MFLOPS:{} time(s):{}'
.format(xmflops2,cpu))
score=xmflops2/82.84
print(' Score based on Pentium III 600MHz : {}'.format(score))
#
# pause
# stop
numbaバージョン
次は、numba版。pure_python版に@numba.jitを付けただけだと、実行時にjacobi関数の中で、pres.pやmtrx.aは知らない！と怒られたので、全て引数で渡すように変更した。結果では、order='C'よりもorder='F'の方が速かったが、ループ順をオリジナルのkji順からijk順に変更するとorder='C'の方が速い。実は、オリジナルのkjiでorder='F'よりも、ijk順でorder='C'の方が1割ほど速かった（複雑になるので、数字は載せていません）。
himenoBMT_numba.py
'''
Created on 2018/06/10
@author: akio
'''
import sys
import time
import numpy as np
import numba
class pres:
pres = []
class mtrx:
a = []
b = []
c = []
class bound:
bnd = []
class work:
wrk1 = []
wrk2 = []
class others:
imax = -1
jmax = -1
kmax = -1
mimax = -1
mjmax = -1
mkmax = -1
omega = 0.8
def readparam(size):
grid_set(size)
others.imax = others.mimax-1
others.jmax = others.mjmax-1
others.kmax = others.mkmax-1
return
def grid_set(size):
if size == "xs":
mimax = 65
mjmax = 33
mkmax = 33
elif size == "s":
mimax = 129
mjmax = 65
mkmax = 65
elif size == "m":
mimax=257
mjmax=129
mkmax=129
elif size == "l":
mimax=513
mjmax=257
mkmax=257
elif size == "xl":
mimax=1025
mjmax=513
mkmax=513
else:
print("Invalid input character !! {}".format(size))
sys.exit("-1")
others.mimax = mimax
others.mjmax = mjmax
others.mkmax = mkmax
def initmt():
pres.p[:,:,:] = 0
mtrx.a[:,:,:,:] = 0
mtrx.b[:,:,:,:] = 0
mtrx.c[:,:,:] = 0
bound.bnd[:,:,:] = 0.0
work.wrk1[:,:,:] = 0.0
work.wrk2[:,:,:] = 0.0
mtrx.a[0:others.imax,0:others.jmax,0:others.kmax,0:3] = 1.0
mtrx.a[0:others.imax,0:others.jmax,0:others.kmax,3] = 1.0/6.0
mtrx.c[0:others.imax,0:others.jmax,0:others.kmax,:] = 1.0
bound.bnd[0:others.imax,0:others.jmax,0:others.kmax] = 1.0
for k in range(others.kmax):
pres.p[:,:,k] = ((k-1)*(k-1))/((others.kmax-1)*(others.kmax-1))
def initmem():
dim = (others.mimax, others.mjmax, others.mkmax)
dim4 = (others.mimax, others.mjmax, others.mkmax,4)
dim3 = (others.mimax, others.mjmax, others.mkmax,3)
atype = np.float32
aorder = 'F'
pres.p = np.ndarray(dim, dtype=atype, order=aorder)
mtrx.a = np.ndarray(dim4, dtype=atype, order=aorder)
mtrx.b = np.ndarray(dim3, dtype=atype, order=aorder)
mtrx.c = np.ndarray(dim3, dtype=atype, order=aorder)
bound.bnd = np.ndarray(dim, dtype=atype, order=aorder)
work.wrk1 = np.ndarray(dim, dtype=atype, order=aorder)
work.wrk2 = np.ndarray(dim, dtype=atype, order=aorder)
#@numba.jit('(i8,f8,i8,i8,i8,i8,i8,i8,f8,f4[:,:,:],f4[:,:,:,:],f4[:,:,:,:],f4[:,:,:,:],f4[:,:,:],f4[:,:,:],f4[:,:,:])',nopython=True)
@numba.jit(nopython=True)
def jacobi(nn,mimax,mjmax,mkmax,imax,jmax,kmax,omega \
,p,a,b,c,bnd,wrk1,wrk2):
gosa = np.zeros((1), dtype=np.float32)
for loop in range(nn):
gosa[0] = 0.0
#ss = np.array(0.0, np.float32)
for k in range(1,kmax-1):
for j in range(1,jmax-1):
for i in range(1,imax-1):
s0 = a[i,j,k,0]*p[i+1,j,k] \
+a[i,j,k,1]*p[i,j+1,k] \
+a[i,j,k,2]*p[i,j,k+1] \
+b[i,j,k,0]*(p[i+1,j+1,k]-p[i+1,j-1,k] \
-p[i-1,j+1,k]+p[i-1,j-1,k]) \
+b[i,j,k,1]*(p[i,j+1,k+1]-p[i,j-1,k+1] \
-p[i,j+1,k-1]+p[i,j-1,k-1]) \
+b[i,j,k,2]*(p[i+1,j,k+1]-p[i-1,j,k+1] \
-p[i+1,j,k-1]+p[i-1,j,k-1]) \
+c[i,j,k,0]*p[i-1,j,k] \
+c[i,j,k,1]*p[i,j-1,k] \
+c[i,j,k,2]*p[i,j,k-1]+wrk1[i,j,k]
ss=(s0*a[i,j,k,3]-p[i,j,k])*bnd[i,j,k]
gosa[0] = gosa[0]+ss*ss
wrk2[i,j,k]=p[i,j,k]+ omega*ss
#if loop < 4: print(f'loop: {loop} gosa: {gosa}')
#print('loop: {} gosa: {}'.format(loop,gosa))
#print(gosa)
p[1:imax-1,1:jmax-1,1:kmax-1] = \
wrk2[1:imax-1,1:jmax-1,1:kmax-1]
return float(gosa[0])
def second():
return time.time()
if __name__ == '__main__':
ttarget = 60.0 # 60
xmflops2 = 0.0
gosa = 0.0
if len(sys.argv) > 1:
size = sys.argv[1].lower()
else:
size = 's'
readparam(size)
initmem()
initmt()
''' numba jit'''
nn=1
gosa = jacobi(nn,others.mimax,others.mjmax,others.mkmax \
,others.imax,others.jmax,others.kmax,others.omega \
,pres.p,mtrx.a,mtrx.b,mtrx.c,bound.bnd,work.wrk1,work.wrk2)
initmt()
print(f' mimax={others.mimax} mjmax={others.mjmax} mkmax={others.mkmax}')
print(f' imax={others.imax} jmax={others.jmax} kmax={others.kmax}')
dt= time.get_clock_info('time').resolution
print(' %a , %10.5e' % ('Time measurement accuracy : ',dt))
## Start measuring
nn=3
print(' Start rehearsal measurement process.')
print(' Measure the performance in 3 times.')
cpu0=second()
## Jacobi iteration
gosa = jacobi(nn,others.mimax,others.mjmax,others.mkmax \
,others.imax,others.jmax,others.kmax,others.omega \
,pres.p,mtrx.a,mtrx.b,mtrx.c,bound.bnd,work.wrk1,work.wrk2)
cpu1=second()
cpu = (cpu1 - cpu0)
flop=(others.kmax-2)*(others.jmax-2)*(others.imax-2)*34.0
if (cpu != 0.0): xmflops2 = flop/cpu*1.0e-6*nn
print(f' MFLOPS:{xmflops2} time(s):{cpu} {gosa}')
#
## end the test loop
nn=int(ttarget/(cpu/3.0))+1
print('Now, start the actual measurement process.')
print(f'The loop will be excuted in {nn} times.')
print('This will take about one minute.')
print('Wait for a while.')
#
## Jacobi iteration
cpu0= second()
## Jacobi iteration
gosa = jacobi(nn,others.mimax,others.mjmax,others.mkmax \
,others.imax,others.jmax,others.kmax,others.omega \
,pres.p,mtrx.a,mtrx.b,mtrx.c,bound.bnd,work.wrk1,work.wrk2)
#
cpu1= second()
cpu = (cpu1 - cpu0)
if (cpu != 0.0): xmflops2 = flop*1.0e-6/cpu*nn
#
### xmflops2=nflop/cpu*1.0e-6*float(nn)
#
print(f' Loop executed for {nn} times')
print(f' Gosa : {gosa}')
print(f' MFLOPS: {xmflops2} time(s):{cpu}')
score=xmflops2/82.84
print(f' Score based on Pentium III 600MHz : {score}')
#
# pause
# stop
numpyバージョン
下にソースのコア部分の抜粋を示します。コメントを参照してください。
def jacobi(nn,mimax,mjmax,mkmax,imax,jmax,kmax,omega \
,p,a,b,c,bnd,wrk1,wrk2,s0,ss):
p111 = p[1:imax-1,1:jmax-1,1:kmax-1]
p211 = p[2:imax ,1:jmax-1,1:kmax-1]
p121 = p[1:imax-1,2:jmax ,1:kmax-1]
p112 = p[1:imax-1,1:jmax-1,2:kmax ]
'''
中略。
上記のように、添字をずらしたndarrayのviewを作り、
'''
for loop in range(nn):
```
k,j,iのループを無くして、
numpy配列の演算に置き換えた。
※注意事項（はまったところ）
s0[:,:,:]=...の "[:,:,:]"を忘れて、s0=とすると、
s0を計算結果の新しいndarrayで上書きしてしまうので注意。
```
s0111[:,:,:] = a0*p211 \
+a1* p121 \
+a2* p112 \
+b0*( p221 - p201 \
-p021 + p001) \
+b1*( p122 - p102 \
-p120 + p100) \
+b2*( p212 - p012 \
-p210 + p010) \
+c0*p011 \
+c1*p101 \
+c2*p110 + w1
ss111[:,:,:] = (s0111*a3-p111)*bnd111
以下、numpy版の全体です。
himenoBMT_numpy2.py
'''
Created on 2018/06/10
@author: akio
'''
import sys
import time
from ctypes import c_float
import numpy as np
import numba
class pres:
pres = []
class mtrx:
a = []
b = []
c = []
class bound:
bnd = []
class work:
wrk1 = []
wrk2 = []
class others:
imax = -1
jmax = -1
kmax = -1
mimax = -1
mjmax = -1
mkmax = -1
omega = 0.8
def readparam(size):
grid_set(size)
others.imax = others.mimax-1
others.jmax = others.mjmax-1
others.kmax = others.mkmax-1
return
def grid_set(size):
if size == "xs":
mimax = 65
mjmax = 33
mkmax = 33
elif size == "s":
mimax = 129
mjmax = 65
mkmax = 65
elif size == "m":
mimax=257
mjmax=129
mkmax=129
elif size == "l":
mimax=513
mjmax=257
mkmax=257
elif size == "xl":
mimax=1025
mjmax=513
mkmax=513
else:
print("Invalid input character !! {}".format(size))
sys.exit("-1")
others.mimax = mimax
others.mjmax = mjmax
others.mkmax = mkmax
def initmt():
pres.p[:,:,:] = 0
mtrx.a[:,:,:,:] = 0
mtrx.b[:,:,:,:] = 0
mtrx.c[:,:,:] = 0
bound.bnd[:,:,:] = 0.0
work.wrk1[:,:,:] = 0.0
work.wrk2[:,:,:] = 0.0
work.s0[:,:,:] = 0.0
work.ss[:,:,:] = 0.0
mtrx.a[0:others.imax,0:others.jmax,0:others.kmax,0:3] = 1.0
mtrx.a[0:others.imax,0:others.jmax,0:others.kmax,3] = 1.0/6.0
mtrx.c[0:others.imax,0:others.jmax,0:others.kmax,:] = 1.0
bound.bnd[0:others.imax,0:others.jmax,0:others.kmax] = 1.0
for k in range(others.kmax):
pres.p[:,:,k] = ((k-1)*(k-1))/((others.kmax-1)*(others.kmax-1))
def initmem():
dim = (others.mimax, others.mjmax, others.mkmax)
dim4 = (others.mimax, others.mjmax, others.mkmax,4)
dim3 = (others.mimax, others.mjmax, others.mkmax,3)
atype = np.float32
aorder = 'F'
pres.p = np.ndarray(dim, dtype=atype, order=aorder)
mtrx.a = np.ndarray(dim4, dtype=atype, order=aorder)
mtrx.b = np.ndarray(dim3, dtype=atype, order=aorder)
mtrx.c = np.ndarray(dim3, dtype=atype, order=aorder)
bound.bnd = np.ndarray(dim, dtype=atype, order=aorder)
work.wrk1 = np.ndarray(dim, dtype=atype, order=aorder)
work.wrk2 = np.ndarray(dim, dtype=atype, order=aorder)
work.s0 = np.ndarray(dim, dtype=atype, order=aorder)
work.ss = np.ndarray(dim, dtype=atype, order=aorder)
#@numba.jit('(i8,f8,i8,i8,i8,i8,i8,i8,f8,f4[:,:,:],f4[:,:,:,:],f4[:,:,:,:],f4[:,:,:,:],f4[:,:,:],f4[:,:,:],f4[:,:,:])',nopython=True)
#@numba.jit(nopython=True)
def jacobi(nn,mimax,mjmax,mkmax,imax,jmax,kmax,omega \
,p,a,b,c,bnd,wrk1,wrk2,s0,ss):
#gosa = np.zeros((1), dtype=np.float32)
a0 = a[1:imax-1,1:jmax-1,1:kmax-1,0]
a1 = a[1:imax-1,1:jmax-1,1:kmax-1,1]
a2 = a[1:imax-1,1:jmax-1,1:kmax-1,2]
a3 = a[1:imax-1,1:jmax-1,1:kmax-1,3]
b0 = b[1:imax-1,1:jmax-1,1:kmax-1,0]
b1 = b[1:imax-1,1:jmax-1,1:kmax-1,1]
b2 = b[1:imax-1,1:jmax-1,1:kmax-1,2]
c0 = c[1:imax-1,1:jmax-1,1:kmax-1,0]
c1 = c[1:imax-1,1:jmax-1,1:kmax-1,1]
c2 = c[1:imax-1,1:jmax-1,1:kmax-1,2]
p111 = p[1:imax-1,1:jmax-1,1:kmax-1]
p211 = p[2:imax ,1:jmax-1,1:kmax-1]
p121 = p[1:imax-1,2:jmax ,1:kmax-1]
p112 = p[1:imax-1,1:jmax-1,2:kmax ]
p221 = p[2:imax ,2:jmax ,1:kmax-1]
p201 = p[2:imax ,0:jmax-2,1:kmax-1]
p021 = p[0:imax-2,2:jmax ,1:kmax-1]
p001 = p[0:imax-2,0:jmax-2,1:kmax-1]
p122 = p[1:imax-1,2:jmax ,2:kmax ]
p102 = p[1:imax-1,0:jmax-2,2:kmax ]
p120 = p[1:imax-1,2:jmax ,0:kmax-2]
p100 = p[1:imax-1,0:jmax-2,0:kmax-2]
p212 = p[2:imax ,1:jmax-1,2:kmax ]
p012 = p[0:imax-2,1:jmax-1,2:kmax ]
p210 = p[2:imax ,1:jmax-1,0:kmax-2]
p010 = p[0:imax-2,1:jmax-1,0:kmax-2]
p011 = p[0:imax-2,1:jmax-1,1:kmax-1]
p101 = p[1:imax-1,0:jmax-2,1:kmax-1]
p110 = p[1:imax-1,1:jmax-1,0:kmax-2]
w1 = wrk1[1:imax-1,1:jmax-1,1:kmax-1]
w2 = wrk2[1:imax-1,1:jmax-1,1:kmax-1]
s0111 = s0[1:imax-1,1:jmax-1,1:kmax-1]
ss111 = ss[1:imax-1,1:jmax-1,1:kmax-1]
bnd111 = bnd[1:imax-1,1:jmax-1,1:kmax-1]
for loop in range(nn):
s0111[:,:,:] = a0*p211 \
+a1* p121 \
+a2* p112 \
+b0*( p221 - p201 \
-p021 + p001) \
+b1*( p122 - p102 \
-p120 + p100) \
+b2*( p212 - p012 \
-p210 + p010) \
+c0*p011 \
+c1*p101 \
+c2*p110 + w1
ss111[:,:,:] = (s0111*a3-p111)*bnd111
w2[:,:,:] = p111+ omega*ss111
ss111[:,:,:] = ss111*ss111
gosa = float(np.sum(ss111, dtype=np.float32))
#if loop < 4: print(f'loop: {loop} gosa: {gosa}')
#print('loop: {} gosa: {}'.format(loop,gosa))
if loop < 3: print(gosa)
p111[:,:,:] = \
w2
return float(gosa)
def second():
return time.time()
if __name__ == '__main__':
ttarget = 60.0 # 60
xmflops2 = 0.0
gosa = 0.0
if len(sys.argv) > 1:
size = sys.argv[1].lower()
else:
size = 's'
readparam(size)
initmem()
initmt()
print(f' mimax={others.mimax} mjmax={others.mjmax} mkmax={others.mkmax}')
print(f' imax={others.imax} jmax={others.jmax} kmax={others.kmax}')
dt= time.get_clock_info('time').resolution
print(' %a , %10.5e' % ('Time measurement accuracy : ',dt))
## Start measuring
nn=3
print(' Start rehearsal measurement process.')
print(' Measure the performance in 3 times.')
cpu0=second()
## Jacobi iteration
gosa = jacobi(nn,others.mimax,others.mjmax,others.mkmax \
,others.imax,others.jmax,others.kmax,others.omega \
,pres.p,mtrx.a,mtrx.b,mtrx.c,bound.bnd \
,work.wrk1,work.wrk2,work.s0,work.ss)
cpu1=second()
cpu = (cpu1 - cpu0)
flop=(others.kmax-2)*(others.jmax-2)*(others.imax-2)*34.0
if (cpu != 0.0): xmflops2 = flop/cpu*1.0e-6*nn
print(f' MFLOPS:{xmflops2} time(s):{cpu} {gosa}')
#
## end the test loop
nn=int(ttarget/(cpu/3.0))+1
print('Now, start the actual measurement process.')
print(f'The loop will be excuted in {nn} times.')
print('This will take about one minute.')
print('Wait for a while.')
#
## Jacobi iteration
cpu0= second()
## Jacobi iteration
gosa = jacobi(nn,others.mimax,others.mjmax,others.mkmax \
,others.imax,others.jmax,others.kmax,others.omega \
,pres.p,mtrx.a,mtrx.b,mtrx.c,bound.bnd \
,work.wrk1,work.wrk2,work.s0,work.ss)
#
cpu1= second()
cpu = (cpu1 - cpu0)
if (cpu != 0.0): xmflops2 = flop*1.0e-6/cpu*nn
#
### xmflops2=nflop/cpu*1.0e-6*float(nn)
#
print(f' Loop executed for {nn} times')
print(f' Gosa : {gosa}')
print(f' MFLOPS: {xmflops2} time(s):{cpu}')
score=xmflops2/82.84
print(f' Score based on Pentium III 600MHz : {score}')
#
# pause
# stop
後書き
ソースコードを添付すると無駄に長くなってしまい、読み難くすみません。