非同期デリゲートを用いた並列計算

背景

普段の研究では、光学素子の集光シミュレーションにはフレネル・キルヒホッフの回折積分を用いているのですが、諸所の事情から最近計算量が大幅に増え、C#でベタに書いたコードでは実行速度がネックに･･･。
というわけで、とりあえず計算部分をマルチスレッド化させて実行速度を稼いでみました。

実装に関して

C#でマルチスレッドを実現する方法はいくつかありますが、今回は処理を実行するメソッド(積分計算部分)には引数を与えたかったので、以下の記事を参考に非同期デリゲートを用いてコードを書きました。

C# でマルチスレッド非同期デリゲート編: http://sonic64.com/2005-04-12.html

非同期デリゲートは内部でスレッドプールを用いて実現されているようで、スレッドプールにはMSDNに

同時にアクティブにできるスレッドの数が制限されます。既定では、ワーカースレッドが CPU ごとに 25、I/O 完了スレッドが 1,000 に制限されています。

とあるように若干制限がありますが今回の場合はなんら問題ないので、特に何も設定していません。参考までに、これらの制限が問題となる場合は解決策がMSDNに載っています。

マネージスレッドプール: http://msdn.microsoft.com/ja-jp/library/0ka9477y(VS.80).aspx

また、デフォルトではCPUコア数に合わせて自動的にスレッド分割数を設定されるように、下記記事のサンプルコードを参考に、GetCPUCoreCount()って名前でkernel32.dllのGetSystemInfo()関数のラッパーを作成して呼び出しています。

【C#】ＣＰＵのコア数を取得する方法: http://gurizuri0505.halfmoon.jp/20080806/2423

解説みたいなもの

大雑把に今回の計算部分について図(クリックで拡大)で説明すると

という感じになります。ざっくりしすぎていますが…、雰囲気だけでも汲み取っていただければ。

実装コード

何かの参考になれば、、、ということで今回書いたコード(計算部分のみ)を以下に載せておきます。

/*
 * WaveFieldクラスは、位置(x,y)と波動場(u)情報を持ったクラス
 * thisは、伝播元
 * destinationは、伝播先
 * threadNumberは、分割スレッド数
 */
public WaveField PropagetoTo(WaveField destination)
{
    //伝播先の配列数以上には分割しない
    if (threadNumber == 0 || threadNumber > destination.Length)
    {
        //CPU論理コア数(だったっけ･･･)を取得
        threadNumber = SystemInfo.GetCPUCoreCount();
    }

    //分割スレッド数分だけ非同期delegate作成
    FresnelIntegralDelegate[] delegateList
      = new FresnelIntegralDelegate[threadNumber];
    IAsyncResult[] asyncResultList = new IAsyncResult[threadNumber];
    for (int i = 0; i < threadNumber; i++)
    {
        //各スレッドに計算させる伝播先の配列の範囲、を計算
        int start = i * destination.Length / (int)threadNumber;
        int end = (i + 1) * destination.Length / (int)threadNumber;

        FresnelIntegralDelegate integral
         = new FresnelIntegralDelegate(FresnelIntegral);
        //EndInvokeを呼び出すために戻り値を保存(いらないかも？)
        asyncResultList[i] = integral.BeginInvoke(this.x, this.y, this.u,
                    destination.x, destination.y, start, end, null, null);
        delegateList[i] = integral;
    }

    for (int i = 0; i < threadNumber; i++)
    {
        int start = i * destination.Length / (int)threadNumber;
        int end = (i + 1) * destination.Length / (int)threadNumber;
        //各スレッドでの計算結果を受け取り
        Complex[] ux = delegateList[i].EndInvoke(asyncResultList[i]);

        for (int j = start; j < end; j++)
        {
            //伝播先の波動場情報に計算結果をマージ
            destination.u[j] += ux[j - start];
        }
    }

    return destination;
}

/// フレネル積分計算デリゲート
private delegate Complex[] FresnelIntegralDelegate
                           (double[] x0, double[] y0, Complex[] u0,
                             double[] x, double[] y, int start, int end);

/// フレネル積分計算
private Complex[] FresnelIntegral
                   (double[] x0, double[] y0,Complex[] u0,
                    double[] x, double[] y, int start, int end)
{
    /* 与えられた範囲の伝播先の波動場uを計算 */
    return u;
}

結果

伝播元を1000分割・伝播先を5000分割して、研究室・手元にあったPCでテストした結果を載せます(クリックで拡大)。
1000x5000回のフレネル・キルヒホッフ積分の計算にかかった時間を、Enviorment.TickCountで計りました。計測の際には、CnQやEISTによるクロックの変動を除くため本計測の前に余分に計算を行ってCPUを最高速にもって行っています。
上二つはシングルコアになります。デュアルコアでは、どのアーキテクチャでもいい感じに高速化できました。

普段Athlonを使っている自分としては、Core2の場合に最もいい感じに倍速化してしまったのが残念すぎるのですが･･･。この辺、個人的にいろいろ興味深かったので、そのうち各CPUコアのアーキテクチャ(主にL2キャッシュの構造)を踏まえて考察を書こうかなと。

課題

デリゲートを作成せずに、引数ありの匿名メソッド作ってBeginInvokeする方法がわからなかった。まぁいっか。