再帰というのは、ここを読んでいる人が知らない、ということはないと思いま
すが、ある関数なりサブルーチンが自分自身を呼び出すことです。例えば、ツ
リーコードのような、データ構造自体が自分と同じ構造をもつ、という意味で
再帰的なものでは、その処理も再帰を使って書くのが自然です。
ツリーコードの基本的なステップは
ツリーのデータ構造を作る
ツリーノードの物理量(位置、質量等)を計算する
ツリー構築を使って各粒子への力を計算する
です。この中でもっとも単純な、物理量を計算するところを再帰を使って書い
た例を以下に示します。
void bhnode::set_cm_quantities()
{
int i;
cmpos = 0.0;
cmmass = 0.0;
if (isleaf){
bhparticle * bp = bpfirst;
for(i = 0; i < nparticle; i++){
real mchild = (bp+i)->get_rp()->get_mass();
cmpos += mchild*(bp+i)->get_rp()->get_pos();
cmmass += mchild;
}
}else{
for(i=0;i<8;i++){
if (child[i] != NULL){
child[i]->set_cm_quantities();
real mchild = child[i]->cmmass;
cmpos += mchild*child[i]->cmpos;
cmmass += mchild;
}
}
}
cmpos /= cmmass;
}
これは C++ を使っていて、3次元ベクトルをクラスにしているので、重心の計
算のために質量掛ける位置を積算するのが
cmpos += mchild*(bp+i)->get_rp()->get_pos();
ですむ、といったことがあって単純ですが、基本的な構造が
ノードの物理量の計算:
子供が全部粒子の時
粒子の位置、質量を計算
子供がノードの時
各子供について
物理量を計算 <---- ここが再帰
位置、質量を積算
位置と質量の積の積算を質量でわって、重心を求める
というもので、要するに、ノードの重心はその下のノードなり粒子の重心の重
心である、という定義を、素直にプログラムに表現しています。
これを C なのですが再帰でなくて書いてあるのが以下の例です。
void calc_com_tree (prm_struct *pprm, ptcl_struct *pptcl, tree_struct *ptree)
{
int nnl, nnd, nst, ilev, inl, jnl, ind, ipn, ip;
int *lhp, *lnp, *lnd, *lhn, *ltn, *lnl, *lst, *lpn, *lfg;
real invn, x0, y0, z0, x1, y1, z1, ds, l0, lbox, invtheta, rc;
real *ppos, *com, *rc2, *minpos, lln[MAXTREELEV], lil[MAXTREELEV];
ppos = pptcl->ppos;
lbox = ptree->lbox;
minpos=ptree->minpos;
lhp = ptree->lhp;
lnp = ptree->lnp;
lnd = ptree->lnd;
lhn = ptree->lhn;
ltn = ptree->ltn;
nnl = ptree->nnl;
nnd = ptree->nnd;
invtheta = 1.0 / pprm->theta;
MALLOC (com, real, nnd * NDIM); ptree->com = com;
for (ind=0; ind<nnd*NDIM; ind++) com[ind] = 0.0;
MALLOC (rc2, real, nnd); ptree->rc2 = rc2;
MALLOC (lnl, int, nnd);
MALLOC (lst, int, nnl * 3);
lfg = lst + nnl;
lpn = lfg + nnl;
/* push the root node to the stack */
nst = 0;
lst[nst] = 0;
lfg[nst] = FALSE;
lpn[nst] = NULLNODE;
nst++;
lnl[0] = 0;
/* main loop */
while (nst){
/* pop from stack */
nst--;
inl = lst[nst];
ipn = lpn[nst];
ind = lnd[inl];
if (lfg[nst]){
/* ascend tree */
if (ipn != NULLNODE){
com[ipn*NDIM] += com[ind*NDIM];
com[ipn*NDIM+1] += com[ind*NDIM+1];
com[ipn*NDIM+2] += com[ind*NDIM+2];
}
} else {
/* descend tree */
if (ind == NULLNODE){
ip = lhp[inl];
com[ipn*NDIM] += ppos[ip*NDIM];
com[ipn*NDIM+1] += ppos[ip*NDIM+1];
com[ipn*NDIM+2] += ppos[ip*NDIM+2];
} else {
if (ipn != NULLNODE) lnl[ind] = lnl[ipn] + 1;
lfg[nst++] = TRUE;
for (jnl=ltn[ind]; jnl>=lhn[ind]; jnl--){
lst[nst] = jnl;
lfg[nst] = FALSE;
lpn[nst] = ind;
nst++;
}
}
}
}
FREE_MEM (lst);
/* calculate cell size */
for (ilev=0; ilev<MAXTREELEV; ilev++){
l0 = lbox / ((real)( 1 << ilev));
lln[ilev] = l0;
lil[ilev] = 1.0 / l0;
}
/* calculate center of mass */
for (ind=0; ind<nnd; ind++){
invn = 1.0 / (real) lnp[ind];
com[ind*NDIM] *= invn;
com[ind*NDIM+1] *= invn;
com[ind*NDIM+2] *= invn;
}
/* calculate and center of node and critical radius */
for (ind=0; ind<nnd; ind++){
/* center of mass */
x0 = com[ind*NDIM] - minpos[0];
y0 = com[ind*NDIM+1] - minpos[1];
z0 = com[ind*NDIM+2] - minpos[2];
/* center of node */
ilev=lnl[ind];
x1 = lln[ilev]* (0.5+ (real)((int) (x0 * lil[ilev])));
y1 = lln[ilev]* (0.5+ (real)((int) (y0 * lil[ilev])));
z1 = lln[ilev]* (0.5+ (real)((int) (z0 * lil[ilev])));
/* check con is consistent with com of children (see tree.h) */
CHECK_COMCON;
/* distance between center of mass and center of node */
x1 -= x0;
y1 -= y0;
z1 -= z0;
ds = sqrt(x1 * x1 + y1 * y1 + z1 * z1);
rc = lln[ilev] * invtheta + ds;
rc2[ind] = rc * rc;
}
FREE_MEM (lnl);
}
すみません、何をしているかよくわからないので、解説とか変更の試みは今回
省略します。nst という変数があって、これがスタックポインタになっていて、
色々なものがスタックに積まれたり取り出されたりしているようですが、私の
理解力を超えています。再帰ができない Fortran で、ベクトル化できるよう
に書いたものの例を以下に示します。
subroutine hacCF2
#include "tdefs.F"
integer top(nbits), last(nbits), nlevel
integer i, j, k, r, lastnd, level
REAL dterm, mr, dx1, dx2, dx3
lastnd=nxtnode+ncell-1
do 200 i=nxtnode,lastnd
mass(i)=0.0
do 210 k=1,ndim
pos(i,k)=0.0
210 continue
200 continue
top(1)=nxtnode
last(1)=nxtnode
nlevel=1
do 230 i=nxtnode+1, lastnd
if(sizetol2(i) .lt. sizetol2(i-1)*0.9) then
nlevel=nlevel+1
top(nlevel)=i
last(nlevel-1)=i-1
endif
230 continue
last(nlevel)=lastnd
do 240 level=nlevel, 1, -1
do 250 j=1,nsubcell
do 260 i=top(level), last(level)
r=subp(i,j)
if(r .ne. null) then
mr=mass(r)
mass(i)=mass(i)+mr
pos(i,1)=pos(i,1)+pos(r,1)*mr
pos(i,2)=pos(i,2)+pos(r,2)*mr
pos(i,3)=pos(i,3)+pos(r,3)*mr
endif
260 continue
250 continue
do 280 i=top(level), last(level)
do 270 k=1,ndim
pos(i,k)=pos(i,k)/mass(i)
270 continue
280 continue
240 continue
end
ここでは、再帰によってツリーを辿るのをスタックを使ってエミュレーション
する、といった面倒なことはしないで、ツリーの階層の下から順番に作る、
同じレベルのノードの間でベクトル化する、という方法をとっています。ここ
では、そもそもノードがレベルの順番に並んでいて、上のレベルのノードは添
字が小さい、と仮定しています。このことは、ツリーの作り方によって保証さ
れています。
230 のループでは、レベルの境界を捜しています。これは無能なコードで、
の計算量で全ノードをみていますが、本当は2分探索を繰り返すことで
でできるし、そもそもツリーを作った時に憶えておけば計算
する必要もないのですが、まあ、20年以上前に書いたコードなので、、、
レベルの境界が決まれば、後は下から作るだけで、レベルのループが 240 で
す。250、260 は奇妙な順番で、 250 で8個ある子供についてのループを回し、
その内側の 260 で同じレベルにある全ノードを回るループにしていますが、
これは大抵のベクトル化 fortran コンパイラが最内側ループしかベクトル化
しないので自然な順番と入れ替えているものです。
再帰はアルゴリズムを簡潔で、従って間違いが少なくデバッグも容易な形に表
現する優れた方法ですが、並列化やベクトル化とは必ずしも適合しないことも
あります。特に、コンパイラの能力として最内側ループだけしか並列化(ベク
トル化)しないベクトル化コンパイラの場合には、仮に再帰が扱える言語やコ
ンパイラであったとし、また再帰の中に自然な並列性が表現されていたとして
も、再帰の形に書いたものがハードウェアで効率良く走る可能性は殆どありま
せん。
例えば、上の再帰の例では、8個の子供について重心を計算するのは完全に並
列に行えますし、さらにその8個の子供の64個くらいの子供についても並列に
行えるわけです。従って、MIMD で共有メモリ計算機用の良くできた並列化コ
ンパイラなら勝手に並列化してくれてもよいし、実際にデータフローマシン
用のコンパイラなら自然に並列化されます。
なお、例えば OpenMP でも、task や taskwait プラグマを使うことで並列実
行できるようです。但し、なにも考えずにこの辺を使うと、ノード数だけのタ
スクが生成されてしまうので実行効率は落ちるかもしれません。途中から並列
化しない再帰ルーチンを呼ぶといった工夫が必要になる可能性もあります。
ベクトル化と例えば OpenMP を使った並列化は、「どのように並列化されるべ
きか」という高レベルの概念は同じであっても、「どうやってコーディングす
るべきか」という観点では全く違ってきてしまうことが普通です。例えば、
上のベクトル化できる fortran ルーチンは、もちろん最内側ループを
omp parallel で並列化できますが、あまり効率は良くありません。これは、
このループが、ツリーのレベル数×8回起動されるもので回数が多いので、
起動のオーバーヘッドが大きくなるからです。もちろん、この場合には
250 と 260 のループを入れ替えるだけで回数が 1/8 になるのでだいぶましで、
十分に粒子数が多い時にはその程度でよいかもしれません。