分散処理で万有引力のシュミレーション

宇宙的なもの作ります。

DemoMPI.zip(プロジェクトファイル一式)

万有引力の法則は
$F=G\frac{Mm}{r^2}$
なので、物体にかかる力は要するに距離の二乗に反比例します。存在する物体それぞれに対してこの式を適用すると計算量は $\frac{n^2}{2}$ になります。nが増えると計算量も指数爆発(?)を起こすので、分散処理しようという話。

MPICH2のインストール
- インストール
- パスを通す
- ローカルでエミュレート
- 複数台で実行
ハローワールド
- Visual C++のインストール
- Platform SDKのインストール
- パスの設定
- 実行
XNAで描画プログラムを作る
- 物体を表すクラス
- メインクラス
MPICH2で計算プログラムを作る
- vector3
- main
検証
まとめ
参考

MPICH2のインストール

分散処理のシステムは色々あるみたいですが今回はMPICH2を使います。ICTスクールではMPI.NETという.NETでのMPI実装を使ったので出来ればそっちが良かったのですが、複数台のマシンで動かすにはWindowsServerの何とかエディションみたいなのが必要なんでこっちにしました。ただ、その何とかエディション自体は180日評価版が落としてこれ、開発だけなら普通のXPでも出来ます。他にもmpiJavaとかMapReduce系もあるみたいです。

インストール

まずは、公式サイトからインストーラーをダウンロードしてきます。

実行するとゴチャゴチャと聞いてきますが全部デフォルトで大丈夫です。また、インストールにはMicrosoft Visual C++ 2005 SP1 再頒布可能パッケージ (x86)が必要なのでなければこちらもインストールしておきます。

パスを通す

MPIの実行ファイルにパスを通します。マイコンピュータを右クリック -> [プロパティ] -> [詳細設定] ->システム環境変数のPATHを編集。末尾に;C:\Program Files\MPICH2\binを追加。コマンドプロンプトからmpiexecを実行してUsageが出るのを確認します。
(※)MPIはどの実装でも大体mpiexecが存在するようです。なので、MPICH2とMPI.NETを共存させたりするときは名前の衝突に気をつけましょう。

ローカルでエミュレート

次のコマンドを実行します

mpiexec -localonly 4 C:\Program Files\MPICH2\examples\cpi.exe

Enter the number of intervals:と聞かれるので適当に数字を入れると円周率の計算が始まり、かかった時間と誤差が表示されます。-localonly(4はCPUの数)で分散処理をエミュレートできます。これでインストールは完了です。

複数台で実行

実際に複数のマシンを繋いで実行してみます。今回はMacBook(Intel(R) Core(TM)2 CPU T7400 @ 2.16GHz (2 CPUs))とLetsNoteR5(Genuine Intel(R) CPU U1300 @ 1.06GHz)の2台を使用しました。

まず、全てのマシンにメインのマシンと同じユーザー名とパスワードで管理者権限のアカウントを作ります。次に、MPICH2をインストールし、パスを通します。マシンを同じネットワーク上に接続し、メインマシンで[すべてのプログラム] -> [MPICH2] -> [wmpiconfig.exe]を起動します。Domainをネットワークのドメインにし、Get Hostsを押すとノード一覧が表示されます。ApplyAllを押し、コマンドでmpiexec -registerと実行すると準備完了です。

全てのマシンの同じディレクトリ(ここではC:\tmpとします)に実行したいプログラムを置き、メインマシンでコマンドを実行します。

mpiexec -hosts 2 macbook letsnoter5 C:\tmp\cpi.exe

-hostsはマシンの数で、macbookとletsnoter5はコンピューター名です。各マシンのCPU使用率を眺めていると分散されている様子が良くわかります。

ハローワールド

MPICH2を実行する環境が出来たので、実際にプログラムを開発する環境も設定します。

Visual C++のインストール

VisualC++はExpressEditionが無償で入手できます。ここでは2005で説明しますが、最新版の2008でも多分あんまり変わらないと思います。
公式サイトからインストーラーをダウンロードしてインストール。再びゴチャゴチャ聞かれますが気にせずデフォルトで。

Platform SDKのインストール

今回はWindowsのネットワーク周りのAPIを使うのでPlatformSDKもインストールします。手順はhttp://www.microsoft.com/japan/msdn/vstudio/express/2005/visualc/usingpsdk/に詳しく書いてあるのでここでは割愛します。

パスの設定

Visual C++にMPIライブラリのパスを設定します。[ツール] -> [オプション] -> [プロジェクトおよびソリューション] -> [VC++ディレクトリ]のインクルードファイルにC:\Program Files\MPICH2\includeを、ライブラリファイルにC:\Program Files\MPICH2\libをそれぞれ追加します。

適当な空のコンソールアプリケーション(HelloMPIとします)を新規作成し、[プロジェクト] -> [HelloMPIのプロパティ] -> [構成プロパティ] -> [リンカ] -> [入力] -> [追加の依存ファイル]にmpi.lib cxxd.lib(Releaseの場合はcxx.lib)を追加します。

実行

次のソースファイルを追加します。

#include <iostream>
#include "mpi.h"

int main(int argc,char **argv){
    MPI::Init(argc,argv);
    int size = MPI::COMM_WORLD.Get_size();
    int rank = MPI::COMM_WORLD.Get_rank();
    std::cout << "Hello MPI! rank:" << rank << " size:" << size << std::endl;
    MPI::Finalize();
    std::cin.get();
}

MPIの開始にはInitに実行引数を渡し、終了はFinalizeを呼び出します。Get_sizeとGet_rankはそれぞれノード数と自身のノード番号を表します。実行すると"Hello MPI! rank:0 size:1"と表示されます。

XNAで描画プログラムを作る

環境が設定できたので実際のプログラムを作り始めます。シュミレートの様子を描画するプログラムと、実際にシュミレートを行うプログラムを作り、それらを通信させます。便宜上これらをXNA側、MPI側と呼ぶことにします。まずはXNA側を作ります。Processingにしようか迷いましたが、XNAのが速い(多分)のと大体出来てるのがあるというのでこっちにしました。MPI側の方をmpiJavaにするとProcessingだといいことが起きそうな気がしなくも無いです。XNAの開発環境は構築済みとします。

物体を表すクラス

XNA側では物体の座標情報のみを保持し、毎ステップMPI側から新しい座標を受信します。送られてくる情報を描画するだけなのでこちらは非常に単純です。

public class Cube
{
    public Cube(Model model)
    {
        this._position = Vector3.Zero;
        this._model = model;
        this._transforms = new Matrix[this._model.Bones.Count];
        this._model.CopyAbsoluteBoneTransformsTo(this._transforms);
    }

    public void Update(float x, float y, float z)
    {
        this._position.X = x;
        this._position.Y = y;
        this._position.Z = z;
    }

    public virtual void Draw()
    {
        foreach (ModelMesh mesh in this._model.Meshes)
        {
            foreach (BasicEffect effect in mesh.Effects)
            {
                effect.EnableDefaultLighting();
                effect.World = this._transforms[mesh.ParentBone.Index] * Matrix.CreateTranslation(this._position);
                effect.View = Camera.Instance.View;
                effect.Projection = Camera.Instance.Projection;
            }
            mesh.Draw();
        }
    }

    private Vector3 _position;
    private Model _model;
    private Matrix[] _transforms;
}

CubeはXNA側に表示する物体を表します。メインループのUpdate内でMPI側から所得した座標を設定します。

メインクラス

メインクラスはCubeのリストと通信関連のフィールドを保持します

private List<Cube> _nodes;
private int N = 10;
private TcpClient _tcp;
private NetworkStream _ns;
private SpriteFont _font;

Nはノードの数です。初期化時にMPI側と接続を確立し、ノード数を所得します。

public Game1()
{
    graphics = new GraphicsDeviceManager(this);
    Content.RootDirectory = "Content";
    this.IsMouseVisible = true;
    this._nodes = new List<Cube>();
    this._tcp = new TcpClient("localhost", 12345);
    this._ns = this._tcp.GetStream();
    byte[] buff = new byte[4];
    this._ns.Read(buff, 0, 4);
    this.N = BitConverter.ToInt32(buff, 0);
}

LoadContentでノードを初期化します。デバッグ用のフォントも用意しておきます。

protected override void LoadContent()
{
    spriteBatch = new SpriteBatch(GraphicsDevice);
    this._font = Content.Load<SpriteFont>("SpriteFont1");
    for (int i = 0; i < N; i++)
    {
        this._nodes.Add(new Cube(Content.Load<Model>("blue")));
    }
}

protected override void UnloadContent()
{
    base.UnloadContent();
    this._ns.Close();
    this._tcp.Close();
}

Unloadの方ではネットワーク周りのリソースを開放します。
Drawではノードの描画とFPSの表示を行います。

protected override void Draw(GameTime gameTime)
{
    graphics.GraphicsDevice.Clear(Color.CornflowerBlue);
    graphics.GraphicsDevice.RenderState.DepthBufferEnable = true;
    foreach (Cube n in this._nodes) n.Draw();
    spriteBatch.Begin();
    spriteBatch.DrawString(this._font, "" + 1000.0 / (0.0000001 + gameTime.ElapsedGameTime.Milliseconds), Vector2.Zero, Color.White);
    spriteBatch.End();

    base.Draw(gameTime);
}

最後に、Updateでネットワーク越しに座標を所得し、キューブの位置を更新します。

protected override void Update(GameTime gameTime)
{
    InputManager.PreUpdate();
    Camera.Instance.Update();
    for (int i = 0; i < N; i++)
    {
        byte[] buff = new byte[4];
        this._ns.Read(buff, 0, 4);
        float x = BitConverter.ToSingle(buff, 0);
        this._ns.Read(buff, 0, 4);
        float y = BitConverter.ToSingle(buff, 0);
        this._ns.Read(buff, 0, 4);
        float z = BitConverter.ToSingle(buff, 0);
        this._nodes[i].Update(x, y, z);
    }
    if (InputManager.CurrentKey.IsKeyDown(Keys.Escape))
    {
        this._ns.WriteByte(1);
        this.Exit();
    }
    else
    {
        this._ns.WriteByte(0);
    }
    InputManager.PostUpdate();
    base.Update(gameTime);
}

InputManagerとCameraはそれぞれユーティリティクラスです。また、キューブの更新が終了した後MPI側に応答メッセージを送信します。Escキーが押されていた場合、終了メッセージを送信しプログラムを終了します。

XNA側は以上です。次にMPI側を作ります。

MPICH2で計算プログラムを作る

こちらのプログラムでは各ノードが座標を共有し、更新処理を分割して行います。各ノードが計算した結果をヘッドノードに集めてXNA側に送信します。

vector3

まずは3次元の座標を表す構造体と必要な演算用の関数を定義します。

struct Vector3{
  float x;
  float y;
  float z;
};

Vector3 subVector(Vector3 a, Vector3 b);
Vector3 multVector(Vector3 a, float b);
float squareMagnitude(Vector3 a);
Vector3 Normalize(Vector3 a);

MPIでは基本的にプリミティブか構造体しか送受信できないので、クラスではなく構造体を使っています。vector3.hというヘッダファイルを作りこれらを定義し、vector3.cppに実装部を書きます。

main

まずは必要な変数を用意します

int main(int argc,char **argv)
{
    int size, rank, N=100;
    if(argc == 2) N = atoi(argv[1]);
    Vector3 *positions = new Vector3[N];
    Vector3 *velocities = new Vector3[N];

    // for network
    WSADATA wsaData;
    SOCKET sock0;
    struct sockaddr_in addr;
    struct sockaddr_in client;
    int len;
    SOCKET sock;

positionsは座標を表し、velocitiesは速度を表します。
次に計算するノードのペアの配列を用意します。

    // init tasks
    int jobLen = 0;
    for (int i = 0; i < N; i++) jobLen += i;
    std::pair<int, int> *jobs = new std::pair<int, int>[jobLen];
    int index=0;
    for (int i = 0; i < N; i++){
        for (int j = i + 1; j < N; j++){
            jobs[index].first = i;
            jobs[index].second = j;
            index++;
        }
    }

    // mpi init
    MPI::Init(argc,argv);
    size = MPI::COMM_WORLD.Get_size();
    rank = MPI::COMM_WORLD.Get_rank();
    int start = (jobLen / size) * rank + (std::min)(rank, jobLen % size);
    int end = (jobLen / size) * (rank + 1) - 1 + (std::min)(rank + 1, jobLen % size);

jobsは全ての組み合わせのリストで、各ノードはjobsのstartからendまでのペアについての計算を行います。
次にvector3のデータタイプを定義します

    // create datatype
    int counts[] = {1, 1, 1};
    MPI::Aint byte[3] = {0, 4, 8};
    MPI::Datatype mpiVec3, types[] = {MPI::FLOAT, MPI::FLOAT, MPI::FLOAT};
    mpiVec3 = MPI::Datatype::Create_struct(3, counts, byte, types);
    mpiVec3.Commit();

MPI上で構造体を送受信するためにはその構造体のDatatypeを定義する必要があります。構造体からDatatypeを作成するにはMPI::Datatype::Create_structを使います。引数は順に構造体内のデータの種類、それぞれのデータの次元、開始バイト、MPIでのDatatypeです。作成したDatatypeはCommitして登録します。
次に、キューブの座標とネットワークを初期化します。

    if (rank == 0){
        // init positions
        int max = 50;
        for (int i = 0; i < N; i++){
            positions[i].x = (float)(rand() % max - max / 2);
            positions[i].y = (float)(rand() % max - max / 2);
            positions[i].z = (float)(rand() % max - max / 2);

            velocities[i].x = (float)rand()/RAND_MAX - 0.5f;
            velocities[i].y = (float)rand()/RAND_MAX - 0.5f;
            velocities[i].z = (float)rand()/RAND_MAX - 0.5f;
        }

        // init network
        if (WSAStartup(2, &wsaData) != 0) { return 1; }
        sock0 = socket(AF_INET, SOCK_STREAM, 0);
        addr.sin_family = AF_INET;
        addr.sin_port = htons(12345);
        addr.sin_addr.S_un.S_addr = INADDR_ANY;
        bind(sock0, (struct sockaddr *)&addr, sizeof(addr));
        listen(sock0, 5);
        len = sizeof(client);
        sock = accept(sock0, (struct sockaddr *)&client, &len);
        send(sock, (char*)&N, 4, 0);
    }

初期化はそれぞれ代表してヘッドノードが行います。接続が確立されるとノード数を送信します。ちなみにこの間、他のノードはメインループ先頭の座標を同期する部分で待機しています。
次にメインループです。

    clock_t prevTime = clock();
    while (true){
        // sync position
        MPI::COMM_WORLD.Bcast(positions, N, mpiVec3, 0);

        // calc accele
        float *acc = new float[N * 3];
        for(int i=0; i<N*3; i++) acc[i] = 0;
        for (int i = start; i <= end; i++){
            Vector3 sub = subVector(positions[jobs[i].second], positions[jobs[i].first]);
            float f = 1 / squareMagnitude(sub);
            Vector3 F = multVector(Normalize(sub), f);
            acc[jobs[i].first * 3] += F.x;
            acc[jobs[i].first * 3 + 1] += F.y;
            acc[jobs[i].first * 3 + 2] += F.z;
            acc[jobs[i].second * 3] -= F.x;
            acc[jobs[i].second * 3 + 1] -= F.y;
            acc[jobs[i].second * 3 + 2] -= F.z;
        }

ループの最初にBcastでヘッドノードから各ノードに座標を送信します。引数は順に同期するオブジェクト、次元数、データタイプ、送信元ランクです。座標を同期したらそれを元に加速度を計算します。各ノードは担当する範囲の計算をし、結果をfloatの配列に格納します。ここで、float配列はキューブ数の3倍の大きさで、順にx, y, z座標を格納しています。
次に計算結果をヘッドノードに集めます。

        // reduce accele
        float *newAcc = new float[N * 3];
        MPI::COMM_WORLD.Reduce(acc, newAcc, N * 3, MPI::FLOAT, MPI::SUM, 0);

Reduceは一箇所にデータを集めつつ、それらを演算するためのメソッドです。引数は送信するデータ、収集結果を格納するデータ、次元数、データタイプ、演算子、集めるノードです。ここでは全ノードがnewAcc += accを行っている様なイメージです。
次にデータを送信します。

        char retval;
        if (rank == 0){
            clock_t currentTime = clock();
            float elapsed = (float)(currentTime - prevTime)/CLOCKS_PER_SEC;
            for (int i = 0; i < N; i++){
                velocities[i].x += newAcc[i * 3] * elapsed;
                velocities[i].y += newAcc[i * 3 + 1] * elapsed;
                velocities[i].z += newAcc[i * 3 + 2] * elapsed;
                positions[i].x += velocities[i].x * elapsed;
                positions[i].y += velocities[i].y * elapsed;
                positions[i].z += velocities[i].z * elapsed;

                send(sock, (char*)&positions[i].x, 4, 0);
                send(sock, (char*)&positions[i].y, 4, 0);
                send(sock, (char*)&positions[i].z, 4, 0);
            }

            recv(sock, &retval, 1, 0);
            prevTime = currentTime;
        }

速度に加速度を加算し、位置に速度を加算します。それらのデータを送信し、応答メッセージを受け取ります。

        MPI::COMM_WORLD.Bcast(&retval, 1, MPI::CHAR, 0);
        if(retval == 1) break;
    }

    if(rank == 0){
        closesocket(sock);
        WSACleanup();
    }

    MPI::Finalize();
}

最後に応答メッセージを同期し、終了判定を行います。MPI側も完成です。

検証

ではプログラムを動かしてみます。MPI側、XNA側の順に立ち上げます。

デフォルトでキューブの数は100個ですが、MPI側の実行引数で数を変えれます。僕の環境では200個程度までなら30FPSはキープできました。複数マシンではどうなるでしょう。同じように200個でやってみます。

ほとんど変わりません。試しに円周率計算プログラムも比較してみると

C:\Documents and Settings\NaoyukiTotani>mpiexec -localonly 2 c:\tmp\cpi
Enter the number of intervals: (0 quits) 1000000000
pi is approximately 3.1415926535900072, Error is 0.0000000000002141
wall clock time = 7.675634
Enter the number of intervals: (0 quits) ^D
No number entered; quitting

C:\Documents and Settings\NaoyukiTotani>mpiexec -hosts 2 macbook letsnoter5 c:\tmp\cpi
Enter the number of intervals: (0 quits) 1000000000
pi is approximately 3.1415926535900072, Error is 0.0000000000002141
wall clock time = 15.524865
Enter the number of intervals: (0 quits) ^D
No number entered; quitting

2台だと倍近くの時間がかかっています。並列化して遅くなる原因はプログラムの書き方やネットワークの設定等、様々な要素があります。マシンに性能差がある場合単純に半分づつにするとこうなるようです。

まとめ

C++難しい。