fract16型の二つのベクトルの対応する要素同士を固定小数点乗算を行います。
fract16は実際にはshort型の別名であるため、C++言語で普通に乗算を行っても固定小数点演算になりません。アセンブリ言語で記述することを思いつくわけですが、ここでBlackfinならではの難しさにぶつかります。できればデュアルMAC機能を使って最高性能で演算を行いたいところですが、そのためには以下の条件を満たさなければなりません。
これらの条件が常に成り立つように制限すると使いにくいルーチンになります。一方でアセンブリ言語でこれらの場合分けを行うと、記述量が増えて面倒です。なにより、私の技量では時間ばかりかかってバグだらけになることが見えています。
そこで、C++言語の拡張機能を使います。VisualDSP++の拡張機能には固定小数点乗算を行う組み込み関数があります。これを使ってどうなるか見てみましょう。
void stk::mul( fract16 dst[], fract16 src1[], fract16 src2[], int count )
{
#pragma no_alias
#pragma different_banks
for ( int i=0; i<count; i++ )
dst[i] = multr_fr1x16( src1[i], src2[i]);
}
上のプログラム自身は単純なベクトル同士の乗算です。ただし、プラグマを使ってコンパイラにヒントを与えているところがポイントです。
no_aliasプラグマは、続くループの中で使われるポインタ同士が互いの別名ではないことを保証します。これはパラメータとして渡す配列同士が重なっていないことに対して、ユーザーが責任を持つことを意味します。具体的にはdstと他の配列の間に重なりがあってはならないことになります。この制限は演算順序を変えるような強い最適化を行うために必要です。
different_banksプラグマはポインタがそれぞれ違うバンクにあってデュアル・アクセスが有効であることを宣言します。コンパイラはこれを受けて強い最適化を行います。実際にはこのプラグマは存在自身が奇妙です。最終的にデータがどこに置かれようが、Blackfinアーキテクチャは正しい実行結果を保証します。ですから、コンパイラはdifferent_banksを指定されなくても黙って最適化すればいいのです。不思議なプラグマです。
さて、このような指定を行った上でVisualDSP++3.5( February Update )が生成するコードを以下に示します。
生成されたコードは非常に優秀です。コンパイラは上に上げた条件をチェックして可能ならばSIMD演算を行い、それができないならSISD演算を行うコードを生成します。このコードは繰り返しが奇数の場合も可能な限りSIMD乗算を行うよう努力します。
また、肝心のループ部分ですが、赤い部分がループボディーです。2行で乗算2回というのは人間が書くコードに比べてまったく遜色ありません。きっちり性能が出ていますし、場所によっては機械ならではの執拗なスケジューリングを行っています。ループのセットアップをしながらそのディレイスロットで条件判定を行い、ループを実行せずに出て行く最適化（[FFA009AE]）など白眉といえるでしょう。
Dumping range 0xffa00970 - 0xffa00a6f by 1 (Assembly)...
[FFA00970] [ -- SP ] = ( R7:7 , P5:5 ) ;
[FFA00972] P2 = R0 ;
[FFA00974] P1 = R2 ;
[FFA00976] P0 = R1 ;
[FFA00978] R7 = [ SP + 0x14 ] ;
[FFA0097A] CC = R7 <= 0 ; // count < 0 なら
[FFA0097C] IF CC JUMP 132 /*0xFFA00A00*/ ; // 終わり
[FFA0097E] R3 = R1 | R2 ;
[FFA00980] R0 = R3 | R0 ;
[FFA00982] R3 = 3 ;
[FFA00984] R0 = R0 & R3 ;
[FFA00986] CC = R0 == 0 ; // 32ビット整列でないなら
[FFA00988] IF ! CC JUMP SISD /*0xFFA009D8*/ ; // SISD処理
[FFA0098A] CC = R7 <= 1 ; // count = 1なら、
[FFA0098C] IF CC JUMP CALC_ONE /*0xFFA009CC*/ ; // 一回だけ実行
[FFA0098E] R0 = 1 ;
[FFA00990] R3 = R7 >>> 1 ;
[FFA00994] R0 = MAX ( R0 , R3 ) ;
[FFA00998] R0 += -1 ;
[FFA0099A] I0 = R1 ;
[FFA0099C] I1 = R2 ;
[FFA0099E] P0 = R0 ;
[FFA009A0] LSETUP ( 16 /*0xFFA009B0*/ , 24 /*0xFFA009B8*/ ) LC0 = P0 ;
[FFA009A4] CC = R0 == 0 ;
[FFA009A6] MNOP || R1 = [ I1 ++ ] || R0 = [ I0 ++ ] ;
[FFA009AE] IF CC JUMP CALC_TWO /*0xFFA009BA*/ ; // count=2ならループは使わない
[FFA009B0] R2.H = R0.H * R1.H , R2.L = R0.L * R1.L || R1 = [ I1 ++ ] || R0 = [ I0 ++ ] ;
[FFA009B8] [ P2 ++ ] = R2 ;
CALC_TWO:
[FFA009BA] R2.H = R0.H * R1.H , R2.L = R0.L * R1.L ;
[FFA009BE] R0 = 1 ;
[FFA009C0] R0 = R7 & R0 ;
[FFA009C2] CC = R0 == 1 ; // 奇数か
[FFA009C4] [ P2 ++ ] = R2 ;
[FFA009C6] P1 = I1 ;
[FFA009C8] P0 = I0 ;
[FFA009CA] IF ! CC JUMP 54 /*0xFFA00A00*/ ( BP ) ; // 偶数なら終わり
CALC_ONE:
[FFA009CC] R0 = W [ P0 ++ ] ( X ) ;
[FFA009CE] R1 = W [ P1 ++ ] ( X ) ;
[FFA009D0] R0.L = R0.L * R1.L ;
[FFA009D4] W [ P2 ++ ] = R0 ;
[FFA009D6] JUMP.S 42 /*0xFFA00A00*/ ;
SISD:
[FFA009D8] R0 = 1 ;
[FFA009DA] R0 = MAX ( R0 , R7 ) || R1 = W [ P1 ++ ] ( X ) || NOP ;
[FFA009E2] P5 = R0 ;
[FFA009E4] I1 = P0 ;
[FFA009E6] R0.L = W [ I1 ++ ] ;
[FFA009E8] I0 = P2 ;
[FFA009EA] LSETUP ( 6 /*0xFFA009F0*/ , 14 /*0xFFA009F8*/ ) LC0 = P5 ;
[FFA009EE] NOP ;
[FFA009F0] R0.L = R0.L * R1.L || R1 = W [ P1 ++ ] ( X ) || NOP ;
[FFA009F8] MNOP || R0.L = W [ I1 ++ ] || W [ I0 ++ ] = R0.L ;
EXIT:
[FFA00A00] ( R7:7 , P5:5 ) = [ SP ++ ] ;
[FFA00A02] RTS ;
この程度のコードならば、高級言語で書いても十分な性能が出ることがわかります。活用しない手はありません。