Richmond のモーメント法解析におけるベクトル化率の向上
東北大学大学院工学研究科
陳

強

澤谷 邦男

1. はじめに
Richmond のモーメント法[1]はアンテナや電磁界の数値解析によく用いられる手法である．
この手法は分割セグメント数が少なくても計算精度がよいというすぐれた特徴[2]を持っているが，大型のアレーアンテナや平面アンテナ解析などの大規模な数値解析において CPU時間が長いという問題点がある．そのため，大規模な数値解析に対して高速な計算を行えるアルゴリズムの開発が必要である．スーパーコンピュータはベクトル演算の機能[3]を持っており，ベクトル化率の高いプログラムに対して高速な計算が行える．本報告ではスーパーコンピュータのベクトル演算機能を利用して Richmond のモーメント法による数値解析を高速に行うためのベクトル化率の高いアルゴリズムを開発した結果を述べる．
2.

高いベクトル化率のアルゴリズムの開発
Richmond のモーメント法は，基底関数と重み関数として正弦状関数を用いた手法である．

そのプログラムのフローチャートを図 1に示す．まずアンテナの形状をダイポールセグメントに分割して，各セグメントの位置情報を計算しておく．その後位置情報に基づいて，各セグメント間の自己，相互インピーダンスを計算してインピーダンス行列を生成する．さらにインピーダンス行列の逆行列を求めることにより，アンテナ上の電流分布，入力インピーダンス，放射電磁界などを求めることができる[1]．
モーメント法による数値解析の一連の流れで，計算時間の大分部はダイポールセグメント間のインピーダンス行列を求めるための積分計算と，インピーダンス行列の逆行列計算で占められている．逆行列計算は東北大学大型計算機センターにスーパーコンピュータ用のベクトル化率の高いサブルーチンが用意されているので，それを用いることにより高速化が図ることができる．プログラム全体の計算時間の高速化を図るには，いかにインピーダンス行列を生成するサブルーチンに対してベクトル化率を向上させて，高速化するかが問題となってくる．
一般的にプログラムのベクトル化率を上げるにはプログラム中にベクトル演算対象となるDO ループ文を多くつくることと，DO ループの振り返し数（ループ長）を大きくすることが考えられる．繰り返し数を大きくする方法としては･多重 DO ループは最深の DO ループのループ長が大きくなるようにループを入れかえる．
･小さなループ長の DO ループは分離させる．
･多重 DO ループを一重の DO ループに変換する．という手法がある．

START
アンテナ形状をダイポールセグメントに分割し，座標を配列に代入
DO
I=1,N
J=I,N（Ｎ:セグメント数)

Z行列要素を求める
座標変換
積分計算
RETURN

Z行列の逆行列を計算する
電流分布
入力インピーダンス
放射電磁界を計算する
STOP

図 1

Richmond のモーメント法のフローチャート

図 2に従来のモーメント法のアルゴリズムを示す．外側のループは分割セグメントの相互インピーダンスを計算するものであり，計算の規模が大きいほど，ループ長が大きい．それに対して，内側のループはインピーダンス行列要素を計算する積分分点数であり，4〜6 ぐらいの短いものである．外側の大きなループはベクトル演算の対称とはならず，最深の小さなループに対してのみベクトル化されるのでプログラム全体のベクトル化率はわずか 10〜30%である．そのため，ベクトル化率が上がるような新しいアルゴリズムの開発が必要である．そこで，図 3に示すように，まず繰り返し数の小さい積分計算のループを分離し，積分計算に必要なセグメント間の位置情報をあらかじめ計算して配列に入れ，その後，ループ間の配列の受け渡しにより大きな繰り返し数のループ内で積分計算を行うようにした．さらに，積分計算のループの繰り返し数を大きくするために，二重ループを一重ループに変換した．
これによりインピーダンス行列要素を求める部分のベクトル化率が 99%以上に向上し，プログラム全体のベクトル化率も 99%以上になった．

N×Nのインピーダンス行列を計算するサブルーチン
N:セグメント数（大きい数）

M:積分分点数（小さい数）

DO I=1,N
DO J=I,N[座標変換]
DO K=1,M[積分計算]ENDDO
Z(I,J)=Z
ENDDO
ENDDO

図 2

ベクトル演算対象外ループベクトル演算対象ループ

従来のアルゴリズム

N×Nのインピーダンス行列を計算するサブルーチン
N:セグメント数（大きい数）
DO I=1,N
DO J=I,N[座標変換](変換した位置情報を配列に入れる)
ENDDO
ENDDO

ベクトル演算対象ループ

DO K=1,N(N+1)/2[積分計算]
Z(K)=Z
ENDDO

ベクトル演算対象ループ

図 3

開発したアルゴリズム

分かりやすく説明するために，図 4と図 5に９行×９列のインピーダンス行列の計算例を示す．白い部分は計算されていない行列要素を表し，黒く塗られた部分は計算された行列要素を表す．従来のアルゴリズムにおけるインピーダンス行列計算では，図 4に示すように行列の各要素を一つづつ計算している．一方，開発したアルゴリズムにおけるインピーダンス計算では，プログラムのベクトル化によって 10 倍高速になったとすると図 5に示すように一命令で 10 個要素を計算できる．つまり，従来のアルゴリズム計算ではインピーダンス行列を計算するのに 15 命令必要だったのに，開発したアルゴリズムによる計算ではわずか 2命令で計算できることがわかる．
Richmond のモーメント法のプログラムがベクトル化によって高速化を実現できるとすると，多くのセグメント数で分割された大型のアンテナ解析でもインピーダンス行列が少ないベクトル計算命令で求められ，モーメント法によるアンテナ解析の計算時間を大幅に減少することが期待できる．

Z１１ Z12 Z13 Z14 Z15

Z１１ Z12 Z13 Z14 Z15

Z１１ Z12 Z13 Z14 Z15Z21 Z22 Z23 Z24 Z25

Z21 Z22 Z23 Z24 Z25

Z21 Z22 Z23 Z24 Z25

Z31 Z32 Z33 Z34 Z35

Z31 Z32 Z33 Z34 Z35

Z31 Z32 Z33 Z34 Z35

Z41 Z42 Z43 Z44 Z45

Z41 Z42 Z43 Z44 Z45

Z41 Z42 Z43 Z44 Z45

Z51 Z52 Z53 Z54 Z55

Z51 Z52 Z53 Z54 Z55

Z51 Z52 Z53 Z54 Z55

Z１１ Z12 Z13 Z14 Z15

Z１１ Z12 Z13 Z14 Z15

Z21 Z22 Z23 Z24 Z25

Z21 Z22 Z23 Z24 Z25

Z31 Z32 Z33 Z34 Z35

・・ ・・ ・・

Z31 Z32 Z33 Z34 Z35

Z41 Z42 Z43 Z44 Z45

Z41 Z42 Z43 Z44 Z45

Z51 Z52 Z53 Z54 Z55

Z51 Z52 Z53 Z54 Z55

図 4従来のアルゴリズムにおけるインピーダンス行列計算

Z１１ Z12 Z13 Z14 Z15

Z１１ Z12 Z13 Z14 Z15

Z１１ Z12 Z13 Z14 Z15

Z21 Z22 Z23 Z24 Z25

Z21 Z22 Z23 Z24 Z25

Z21 Z22 Z23 Z24 Z25

Z31 Z32 Z33 Z34 Z35

Z31 Z32 Z33 Z34 Z35

Z31 Z32 Z33 Z34 Z35

Z41 Z42 Z43 Z44 Z45

Z41 Z42 Z43 Z44 Z45

Z41 Z42 Z43 Z44 Z45

Z51 Z52 Z53 Z54 Z55

Z51 Z52 Z53 Z54 Z55

Z51 Z52 Z53 Z54 Z55

図 53.

開発したアルゴリズムにおけるインピーダンス行列計算

数値解析結果従来のアルゴリズムのプログラムを WS と SX で実行し，計算時間の比較を行った．その結果を図 6に示す．横軸がアンテナのダイポールセグメント分割数 N,縦軸が CPU TIME[sec]を表す．図から WS の計算時間が SX の約 10 分の 1 であることがわかる．この結果から，SX を使用して計算の高速化を図るためには高いベクトル化率のプログラムが不可欠でCPU Time [sec]あることがわかる．

10310210110010-110-210-310-410SXWS0101102

Number of Segment
図 6103N

従来のアルゴリズムにおける WS と SX での計算時間比較(WS:Alpha21164A-600MHz，SX:Super Computer SX4 in Tohoku University)

開発したアルゴリズムと従来のアルゴリズムのプログラムを SX で実行し，インピーダンス行列計算時間の比較を行った．その結果を図 7に示す．横軸がアンテナのダイポールセグ
メント分割数 N,縦軸が CPU TIME [sec]を表す．この図から，開発したアルゴリズムの方が従来のアルゴリズムより約 40 倍計算時間が高速化できたことがわかる．

103[sec]

102101CPU Time

10010-110-210-310-410OldNew0101102

Number of Segment

103N

図 7 開発したアルゴリズムと従来のものによるインピーダンス行列計算に要する CPU Time の比較(Old： Old Algorithm，N=601 の時ベクトル化率：0.0%，
New：New Algorithm，N=601 の時ベクトル化率:：99.69%．)

次に開発したアルゴリズムと従来のアルゴリズムのプログラムを SX で実行し，逆行列計算時間を含めたプログラム全体の計算時間の比較を行った．その結果を図 8に示す．横軸がアンテナのダイポールセグメント分割数 N,縦軸が CPU TIME [sec]を表す．この図から，開発したアルゴリズムの方が従来のアルゴリズムより約 50 倍計算時間が高速化できたことがわかる．

CPU Time [sec]10101010321010-110-210-310-410Old

New
0101102

Number of Segment
図 8103N開発したアルゴリズムと従来のものによる MoM の CPU Time の比較（Old：Old Algorithm，N=601 の時ベクトル化率:34.18%，New：New Algorithm，N=601 の時ベクトル化率:99.55%)

4.

まとめRichmond のモーメント法におけるインピーダンス行列を計算するサブルーチンに対して高いベクトル化率で計算するためのアルゴリズムを開発した．その結果，ベクトル化率は99%以上になり，インピーダンス行列の計算時間で従来のアルゴリズムに対して約 40 倍高速になった．また，逆行行列計算時間も含むプログラム全体の計算時間でも約 50 倍高速になった．このような解析結果から，今回開発したアルゴリズムはスーパーコンピュータのハードウェア性能を十分に引き出せたものと考えられる．今後，スーパーコンピュータのハードウェア性能がさらに向上することによって，今回開発したアルゴリズムはさらに高速な計算を実現できるものと考えている．
謝辞
本文は東北大学大学院工学研究科電気･通信工学専攻古屋聡士氏（平成 12 年博士前期課程修了）が在学中に行った研究をまとめたものである．
本研究の一部は東北大学大型計算機センターのスーパーコンピュータを利用し，同センターとの共同研究で行われたものである．また，研究にあたっては同センターの有益な指導と多大な協力をいただいた．
参考文献[1] J. H. RICHMOND and N. H. GEARY, “Mutual Impedance of Nonplanar-Skew
Sinusoidal Dipoles”, IEEE Trans. Antennas Propagat., Vol. AP-23, No. 3, pp. 412-414,
May 1975.
[2] 澤谷, “モーメント法によるアンテナ解析中級コース”, 電子情報通信学会第 12 回アンテナ・伝搬における設計・解析手法ワークショップ, 1998.
[3] 吉田，FORTRAN90/SX の自動ベクトル化について，東北大学大型計算機センター広報，
SENAC, Vol. 30, No. 4, pp. 27-42, 1997.