SSEによる単精度の数値積分の高速化

- SSEによる単精度の数値積分の高速化 -

yosirin-hello / yosi-s@galaxy.ocn.ne.jp

1. はじめに
	数日前, 私は”数値積分をSSEを使って高速化できないだろうか？”と考えました. そこで春休みの時間を利用して, それが実現できるように試行錯誤をしてきました. 最初はSSEを使って積分を４個同時に計算してるのにもかかわらずCで書いたプログラムよりも若干遅い結果しかでませんでした.　しかし, 積分の計算のループ中のプログラムをすべてインラインアセンブリとSSEで書いてみたら, 結果はCだけのものより約４倍高速に積分することが可能になりました. 今回は, きちんと積分できているのかを確かめるためにπ(円周率)を求めて, 実際のπの値と比較することにしました.

2. πを求める関数
	πは上式で求めることができます. なぜこの関数を選んだのかと言うと, [0,1]の間で関数の微分が穏やかで単調減少だからです. もし積分する区間で関数の微分が大きく変化する場合, 積分値に影響が出やすくなります. それを避けるためにこの関数を選びました.

3. Cによる実装

List1. Cによる実装

List1.のポイント

float Calc_PI( int div_count )
{
        float width, x, pi = 0.0f;

        width = 1.0f / ( float )div_count; // (1)

        for( int i = 0; i < div_count; i++ ) {
                x = ( i + 0.5f ) * width; // (2)
                pi += f( x );
        }
        pi *= width; // (3)

        return pi;
}

短形を集めて積分していく方法で実装しました.数値積分の公式で言うと, 中点公式, すなわち0次のニュートン・コーツの公式です. この関数の場合, 1次のニュートン・コーツの公式(台形公式)よりも誤差が半分であることが知られています.

(1) 積分区間 / 分割数で微小間隔を計算しています.
(2) この関数の場合, i + 0.5 をとることで計算の精度を向上させることができます.(Fig1.)
(3) forループの中で毎回widthとf(x)との積をとっていなかったのでここで積をとります.

4. SSEで実装する前に
	まず, List1. をどのように最適化すれば高速化することができるのか考えてみます.　私はFig2.のようにすればCよりも単順に４倍高速化することができると考えました. 次に, SSEで最適化した後に, できるだけ最適化前の４倍高速になるようにするため以下の点を検討してみました.. 1. List1. のforループのiをどうすればよいのか 2. x = ( i + 0.5f ) * widthの部分を並列化する必要がある 3. πを求める関数をSSEで書き直すことができないか 4. πを求める関数に割り算が含まれているので近似的に割り算を実現できないか 1. と2. を検討した結果がFig3. です. XMM5がList1.のforループの i （カウンタ)の働きをしています.　こうすることでカウンタをSSEに最適化することができます. その次の3. はxを並列化したのだから, なんとかしてSSEで式を展開しなければ意味がありません. 幸いにも, 三角関数や指数関数などが含まれていないので簡単にSSEで書くことができます. 最後の4. ですが, RCPPSという近似的に逆数を求めることができる命令があるので高速化のために使うことにしました. 以上の4つのことを考慮してできたのが以下の図です. Fig4. に出てくるSum(x3, x2, x1, x0)はパラメータに与えられた数を合計するというものです. その内部のデータの流れの様子はFig5. ,Fig6. のようにしました.

5. CとインラインアセンブリとSSEによる実装

List2. は関数の部分でRCPPS命令を使っています. List3. は実験のためにRCPPS命令を使わずにDIVPS命令を使って, 積分の精度を向上させています.

List2. CとインラインアセンブリとSSEによる実装(RCPPS版)

List3. CとインラインアセンブリとSSEによる実装(DIVPS版)

#define ALIGNED __declspec(align(16)) // (1)

float Calc_PI_SSE_rcpps( int div_count )
{
        const ALIGNED float _0_5[4] = { 0.5f, 0.5f, 0.5f, 0.5f };
        const ALIGNED float _1_0[4] = { 1.0f, 1.0f, 1.0f, 1.0f };
        const ALIGNED float _4_0[4] = { 4.0f, 4.0f, 4.0f, 4.0f };
        ALIGNED float temp[4]       = { 0.0f, 1.0f, 2.0f, 3.0f };
        ALIGNED float width_arr[4];
        float width, pi;


        width = 1.0f / ( float )div_count;

        width_arr[0] = width_arr[1] = width_arr[2] = width_arr[3] = width;

        __asm {
                movaps xmm5, [temp]
                xorps  xmm7, xmm7  //  pi = 0
                xor    eax, eax    // eax = 0
loop1:
                // x[] = ( tmep[] + _0_5[] ) * width_arr[];
                movaps xmm3, xmm5
                addps  xmm3, [_0_5]
                mulps  xmm3, [width_arr]

                // result[] = 1 / ( 1 + x[]^2 )
                mulps  xmm3, xmm3
                addps  xmm3, [_1_0]
                rcpps  xmm3, xmm3

                // result[0] から result[3]までの合計の計算
                movaps xmm1, xmm3
                shufps xmm1, xmm3, 0B1h
                addps  xmm1, xmm3
                movaps xmm3, xmm1
                shufps xmm1, xmm3, 4Eh
                addps  xmm1, xmm3

                addss  xmm7, xmm1   // pi += result[0]から[3])          
                addps  xmm5, [_4_0] // temp[] += 4

                add    eax, 4
                cmp    eax, div_count
                jl     loop1


                mulss xmm7, [width_arr]
                mulss xmm7, [_4_0]
                movss [pi], xmm7
        }

        return pi;
}

float Calc_PI_SSE_divps( int div_count )
{
        const ALIGNED float _0_5[4] = { 0.5f, 0.5f, 0.5f, 0.5f };
        const ALIGNED float _1_0[4] = { 1.0f, 1.0f, 1.0f, 1.0f };
        const ALIGNED float _4_0[4] = { 4.0f, 4.0f, 4.0f, 4.0f };
        ALIGNED float temp[4]       = { 0.0f, 1.0f, 2.0f, 3.0f };
        ALIGNED float width_arr[4];
        float width, pi;

        width = 1.0f / ( float )div_count;

        width_arr[0] = width_arr[1] = width_arr[2] = width_arr[3] = width;

        __asm {
                prefetcht0  [_1_0]
                movaps xmm5, [temp]
                movaps xmm6, [_1_0]
                xorps  xmm7, xmm7  // pi = 0
                xor    eax, eax    // eax = 0
loop1:
                // x[] = ( tmep[] + _0_5[] ) * width_arr[];
                movaps xmm3, xmm5
                addps  xmm3, [_0_5]
                mulps  xmm3, [width_arr]

                // result[] = 1 / ( 1 + x[]^2 )
                movaps xmm0, [_1_0]
                mulps  xmm3, xmm3
                addps  xmm3, [_1_0]
                divps  xmm0, xmm3

                // result[0] から result[3]までの合計の計算
                movaps xmm1, xmm0
                shufps xmm1, xmm0, 0B1h
                addps  xmm1, xmm0
                movaps xmm0, xmm1
                shufps xmm1, xmm0, 4Eh
                addps  xmm1, xmm0

                addss  xmm7, xmm1 // // pi += result[0]から[3]
                        
                addps  xmm5, [_4_0] // temp[] += 4

                add    eax, 4
                cmp    eax, div_count
                jl     loop1


                mulss xmm7, [width_arr]
                mulss xmm7, [_4_0]
                movss [pi], xmm7
        }

        return pi;
}

(1) メモリーを16byteに境界をするためのマクロです

6. 実験結果
	実験結果から判ったことを以下に示します. ・積分区間の分割数に比例してサイクル数が増加していく (Graph1.) ・計測した積分区間で誤差が最小になったのは, Cの場合は900分割, SSE(RCPPS)の場合は200分割, SSE(DIVPS)の場合は600分割の時だった (Graph1.) ・ RCPPS命令を使うと, やはり誤差が大きくなるが用途によっては十分である (Graph1.) ・コンパイラによる最適化の効果が一番大きかったのはCで, 約3倍高速化した. また, SSEの方は両方とも効果は薄かった(Table1.) ・ RCPPS命令を使うと, 予想以上の約5.4倍の高速化ができた (Table2.) ・ DIVPS命令を使っても約3.8倍の高速化が実現できた (Table2.)

7. 最後に
	SSEの威力はこの数値積分の実験で十分に検証できました. しかし, SSEは単精度(32ビット)の浮動小数点で演算を行うので, FPUの80ビットの浮動小数点にはどうしても精度の面で勝つことができません. また, SSE2の64ビットの精度でも, SSEの浮動小数点よりかはかなり精度が向上しますがそれでもFPUの精度には勝てません. 将来, SIMDで扱える浮動少数点の精度が向上すればFPUよりももっとアドバンテージが高くなるので, それが実現されれば本当にいいですね.

8. 参考文献

1.	C MAGAZINE 2001.10 特集1 プログラムの高速化・最適化, SOFT BANK Publishing
2.	Streaming SIMD Extensions - Matrix Multiplication, Intel Application Note AP-930, 4 Source Code
3.	Intel Pentium4 and Intel Xeon Processor Optimization Reference Manual, Appendex C IA-32 Instruction Latency and Throughput
4.	Using the RDTSC Instruction for Performance Monitoring
5.	Microsoft MSDN Library

9. 改訂履歴

改訂	履歴	日付
0.9	もう少しで完成	02/16/2002
1.0	実験結果を整理して公開	02/16/2002