サイト内の現在位置

NumPy互換数値演算ライブラリNLCPyの新機能SCAインターフェース

no.020ステンシル定義式における係数の取り扱いやステンシル計算式の定義などについて2022.8.1

1. NLCPy SCAインターフェースの使用方法
 前回は「ライフシミュレーション」プログラムをサンプルとして、descriptorやkernelを使ったステンシル計算について説明しました。今回はステンシル定義式の係数の取り扱いや、ループ処理を使ったステンシル計算の定義など、前回触れなかった内容を説明します。

2. ステンシル定義式における係数の取り扱い
 以下のステンシル定義式を例に係数について説明します。

dx(0)=c1×dx(-1)-c2×dx(+1)
係数c1とc2において、Figure-1の入力[3]ではステンシル定義式にスカラ値を直接指定しています。この他、係数をndarray配列に収納して式内で参照することもできます(入力[4])。


fig_1

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倍された係数が計算に反映されます。


fig_2

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を定義しなおすことで、新しい係数を反映したステンシル計算が実行されます。


fig_3

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を再定義することで新しい係数を使ったステンシル計算が行われます。


fig_4

Figure-4


 Figure-5における入力[2]のように係数をcoef=nlcpy.arange(10.) により定義して、入力[3]で記述子dx[0]との積をとるといった使い方もあります。この場合も、係数に対する変更は入力[4]のように直ちに次のステンシル計算に反映されます。


fig_5

Figure-5


 3. ループを使用したステンシル記述定義
 Figure-6のようにループを使用したdesc定義を行う方法があります。

desc=nlcpy.sca.empty_description()
その際は上記のように、先にempty_description()を使って空の記述定義を行います。次いで、ループの中でdesc+=dx[i]を定義します。ループ終了後にSCA kernelを生成します。多くの場合、Pythonのループ処理は低速ですが、このループ処理は記述定義を目的としたもので、ステンシル計算は低速なループ処理と無関係にベクトルエンジンで一括処理されます。
 また、ループ処理によるステンシル記述定義は高次元な配列に対しても実施できます。Figre-6の例では、3x3の2次元平面のステンシル計算を行っています。


fig_6

Figure-6


 4. 計算結果出力先におけるオフセット調整
 SCA kernelを生成する際に、ステンシル計算結果が入るdesc_oに対してオフセットを指定する方法があります。Figure-7に示す入力[2]のように、desc_o=dxout[1]とオフセットする値を指定することで、プラス方向に1オフセットされた結果が返されます。反対に入力[3]のようにdesc_o=dxout[-1]とすると、マイナス方向へオフセットされます。


fig_7

Figure-7


 1より大きい袖領域を確保するには、SCA kernelを生成時に、Figure-8に示す入力[2]のようにdesc_o=dxout[2]と指定します。このとき、dxout[b]のbはdesc=dx[a]で指定したa以下(b<=a)の値を選択します。これは袖領域の大きさdesc=dx[a]を超えたオフセットはエラーとなるためです。


fig_8

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]で同じステンシル記述式を使用しています。


fig_9