next up previous
Next: About this document ... Up: 計算天文学 II 第10回 データ構造とアルゴリズム(1) Previous: 5 木構造と再帰的アルゴリズム

6 各関数のテスト

というわけで、一応木構造をつくるためのデータ構造、アルゴリズムとその実 現を紹介した。これを実際に多体問題の数値解を求めるプログ ラムにするためには、ノードの物理データの計算、重力計算、時間積分、初期 条件の生成/読み込み、結果の表示/解析等いろんなことをするプログラムを 作る必要がある。

それらに入る前に、とりあえず今まで出来たところを動かしてみることにする。

割合いろんな関数を作っていろんなことをやらせるので、全体を動かして時に ちゃんと動いているかどうかをどうやって判断し、さらに動いていない時にど こがおかしいかを見つけ出すのはそんなに簡単ではなくなる。

ある程度大きな、沢山の関数からできたプログラムをテスト・デバッグするた めの基本的な考え方は以下の 2 つである。

  1. それぞれの関数を単体でテストする。
  2. 組み合わせたプログラムをテストするが、動作が正しいかどうかを関数ごとにチェックする。

もちろん、プログラムをデバッグする時の大きな問題は、

ということである。特に後者がデバッグを難しいものにする。

というわけで、まず、ツリー構造を作り、それが正しいかどうか調べる。メインプログラムは以下のような感じになる。

pt minus 1 pt

#ifdef TREETEST
int main()
{
    particle * pp;
    int n;
    cerr << "Enter n:";
    cin >> n ;
    pp = new particle[n];
    double rsize = 1.0;
    create_uniform_sphere(pp, n, 0 , rsize);
    for(int i =0;i <n;i++){
        PRC(i); PRL(pp[i].pos);
    }
    bhnode * bn = NULL;

    int nnodes = n*2;
    bn = new bhnode[nnodes];
    bn->assign_root(myvector(0.0), rsize*2, pp, n);
    int heap_remainder = nnodes-1;
    bhnode * btmp = bn + 1;
    bn->create_tree_recursive(btmp,  heap_remainder);
    PRL(bn->sanity_check());
    bn->dumptree();
    return 0;
}
#endif

プログラム全体は

http://grape.astron.s.u-tokyo.ac.jp/%7Emakino/kougi/keisan_tenmongakuII/programs/simpletree/BHtree.C

にあるので細かい関数については省略するが、関数 sanity_check は ちょっと中を見ておこう。

pt minus 1 pt

int bhnode::sanity_check()
{
    int i;
    int iret = 0;
    if (nparticle == 1 ){
        // this is the lowest level node. Things to check:
        // particle is in the cell
        particle * p = pfirst;
        if(inbox(pos,p->pos,l)){
            cerr << "Error, particle out of box ... \n";
            dump();
            return 1;
        }
    }else{

        // This is the non-leaf node. Check the position and side
        // length of the child cells and then check recursively..
        for(i=0;i<8;i++){
            if (child[i] != NULL){
                int err = 0;
                err = child[i]->sanity_check();
                if (l*0.5 != child[i]->l) err += 2;
                myvector relpos = cpos-child[i]->cpos;
                for (int k = 0 ; k<ndim;k++){
                    if (fabs(relpos[k]) !=l*0.25)err += 4;
                }
                if (err){
                    cerr << "Child " << i << " Error type = " << err << endl;
                    dump();
                }
                iret += err;
            }
        }
    }
    return iret;
}
ここでは、以下のようなチェックを行なっている。

  1. ノードに粒子が1つだけなら、それがちゃんと中にあるかどうか。
  2. そうでなければ、自分の子供は正しい位置、大きさを持つかどうか。
で、さらに再帰的に子供についても同じチェックをする。

inbox はこんなものである。

pt minus 1 pt

int inbox(myvector  & cpos, // center of the box
          myvector  & pos,  // position of the particle
          real l)         // length of one side of the box
    
{
    for(int  i = 0; i< ndim; i++){
        if (fabs(pos[i]-cpos[i]) > l*0.5) return 1;
    }
    return 0;
}

ここではベクトルの各要素についてそれぞれ調べる必要がある。

さて、この関数が「エラー」を返して来たとしても、少なくとも以下の4通りの可能性が ある。

というわけで、最初は粒子 1 つとか 2 つで何ができるかみるといったところ から始めないと後が大変ではある。逆に数が少なければいろいろ書き出して人 間がチェックするのも不可能ではない。 以下の関数はツリー構造そのものを再帰的にプリントアウトする。 spc は単に引数の数だけ空白文字を印刷する関数である。

pt minus 1 pt

void bhnode::dumptree(int indent)
{
    int i;
    spc(indent); cerr << "node center pos " << cpos ;
    if (nparticle == 1){
        cerr << " IS LEAF" ;PRL(nparticle);
        particle * p = pfirst;
        for(i = 0; i < nparticle; i++){
            for(int j=0;j<indent+2;j++)cerr << " ";
            PRL(p->pos);
            p= p->next;
        }
    }else{
        cerr << " IS _not_ LEAF ";PRL(nparticle);
        for(i=0;i<8;i++){
            if (child[i] != NULL){
                child[i]->dumptree(indent + 2);
            }
        }
    }
}

以下は $N=4$ での実行結果である

pt minus 1 pt

simpletree>BHtreetest
Enter n:4
nbody = 4,  power_index = 0
i = 0,  pp[i].pos = -0.815405  -0.0255656  0.0535006
i = 1,  pp[i].pos = -0.0911332  -0.533643  0.662584
i = 2,  pp[i].pos = 0.863463  0.136119  0.112189
i = 3,  pp[i].pos = -0.49528  -0.403606  0.751962
bn->sanity_check() = 0
node center pos 0  0  0 IS _not_ LEAF nparticle = 4
  node center pos -0.5  -0.5  0.5 IS _not_ LEAF nparticle = 3
    node center pos -0.75  -0.25  0.25 IS LEAFnparticle = 1
      p->pos = -0.815405  -0.0255656  0.0535006
    node center pos -0.25  -0.75  0.75 IS LEAFnparticle = 1
      p->pos = -0.0911332  -0.533643  0.662584
    node center pos -0.25  -0.25  0.75 IS LEAFnparticle = 1
      p->pos = -0.49528  -0.403606  0.751962
  node center pos 0.5  0.5  0.5 IS LEAFnparticle = 1
    p->pos = 0.863463  0.136119  0.112189

それらしいツリー構造ができていることがわかる。

来週は、力の計算とか、残りのところをプログラムする。



Jun Makino
平成18年1月16日