AIアクセラレータ・IMAXの紹介 ~ (4) SIMD化とループ効率化

# 前回の記事
[AIアクセラレータ・IMAXの紹介 ~ (3) 基本的なデータの扱い方及び演算](https://qiita.com/Prokuma/items/eab7ab1fb4ce49813714)

# 前回のプログラムの問題点
前回は、トーンカーブの例題の解説を行いました。しかし、この例題、**控えめにいってIMAXの本領を全く発揮できていません。** とりあえず、何が問題なのか見ていきましょう。

## 一度に1ピクセルしか処理できない
このプログラムの一番の問題点、**一回の実行で1ピクセルしか処理できません。** 正直データの行き来を考えるとCPUでやった方が電力効率とパフォーマンスが良いかもしれません。

## ループ実行が非効率的
~~前回の記事でIMAXのループの書き方について書いておきましたが、**決まった書き方じゃないとループ実行をCPUに頼ってしまいます。** つまり、前回のプログラムは処理が終わるといったんCPUの制御に戻ってしまうので、効率が非常に悪いのが問題です。~~
修正: `conv-c2c`に関連記述がありました。IMAX内ループ実行にはなっていました。しかし、二重ループにはこの方法では対応できないので、後術する二重ループの書き方は効率向上に役立っていると思います。

## 改善案
一回の実行で複数のピクセルの処理（SIMD化）・IMAXのループ実行機能をうまく使えばIMAXの中でループの制御が可能となり、一回の実行で複数ピクセルの処理をしながら、CPUを介さず全体の変換結果を得られます。「これが問題だ！」と指摘した時点で改善案は「その問題を無くす」に尽きますからね。今回は、そのやり方について例題を交えて説明していきます。

# SIMD化
```c
//EMAX5A begin tone_curve mapdist=0
while (loop--) {
    mop(OP_LDR, 1, &BR[0][1][1], (Ull)(rr++), 0LL, MSK_D0, (Ull)r, 320, 0, 0, (Ull)NULL,320); /* stage#0 */
    mop(OP_LDBR, 1, &BR[1][1][1], (Ull)t1, BR[0][1][1], MSK_B3, (Ull)t1, 64, 0, 0, (Ull)NULL, 64); /* stage#1 */
    mop(OP_LDBR, 1, &BR[1][1][0], (Ull)t1, BR[0][1][1], MSK_B7, (Ull)t1, 64, 0, 0, (Ull)NULL, 64); /* stage#1 */
    mop(OP_LDBR, 1, &BR[1][2][1], (Ull)t2, BR[0][1][1], MSK_B2, (Ull)t2, 64, 0, 0, (Ull)NULL, 64); /* stage#1 */
    mop(OP_LDBR, 1, &BR[1][2][0], (Ull)t2, BR[0][1][1], MSK_B6, (Ull)t2, 64, 0, 0, (Ull)NULL, 64); /* stage#1 */
    mop(OP_LDBR, 1, &BR[1][3][1], (Ull)t3, BR[0][1][1], MSK_B1, (Ull)t3, 64, 0, 0, (Ull)NULL, 64); /* stage#1 */
    mop(OP_LDBR, 1, &BR[1][3][0], (Ull)t3, BR[0][1][1], MSK_B5, (Ull)t3, 64, 0, 0, (Ull)NULL, 64); /* stage#1 */
    exe(OP_CCAT, &r1, BR[1][1][0], EXP_H3210, BR[1][1][1], EXP_H3210, 0LL, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
    exe(OP_CCAT, &r2, BR[1][2][0], EXP_H3210, BR[1][2][1], EXP_H3210, 0LL, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
    exe(OP_CCAT, &r3, BR[1][3][0], EXP_H3210, BR[1][3][1], EXP_H3210, 0LL, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
    exe(OP_MMRG, &r0, r1, EXP_H3210, r2, EXP_H3210, r3, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
    mop(OP_STR, 3, &r0, (Ull)(dd++), 0LL, MSK_D0, (Ull)d, 320, 0, 0, (Ull)NULL,320); /* stage#2 */
}
//EMAX5A end
```
前回に比べて何か増えています。特に、`mop()`でロードを指定する呼び出しが倍になっています。これは、同じアドレス先の64ビットデータから2ピクセルの情報を同時に取得できるためです。一番最初のロードが`OP_LDR`になっていることに注目してください。６４ビットの2ピクセル分のデータを一気にロードし、各部分にマスクをかけ、カラーマップから対応するデータをロードしています。そして`OP_CCAT`でRGBのそれぞれに該当する部分を合体し、`OP_MMRG`でさらに一つのデータにまとめ、最後でストアを行っています。

しかし、`OP_CCAT`自体仕様からなくなったにも関わらず、相変わらず仕様書の例題には存在するので、`OP_CCAT`が使われている部分とその下は次のように変えましょう。
```c
exe(OP_MMRG, &r1, BR[1][1][0], EXP_H3210, BR[1][2][0], EXP_H3210, BR[1][3][0], EXP_H3210, OP_AND, 0xffffffff, OP_SLL, 32LL);
exe(OP_MMRG, &r2, BR[1][1][1], EXP_H3210, BR[1][2][1], EXP_H3210, BR[1][3][1], EXP_H3210, OP_AND, 0xffffffff, OP_NOP, 0LL);
exe(OP_ADD3, &r0, r1, EXP_H3210, r2, EXP_H3210, 0LL, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
mop(OP_STR, 3, &r0, (Ull)(dd++), 0LL, MSK_D0, (Ull)d, 320, 0, 0, (Ull)NULL,320); 
```
それぞれのピクセルデータを前回のプログラム同様併合させて、1ピクセルのデータを作り、その下でこれらのAND演算を行います。EX1の演算ではAND演算を行えないですが、1ピクセル目のデータをシフト演算でずらしているため、足し算で2ピクセルのデータの併合を行えます。一筋縄ではいかないので、少し工夫がいります。

これで理論上のスループットは二倍になりましたが、これでも問題があります。**このループをCPUを介さず動かすことができないという致命的な問題です。** では、そのやり方について見ていきましょう。

# IMAX内ループ実行
前回の記事で特殊な形のループの書き方があると紹介しました。
```c
//EMAX5A begin name mapdist=0
for (CHIP=0;CHIP<NCHIP;CHIP++){
    for (INIT1=1,LOOP1=N,some_var=some_var_init;LOOP1--;INIT1=0) {
        for (INIT0=1,LOOP0=M,some_var_2=some_var_2_init;LOOP0--;INIT0=0) {
        }
    }
}
//EMAX5A end
//EMAX5A drain_dirty_lmm
```
このようなループの書き方ですが、`conv-c2c`で文法扱いになっているので、**IMAX内ループ実行をしたいのであればこの形式を必ず守ってください。**　現状、最大二重ループまで対応しています。このループ自体、特殊文法扱いで、IMAX上で変数の初期化等を行うなどにより、CPUを介さずループを実行できるようになります。

```c
// WD = Width
// RMGRP = Height of A Group of Image

//EMAX5A begin name mapdist=0
for (CHIP=0;CHIP<NCHIP;CHIP++){
    for (INIT1=1,LOOP1=RMGRP,rofs=0-WD*8;LOOP1--;INIT1=0) {
        for (INIT0=1,LOOP0=WD,cofs=0-4;LOOP0--;INIT0=0) {
            exe(OP_ADD, &cofs, INIT0?cofs:cofs, EXP_H3210, 4, EXP_H3210, 0,EXP_H3210,OP_AND,0x0ffffffffLL,OP_NOP,0);
            exe(OP_ADD, &rofs, rofs, EXP_H3210, INIT0?WD*8:0, EXP_H3210, 0,EXP_H3210,OP_NOP,0, OP_NOP, 0);
            exe(OP_ADD, &pofs, rofs, EXP_H3210, cofs, EXP_H3210, 0,EXP_H3210, OP_AND,0x000000ffffffffLL, OP_NOP, 0);

            mop(OP_LDR, 1, &BR[2][1][1], (Ull)rtop0[CHIP], pofs, MSK_D0, (Ull)rtop0[CHIP], WD*RMGRP, 0, 0, (Ull)NULL, WD*RMGRP);
            mop(OP_LDBR, 1, &BR[3][1][1], (Ull)t1, BR[2][1][1], MSK_B3, (Ull)t1, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][1][0], (Ull)t1, BR[2][1][1], MSK_B7, (Ull)t1, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][2][1], (Ull)t2, BR[2][1][1], MSK_B2, (Ull)t2, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][2][0], (Ull)t2, BR[2][1][1], MSK_B6, (Ull)t2, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][3][1], (Ull)t3, BR[2][1][1], MSK_B1, (Ull)t3, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][3][0], (Ull)t3, BR[2][1][1], MSK_B5, (Ull)t3, 64, 0, 0, (Ull)NULL, 64);
            exe(OP_MMRG, &r1, BR[3][1][0], EXP_H3210, BR[3][2][0], EXP_H3210, BR[3][3][0], EXP_H3210, OP_AND, 0xffffffff, OP_SLL, 32LL);
            exe(OP_MMRG, &r2, BR[3][1][1], EXP_H3210, BR[3][2][1], EXP_H3210, BR[3][3][1], EXP_H3210, OP_AND, 0xffffffff, OP_NOP, 0LL);
            exe(OP_ADD3, &r0, r1, EXP_H3210, r2, EXP_H3210, 0LL, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
            mop(OP_STR, 3, &r0, (Ull)dtop0[CHIP], pofs, MSK_D0, (Ull)dtop0[CHIP], WD*RMGRP, 0, 0, (Ull)NULL, WD*RMGRP); 
        }
    }
}
//EMAX5A end
//EMAX5A drain_dirty_lmm
```
先ほどSIMD化したプログラムをループ内に入れてみました。`rofs`は画像の行、`cofs`は画像の列を指定します。`pofs`は、この二つを足し合わせることによって演算対象のピクセルを指定しています。`rofs`は、`INIT0`が１の時にしか加算されません。つまり、内側ループが終了するまで`WD*8`は加算されないのです。ちなみに、ここで`WD`ではなく`WD*8`を足す理由は、一気に処理するデータの量が64ビット（8バイト）単位だからです。
`INIT0?cofs:cofs`は一見何の意味もなしていない三項演算子のように見えますが、これも`conv-c2c`の特殊文法です。「IMAX内で初期化する変数」という指定です。
また、どちらも最初の値を０から引いていますが、これは一番最初のところで足すからです。ゼロからスタートさせるための、一種の工夫です。
このforループは、IMAXの一番最後のユニットにマッピングされ、ループの制御に使用されます。
`CHIP`は、ユニットを線型に繋いだものが複数ある場合の対応です。IMAXではこれらのユニットをつなげたものを「チップ」と呼びます。チップが複数ある場合、複数チップによる並列化が可能となります。
`BR[][][]`レジスタ番号が少し変わっていますが、ループ実行のために最初の段が使われるため、そのためにずらしています。

# ユニットを使い切る
まだまだ高速化の余地はあります。上のようなプログラムだと、効率よくループを回せますが、普通1チップ64ユニットあるIMAXのユニットを使いきれていません。このユニットをできるだけ使い切る、これによって処理できるピクセルの数も増えます。
```c
//EMAX5A begin name mapdist=0
for (CHIP=0;CHIP<NCHIP;CHIP++){
    for (INIT1=1,LOOP1=RMGRP,rofs=0-WD*8;LOOP1--;INIT1=0) {
        for (INIT0=1,LOOP0=WD,cofs=0-4;LOOP0--;INIT0=0) {
            exe(OP_ADD, &cofs, INIT0?cofs:cofs, EXP_H3210, 4, EXP_H3210, 0,EXP_H3210,OP_AND,0x0ffffffffLL,OP_NOP,0);
            exe(OP_ADD, &rofs, rofs, EXP_H3210, INIT0?WD*8:0, EXP_H3210, 0,EXP_H3210,OP_NOP,0, OP_NOP, 0);
            exe(OP_ADD, &pofs, rofs, EXP_H3210, cofs, EXP_H3210, 0,EXP_H3210, OP_AND,0x000000ffffffffLL, OP_NOP, 0);/

            mop(OP_LDR, 1, &BR[2][1][1], (Ull)rtop0[CHIP], pofs, MSK_D0, (Ull)rtop0[CHIP], WD*RMGRP, 0, 0, (Ull)NULL, WD*RMGRP);
            mop(OP_LDBR, 1, &BR[3][1][1], (Ull)t1, BR[2][1][1], MSK_B3, (Ull)t1, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][1][0], (Ull)t1, BR[2][1][1], MSK_B7, (Ull)t1, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][2][1], (Ull)t2, BR[2][1][1], MSK_B2, (Ull)t2, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][2][0], (Ull)t2, BR[2][1][1], MSK_B6, (Ull)t2, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][3][1], (Ull)t3, BR[2][1][1], MSK_B1, (Ull)t3, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[3][3][0], (Ull)t3, BR[2][1][1], MSK_B5, (Ull)t3, 64, 0, 0, (Ull)NULL, 64);
            exe(OP_MMRG, &r1, BR[3][1][0], EXP_H3210, BR[3][2][0], EXP_H3210, BR[3][3][0], EXP_H3210, OP_AND, 0xffffffff, OP_SLL, 32LL);
            exe(OP_MMRG, &r2, BR[3][1][1], EXP_H3210, BR[3][2][1], EXP_H3210, BR[3][3][1], EXP_H3210, OP_AND, 0xffffffff, OP_NOP, 0LL);
            exe(OP_ADD3, &r0, r1, EXP_H3210, r2, EXP_H3210, 0LL, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
            mop(OP_STR, 3, &r0, (Ull)dtop0[CHIP], pofs, MSK_D0, (Ull)dtop0[CHIP], WD*RMGRP, 0, 0, (Ull)NULL, WD*RMGRP); 

            mop(OP_LDR, 1, &BR[5][1][1], (Ull)rtop1[CHIP], pofs, MSK_D0, (Ull)rtop1[CHIP], WD*RMGRP, 0, 0, (Ull)NULL, WD*RMGRP);
            mop(OP_LDBR, 1, &BR[6][1][1], (Ull)t1, BR[5][1][1], MSK_B3, (Ull)t1, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[6][1][0], (Ull)t1, BR[5][1][1], MSK_B7, (Ull)t1, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[6][2][1], (Ull)t2, BR[5][1][1], MSK_B2, (Ull)t2, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[6][2][0], (Ull)t2, BR[5][1][1], MSK_B6, (Ull)t2, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[6][3][1], (Ull)t3, BR[5][1][1], MSK_B1, (Ull)t3, 64, 0, 0, (Ull)NULL, 64);
            mop(OP_LDBR, 1, &BR[6][3][0], (Ull)t3, BR[5][1][1], MSK_B5, (Ull)t3, 64, 0, 0, (Ull)NULL, 64);
            exe(OP_MMRG, &r1, BR[6][1][0], EXP_H3210, BR[6][2][0], EXP_H3210, BR[6][3][0], EXP_H3210, OP_AND, 0xffffffff, OP_SLL, 32LL);
            exe(OP_MMRG, &r2, BR[6][1][1], EXP_H3210, BR[1][2][1], EXP_H3210, BR[6][3][1], EXP_H3210, OP_AND, 0xffffffff, OP_NOP, 0LL);
            exe(OP_ADD3, &r0, r1, EXP_H3210, r2, EXP_H3210, 0LL, EXP_H3210, OP_NOP, 0LL, OP_NOP, 0LL);
            mop(OP_STR, 3, &r0, (Ull)dtop1[CHIP], pofs, MSK_D0, (Ull)dtop1[CHIP], WD*RMGRP, 0, 0, (Ull)NULL, WD*RMGRP); 
            //...
        }
    }
}
//EMAX5A end
//EMAX5A drain_dirty_lmm
```
`BR[63][1][0]`および`BR[63][0][0]`は一般にループ制御のために使われるので、`BR[62][][]`までが実質的な使用限界となります。
**とにかくユニットを使い切る**、が並列性能において重要となります。

# おわりに
今回は、前回のトーンカーブ例題を改善させながら、IMAXで並列性能を向上させる方法について紹介しました。次回は、この基本思想の元、少々特殊な命令を使う例題を紹介していきます。今回はここまで。

# 次回の記事
[AIアクセラレータ・IMAXの紹介 ~ (5) 行列積](https://qiita.com/Prokuma/items/f69a1720f051afb23f2e)

# 参考文献
[EMAX6/ZYNQ64(IMAX2) Architecture Handbook](http://archlab.naist.jp/proj-arm64/doc/emax6/emax6j.pdf)