Japan
サイト内の現在位置
NumPy互換数値演算ライブラリNLCPyの新機能SCAインターフェース
no.020ステンシル定義式における係数の取り扱いやステンシル計算式の定義などについて2022.8.11. NLCPy SCAインターフェースの使用方法
前回は「ライフシミュレーション」プログラムをサンプルとして、descriptorやkernelを使ったステンシル計算について説明しました。今回はステンシル定義式の係数の取り扱いや、ループ処理を使ったステンシル計算の定義など、前回触れなかった内容を説明します。
2. ステンシル定義式における係数の取り扱い
以下のステンシル定義式を例に係数について説明します。
dx(0)=c1×dx(-1)-c2×dx(+1)
係数c1とc2において、Figure-1の入力[3]ではステンシル定義式にスカラ値を直接指定しています。この他、係数をndarray配列に収納して式内で参照することもできます（入力[4]）。 In [1]:
import nlcpy as vp
import numpy as np
In [2]:
x = vp.arange(10.)
dx = vp.sca.create_descriptor(x)
print(x)
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
In [3]:
desc = 0.5*dx[-1] - 1.5*dx[1] # スカラ値による係数定義
kern = vp.sca.create_kernel(desc)
kern.execute()
Out[3]:
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
In [4]:
coef_vp = vp.array([0.5, 1.5]) # NLCPy ndarray配列を使用した係数定義
desc = coef_vp[0]*dx[-1] - coef_vp[1]*dx[1]
kern = vp.sca.create_kernel(desc)
kern.execute()
Out[4]:
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
In [5]:
vp.sca.destroy_kernel(kern)
Figure-1
なお、NLCPyとNumPyに入れた係数に対して行う変更は異なる結果を返す場合があります。
係数c1及びc2にNLCPyのndarrayを使用する場合には、ステンシル定義式
desc=coef[0]×dx[-1]-coef[1]×dx[1]
の係数の参照先はndarrayのデータ値そのものではなく、ndarray配列のアドレスになります。アドレスを参照先にすることで、NLCPy ndarray配列の係数に対する変更は、SCA kernelを再度定義することなく変更後の係数を使ったkernelの実行が行われます。Figure-2の入力[3]では、NLCPy ndarrayに格納した係数を2倍して、その直後にkernelを実行しています。kernelを改めて定義せずとも出力[3]では2倍された係数が計算に反映されます。 In [1]:
import nlcpy as vp
import numpy as np
In [2]:
x = vp.arange(10.)
dx = vp.sca.create_descriptor(x)
print(x)
coef_vp = vp.array([0.5, 1.5]) # NLCPyのndarray配列に係数代入
desc = coef_vp[0]*dx[-1] - coef_vp[1]*dx[1]
kern = vp.sca.create_kernel(desc)
kern.execute()
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
Out[2]:
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
In [3]:
coef_vp *= 2.0 # 係数配列に対する変更
kern.execute()
Out[3]:
array([ 0., -6., -8., -10., -12., -14., -16., -18., -20., 0.])
In [4]:
vp.sca.destroy_kernel(kern)
Figure-2
係数の変更後にkernelの定義が不要であることから、一度定義したSCA kernelを繰り返し利用することができます。ステンシル計算結果を次の計算に繰り返し反映させるような処理において、その都度SCA kernelを作成することは不必要なリソース消費につながります。SCA kernelを繰り返し使うことによって無駄を省くことができます。
NLCPyの代わりにNumPyのndarrayへ係数を収納しSCA kernel実行した例がFigure-3です。Figure-2と同様に入力[2]でdescを定義した後にkernelを実行しています。入力[3]で係数を2倍してkernelを再度実行していますが出力[3]には新しい係数が反映されていません。NumPyのndarray配列に入った係数はsca.create_kernel(desc)を実行する際に一時的にNLCPyのndarrayに変換されます。この例ではsca.create_kernel(desc)を行っていないため、入力[3]で行ったNumPyへの係数変更はステンシル計算の結果に反映されません。入力[4]のように改めてSCA kernelを定義しなおすことで、新しい係数を反映したステンシル計算が実行されます。
In [1]:
import nlcpy as vp
import numpy as np
In [2]:
x = vp.arange(10.)
dx = vp.sca.create_descriptor(x)
print(x)
coef_np = np.array([0.5, 1.5]) # NumPyのndarray配列に係数代入
desc = coef_np[0]*dx[-1] - coef_np[1]*dx[1]
kern = vp.sca.create_kernel(desc)
kern.execute()
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
Out[2]:
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
In [3]:
coef_np *= 2.0 # 係数配列に対する変更
kern.execute()
Out[3]:
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
In [4]:
desc = coef_np[0] * dx[-1] - coef_np[1] * dx[1] # ステンシル定義式とkernel再定義
kern = vp.sca.create_kernel(desc)
kern.execute()
Out[4]:
array([ 0., -6., -8., -10., -12., -14., -16., -18., -20., 0.])
In [5]:
vp.sca.destroy_kernel(kern)
Figure-3
NLCPy ndarrayの係数においても、次のように
desc=1.0×coef[0]×dx[-1]-1.0×coef[1]×dx[1]
スカラ値（この例では1.0）がステンシル定義式に含まれていると係数に対する変更は直ちに実行結果に反映されません。これは1.0xcoef[0]および1.0xcoef[1]がsca.create_kernel（）実行時に一時的なメモリに保管されてステンシル計算が行われるためです。Figure-4における入力[4]のように、SCA kernelを再定義することで新しい係数を使ったステンシル計算が行われます。 In [1]:
import nlcpy as vp
import numpy as np
In [2]:
x = vp.arange(10.)
coef_vp = vp.array([0.5, 1.5])
dx = vp.sca.create_descriptor(x)
desc = 1.0*coef_vp[0]*dx[-1] - 1.0*coef_vp[1]*dx[1] # スカラ値(1.0)が含まれる定義式
kern = vp.sca.create_kernel(desc)
kern.execute()
Out[2]:
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
In [3]:
coef_vp *= 2
kern.execute() # coefに対する変更は結果に反映されない
Out[3]:
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
In [4]:
desc = 1.0*coef_vp[0]*dx[-1] - 1.0*coef_vp[1]*dx[1]
kern = vp.sca.create_kernel(desc) # sca.create_kernel()を再度実行することで新しいcoefが結果に反映
kern.execute()
Out[4]:
array([ 0., -6., -8., -10., -12., -14., -16., -18., -20., 0.])
In [5]:
vp.sca.destroy_kernel(kern)
Figure-4
Figure-5における入力[2]のように係数をcoef=nlcpy.arange(10.) により定義して、入力[3]で記述子dx[0]との積をとるといった使い方もあります。この場合も、係数に対する変更は入力[4]のように直ちに次のステンシル計算に反映されます。
In [1]:
import nlcpy as vp
import numpy as np
In [2]:
x = vp.arange(10.)
dx = vp.sca.create_descriptor(x)
print(x)
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
In [3]:
coef = vp.arange(10.) # NLCPy arange()を使用した係数の用意
desc = coef * dx[0]
kern = vp.sca.create_kernel(desc)
kern.execute()
Out[3]:
array([ 0., 1., 4., 9., 16., 25., 36., 49., 64., 81.])
In [4]:
coef *= 2.0 # 係数に対する変更は直ちに結果に反映
kern.execute()
Out[4]:
array([ 0., 2., 8., 18., 32., 50., 72., 98., 128., 162.])
In [5]:
vp.sca.destroy_kernel(kern)
Figure-5
3. ループを使用したステンシル記述定義
Figure-6のようにループを使用したdesc定義を行う方法があります。
desc=nlcpy.sca.empty_description()
その際は上記のように、先にempty_description()を使って空の記述定義を行います。次いで、ループの中でdesc+=dx[i]を定義します。ループ終了後にSCA kernelを生成します。多くの場合、Pythonのループ処理は低速ですが、このループ処理は記述定義を目的としたもので、ステンシル計算は低速なループ処理と無関係にベクトルエンジンで一括処理されます。また、ループ処理によるステンシル記述定義は高次元な配列に対しても実施できます。Figre-6の例では、3x3の2次元平面のステンシル計算を行っています。
In [1]:
import nlcpy as vp
In [2]:
x = vp.arange(25, dtype='f4').reshape(5,5)
dx = vp.sca.create_descriptor(x)
desc = vp.sca.empty_description()
for i in range(-1, 2):
for j in range(-1, 2):
desc += dx[i, j] # ループ内でステンシル式を定義
kern = vp.sca.create_kernel(desc)
kern.execute()
Out[2]:
array([[ 0., 0., 0., 0., 0.],
[ 0., 54., 63., 72., 0.],
[ 0., 99., 108., 117., 0.],
[ 0., 144., 153., 162., 0.],
[ 0., 0., 0., 0., 0.]], dtype=float32) In [3]:
vp.sca.destroy_kernel(kern)
Figure-6
4. 計算結果出力先におけるオフセット調整
SCA kernelを生成する際に、ステンシル計算結果が入るdesc_oに対してオフセットを指定する方法があります。Figure-7に示す入力[2]のように、desc_o=dxout[1]とオフセットする値を指定することで、プラス方向に１オフセットされた結果が返されます。反対に入力[3]のようにdesc_o=dxout[-1]とすると、マイナス方向へオフセットされます。
In [1]:
import nlcpy as vp
In [2]:
x = vp.arange(10.)
print(x)
xout = vp.zeros_like(x)
dx, dxout = vp.sca.create_descriptor((x, xout))
desc = dx[-1] + dx[1]
kern = vp.sca.create_kernel(desc, desc_o=dxout[1]) # プラス方向にオフセット
kern.execute()
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
Out[2]:
array([ 0., 0., 2., 4., 6., 8., 10., 12., 14., 16.])
In [3]:
kern = vp.sca.create_kernel(desc, desc_o=dxout[-1]) # マイナス方向にオフセット
kern.execute()
Out[3]:
array([ 2., 4., 6., 8., 10., 12., 14., 16., 14., 16.])
In [4]:
vp.sca.destroy_kernel(kern)
Figure-7
１より大きい袖領域を確保するには、SCA kernelを生成時に、Figure-8に示す入力[2]のようにdesc_o=dxout[2]と指定します。このとき、dxout[b]のbはdesc=dx[a]で指定したa以下(b<=a)の値を選択します。これは袖領域の大きさdesc=dx[a]を超えたオフセットはエラーとなるためです。
In [1]:
import nlcpy as vp
In [2]:
x = vp.arange(1.0,10.)
print(x)
xout = vp.zeros_like(x)
dx, dxout = vp.sca.create_descriptor((x, xout))
desc = dx[2] # 袖領域の確保
kern = vp.sca.create_kernel(desc, desc_o=dxout[2]) # 2 オフセット
kern.execute()
[1. 2. 3. 4. 5. 6. 7. 8. 9.]
Out[2]:
array([0., 0., 3., 4., 5., 6., 7., 8., 9.])
In [3]:
x = vp.arange(1.0,10.)
print(x)
xout = vp.zeros_like(x)
dx, dxout = vp.sca.create_descriptor((x, xout))
kern = vp.sca.create_kernel(desc, desc_o=dxout[1]) # 1 オフセット
kern.execute()
[1. 2. 3. 4. 5. 6. 7. 8. 9.]
Out[3]:
array([0., 3., 4., 5., 6., 7., 8., 9., 0.])
In [4]:
x = vp.arange(1.0,10.)
print(x)
xout = vp.zeros_like(x)
dx, dxout = vp.sca.create_descriptor((x, xout))
kern = vp.sca.create_kernel(desc, desc_o=dxout[0]) # オフセット無し
kern.execute()
[1. 2. 3. 4. 5. 6. 7. 8. 9.]
Out[4]:
array([3., 4., 5., 6., 7., 8., 9., 0., 0.])
In [5]:
vp.sca.destroy_kernel(kern)
Figure-8
5. Simplified Indexingを使用したステンシル記述の汎用化
Ellipsis ...を使ったインデックス簡素化による汎用的なステンシル記述を行うことも可能です。例えば、ステンシル記述を次のように書くことにより
desc=dx[...,-1]+dx[...,1]
4次元以下の多次元配列に対して同じ記述を使ってステンシル計算が実施されます。 1次元dx[-1]、2次元dx[0,-1]、3次元dx[0,0,-1]、そして4次元ではdx[0,0,0,-1]という各次元での個別の記述を、dx[...,-1]一つで対応します。Figure-9の例では、1次元データに対しての処理を入力[3]、そして2次元データでは入力[5]で同じステンシル記述式を使用しています。