Japan
サイト内の現在位置
NumPy互換数値演算ライブラリNLCPyの新機能SCAインターフェース
no.020ステンシル定義式における係数の取り扱いやステンシル計算式の定義などについて2022.8.11. NLCPy SCAインターフェースの使用方法
前回は「ライフシミュレーション」プログラムをサンプルとして、descriptorやkernelを使ったステンシル計算について説明しました。今回はステンシル定義式の係数の取り扱いや、ループ処理を使ったステンシル計算の定義など、前回触れなかった内容を説明します。
2. ステンシル定義式における係数の取り扱い
以下のステンシル定義式を例に係数について説明します。
import nlcpy as vp
import numpy as np
x = vp.arange(10.)
dx = vp.sca.create_descriptor(x)
print(x)
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
desc = 0.5*dx[-1] - 1.5*dx[1] # スカラ値による係数定義
kern = vp.sca.create_kernel(desc)
kern.execute()
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
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()
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
vp.sca.destroy_kernel(kern)
Figure-1
なお、NLCPyとNumPyに入れた係数に対して行う変更は異なる結果を返す場合があります。
係数c1及びc2にNLCPyのndarrayを使用する場合には、ステンシル定義式
import nlcpy as vp
import numpy as np
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.]
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
coef_vp *= 2.0 # 係数配列に対する変更
kern.execute()
array([ 0., -6., -8., -10., -12., -14., -16., -18., -20., 0.])
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を定義しなおすことで、新しい係数を反映したステンシル計算が実行されます。
import nlcpy as vp
import numpy as np
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.]
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
coef_np *= 2.0 # 係数配列に対する変更
kern.execute()
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
desc = coef_np[0] * dx[-1] - coef_np[1] * dx[1] # ステンシル定義式とkernel再定義
kern = vp.sca.create_kernel(desc)
kern.execute()
array([ 0., -6., -8., -10., -12., -14., -16., -18., -20., 0.])
vp.sca.destroy_kernel(kern)
Figure-3
NLCPy ndarrayの係数においても、次のように
import nlcpy as vp
import numpy as np
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()
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
coef_vp *= 2
kern.execute() # coefに対する変更は結果に反映されない
array([ 0., -3., -4., -5., -6., -7., -8., -9., -10., 0.])
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()
array([ 0., -6., -8., -10., -12., -14., -16., -18., -20., 0.])
vp.sca.destroy_kernel(kern)
Figure-4
Figure-5における入力[2]のように係数をcoef=nlcpy.arange(10.) により定義して、入力[3]で記述子dx[0]との積をとるといった使い方もあります。この場合も、係数に対する変更は入力[4]のように直ちに次のステンシル計算に反映されます。
import nlcpy as vp
import numpy as np
x = vp.arange(10.)
dx = vp.sca.create_descriptor(x)
print(x)
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
coef = vp.arange(10.) # NLCPy arange()を使用した係数の用意
desc = coef * dx[0]
kern = vp.sca.create_kernel(desc)
kern.execute()
array([ 0., 1., 4., 9., 16., 25., 36., 49., 64., 81.])
coef *= 2.0 # 係数に対する変更は直ちに結果に反映
kern.execute()
array([ 0., 2., 8., 18., 32., 50., 72., 98., 128., 162.])
vp.sca.destroy_kernel(kern)
Figure-5
3. ループを使用したステンシル記述定義
Figure-6のようにループを使用したdesc定義を行う方法があります。
また、ループ処理によるステンシル記述定義は高次元な配列に対しても実施できます。Figre-6の例では、3x3の2次元平面のステンシル計算を行っています。
import nlcpy as vp
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()
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)
vp.sca.destroy_kernel(kern)
Figure-6
4. 計算結果出力先におけるオフセット調整
SCA kernelを生成する際に、ステンシル計算結果が入るdesc_oに対してオフセットを指定する方法があります。Figure-7に示す入力[2]のように、desc_o=dxout[1]とオフセットする値を指定することで、プラス方向に1オフセットされた結果が返されます。反対に入力[3]のようにdesc_o=dxout[-1]とすると、マイナス方向へオフセットされます。
import nlcpy as vp
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.]
array([ 0., 0., 2., 4., 6., 8., 10., 12., 14., 16.])
kern = vp.sca.create_kernel(desc, desc_o=dxout[-1]) # マイナス方向にオフセット
kern.execute()
array([ 2., 4., 6., 8., 10., 12., 14., 16., 14., 16.])
vp.sca.destroy_kernel(kern)
Figure-7
1より大きい袖領域を確保するには、SCA kernelを生成時に、Figure-8に示す入力[2]のようにdesc_o=dxout[2]と指定します。このとき、dxout[b]のbはdesc=dx[a]で指定したa以下(b<=a)の値を選択します。これは袖領域の大きさdesc=dx[a]を超えたオフセットはエラーとなるためです。
import nlcpy as vp
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.]
array([0., 0., 3., 4., 5., 6., 7., 8., 9.])
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.]
array([0., 3., 4., 5., 6., 7., 8., 9., 0.])
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.]
array([3., 4., 5., 6., 7., 8., 9., 0., 0.])
vp.sca.destroy_kernel(kern)
Figure-8
5. Simplified Indexingを使用したステンシル記述の汎用化
Ellipsis ...を使ったインデックス簡素化による汎用的なステンシル記述を行うことも可能です。例えば、ステンシル記述を次のように書くことにより
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]で同じステンシル記述式を使用しています。