高速な転置

# はじめに
本稿は[拓殖大学ディジタルコンテンツ研究愛好会 Advent Calendar 2020](https://adventar.org/calendars/5528) 17日目の記事です。

## 転置とは

> m行n列の行列Aに対してAの(i, j)要素と(j, i)要素を入れ替えてできるn行m列の行列のことである。

[wikipedia 転置行列](https://ja.wikipedia.org/wiki/%E8%BB%A2%E7%BD%AE%E8%A1%8C%E5%88%97)

具体的には、

```math
A = 
\begin{bmatrix}
a & b \\
c & d \\
e & f
\end{bmatrix}
```

としたとき、

```math
A^T = 
\begin{bmatrix}
a & c & e \\
b & d & f 
\end{bmatrix}
```

となる行列を転置行列$A^T$と言います。

## プログラムにすると

```cpp
#include <iostream>
using namespace std;

const int ROW = 3;
const int COLUMN = 2;

int main() {
    int A[ROW][COLUMN] = {
        {1, 2},
        {3, 4},
        {5, 6}
    };
    int A_T[COLUMN][ROW];
    for (int h = 0;h < ROW;h++) {
        for (int w = 0;w < COLUMN;w++) {
            A_T[w][h] = A[h][w];
        }
    }
}
```

こうなります。
これは以下のような処理を行っています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299962/e2ae9d7a-425d-e5f7-a957-ac04840316e2.png)

$O(n^2)$と計算量では最適ですが、今回はキャッシュを活かしてさらに高速化していきます。

## 目指す場所

二重ループ内の

```cpp
A_T[w][h] = A[h][w];
```

の部分です。
キャッシュの特性を用いて高速化します。
上記の例では配列がキャッシュに乗るサイズなので起こりませんが、配列のサイズが大きくなったときキャッシュミスが起きます。
キャッシュミスとはなんなのかについて説明する前に、まずは記憶装置について軽く解説をします。

# 記憶装置

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299962/f2e7b541-6e4f-bae5-91c0-df943910db87.png)

記憶装置は大まかに上図のようになっています。
計算を行うレジスタ、一時的にデータを保持するキャッシュ、データを格納するメインメモリとあります。右に行くほど高速だけど容量が小さく、左に行くほど容量が大きいけれど速度が出ないです。

# キャッシュミス

キャッシュには局所性という、頻繁にアクセスされるものや近くアクセスがある予想されるものを格納する特性があります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299962/9ce8dc79-ca0d-ba6c-6e50-18012561c56f.png)

まずCPUがデータを求めたとき、キャッシュを見に行きます。
そこに欲しいデータが存在すればそれでいいのですが、なかったときはメインメモリを見に行きます。これをキャッシュミスといいます。メインメモリとCPUの速度は大きく違うため、その差だけ時間が掛かるからです。
その時、アクセスしたデータはキャッシュに格納されます。直近でアクセスしたものは近い将来再びアクセスされる可能性が高いからです。これを**時間の局所性**といいます。

また、キャッシュには**空間の局所性**も持ちます。
これは、連続してアクセスされるデータは隣接したアドレス空間に存在している可能性が高いことを利用した性質です。
例えば配列を予想してもらうと分かりやすいかもしれません。前からインデックス順に使用する場合、空間の局所性が活きています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299962/811c37be-2258-66de-86bd-a99ad4204b93.png)

さらに、二次元配列のメモリ配置はRow-major order方式において、以下のようになっています。

|配列|メモリ|
|---|---|
|A[0][0]|0|
|A[0][1]|1|
|A[0][2]|2|
|A[0][3]|3|
|...|...|
|A[0][w]|w|
|A[1][0]|w+1|
|A[1][1]|w+2|
|A[1][2]|w+3|
|...|...|
|A[h][w]|h*w|

参考：[多次元配列のメモリレイアウト方式について](https://qiita.com/tyoshitake/items/83ce47f169c5b037f06c)

ではこれらを踏まえて最初に載せたコードを再び見てみます。

```cpp
#include <iostream>
using namespace std;

int main() {
    int A[ROW][COLUMN];
    int A_T[COLUMN][ROW];
    for (int h = 0;h < ROW;h++) {
        for (int w = 0;w < COLUMN;w++) {
            A_T[w][h] = A[h][w];
        }
    }
}
```

`A[h][w]`へのアクセスは空間の局所性に乗っかっており、キャッシュミスは起きていませんが、`A_T[w][h]`はキャッシュミスが頻発しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299962/d432f38b-4ff7-2d7e-5e6b-9c54ce46aec5.png)

そこでキャッシュの局所性を活かした高速化をしていきます。

# ブロック化

大きい配列を一度に処理しようとするのではなく、幾つかのブロックに分けて処理していくやり方です。
配列の要素がそれぞれ独立しており、動的計画法的な前の処理結果によって値の変わらないモノに有効です。まさしく転置がそれに当てはまるわけです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299962/2935594b-84e8-7d86-fc9d-8ade64517083.png)

今回の転置の場合、キャッシュに乗るサイズにブロック分けをします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299962/3afe143e-a979-fde8-3c3e-19bbad868869.png)

これでブロックの処理はメインメモリにアクセスすることがなくキャッシュ内で完結し、キャッシュミスが大きく減少しました。

## プログラムにする

```cpp
int block_size = 8;
for (int b = 0; b < WIDTH; b += block_size) {
  for (int y = 0; y < HEIGHT ; y++) {
    for (int i = 0; i < block_size; i++) {
      int x = min(b+i, WIDTH-1);
      A_T[x][y] = A[y][x];
      // printf("[%d, %d]\n", x, y);
    }
  }
}
```

コメントアウトされている出力を見てみることで理解が捗るかと思います。

# まとめ

コードが複雑になってしまうので最適化はほどほどにしましょう。

# おわりに

読んで頂きありがとうございました。。

僕は2日目[タイピングゲームを作りながらTypeScriptを始めよう](https://qiita.com/masayasviel/items/d30d19af7ba77ed9a316)にもエントリーしているので読んで頂ければ幸いです。
また、他のエントリもぜひご一読ください。
以上、[拓殖大学ディジタルコンテンツ研究愛好会 Advent Calendar 2020](https://adventar.org/calendars/5528) 17日目でした。

# 参考

- [プログラムを高速化する話](https://www.slideshare.net/KMC_JP/ss-45855264)
- [高性能計算](https://fukushimalab.github.io/hpc_exercise/)