FFTの処理を行う際には、並列化と共にメモリアクセスのパターンを正しく設計しないと遅延の原因となります。

「図：2次元FFTの手順」（図13.2「図：2次元FFTの手順」）では、列スキャンを行いましたが、このままではメモリストライド（ローカルメモリの場合はバンクコンフリクト）が発生します。そのためカラムにアクセスするのではなく、行列の転置をおこない行としてアクセスを行います。

具体的には以下のような手順が必要となります。

転置については、以下のtransposeカーネル関数を使います。

__kernel void transpose(
                         __global float2* input,
                         __global float2* output,
                         size_t width,
                         size_t height)
{

        __local float2 tile[size * (size+1)];
        size_t x = get_global_id(0);
        size_t y = get_global_id(1);

        size_t lx = get_local_id(0);
        size_t ly = get_local_id(1);

        size_t gx = get_group_id(0);
        size_t gy = get_group_id(1);

        size_t index_input = y * width + x;
    size_t index_tile = ly * (size+1) + lx;
    tile[index_tile]=  input[index_input];
        barrier(CLK_LOCAL_MEM_FENCE);

        size_t ox = gy * size + lx;
        size_t oy = gx * size + ly;

        size_t index_output = oy * height + ox;
        index_tile = lx * (size+1) + ly;
    output[index_output] = tile[index_tile];

この転置カーネルはタイル方式を採用しており、sizeパラメータで指定した整数を縦・横の長さにしたタイルをローカルメモリとして作りこみます。

size変数は以下のようにコンパイル時に設定した値が使われます。

program = clCreateProgramWithSource(context, 1, new String[] {sb.toString()}, null, null);
StringBuffer op = new StringBuffer();
op.append("-Werror -Dsize=");
op.append(LOCAL_SIZE);
String option = op.toString();
clBuildProgram(program, 0, null, option, null, null);

ここではLOCAL_SIZEは4として設定していますので、行列を4x4のタイルに分けてから転置を行います。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

ここからタイルを作りますが、まず一番左上の部分行列に注目します。

1

2

3

4

9

10

11

12

17

18

19

20

25

26

27

28

このタイル行列は以下のコードで作ることができます。

__local float2 tile[size * (size+1)];
size_t x = get_global_id(0);
size_t y = get_global_id(1);

size_t lx = get_local_id(0);
size_t ly = get_local_id(1);
size_t index_input = y * width + x;
size_t index_tile = ly * (size+1) + lx;
tile[index_tile]=  input[index_input];

widthとheightは幅と高さを設定した値とし、input配列は行からスキャンしていきます。

size_t index_input = y * width + x;

それをローカルメモリにあるtile配列に再配置します。

tile配列は転置すると以下のようになります。

1

9

17

25

2

10

18

26

3

11

19

27

4

12

20

28

表をみると行と列が入れ替わっていることがわかります。例えば1行目の（1,2,3,4）は、1列目に配置されています。

index_tile = lx * (size+1) + ly;

1次元目のローカルインデックス、lx、2次元目のローカルインデックス、lyは、lx=0,ly=0は「1」、lx=0,ly=1は「2」、lx=3、ly=0は25というように読み込みをします。

これによりindex_tileは、転置する際の順序でtile配列の要素を読み込む際の添字として使用ができます。

読み込んだ要素はグローバルメモリ空間にある出力配列に再配置しますが、この再配置を行うための添字の算出には以下のコードを使います。

(gx * size + lx) * height + gy * size + ly

このコードを考えるために以下の部分配列を例としてみましょう。

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

まず行列の高さは8なので、height＝8となります。次ぎにグループIDはgx=1、gy=0となります。これは一番上（gy=0）にあり、かつ横に一つ（gx=0）タイルをずらした部分（タイル）行列となるからです。

一番左上の要素「5」では、lx=0、ly=0となるので、

(1 * 4 + 0) * 8 + 0 * 4 + 0 = 32

となります。インデックス（添字）が32の場合、配列の値は33となります。（1,9,17,25,33,41,49,57）が1列目なので、その中の左から5番目の要素が33となります。

もう一つ例をみてみましょう。一番左側（gx=0）の一つ下（gy=0）の部分配列を考えてみましょう。

33

34

35

36

41

42

43

44

49

50

51

52

57

58

59

60

先ほど値5のインデックスが値33と交換されることを確認したので、今度は33が5と交換されるか確認をしてみます。

(0 * 4 + 0) * 8 + (1 * 4 + 0) = 4

インデックスが4の配列要素の値は5となるので、これで相互の関係が確認できました。

最後に対角線上(diagonal)の要素が一致するかを確認してみましょう。今回は縦横一つずらし（gx=1,gy=1）のタイル行列を見てみます。

37

38

39

40

45

46

47

48

53

54

55

56

61

62

63

64

例えば37の値（インデックスは36）の配置が変わらないことを確認してみましょう。

(1 * 4 + 0) * 8 + (1 * 4 + 0) = 36

これにより対角線上の要素についても転置対応が確認できました。

上記の例の場合では、ホストコードでは以下のように、設定を行います。

private static final int width = 8;
private static final int height = 8;

private static final int DATA_SIZE = width*height;
private static final int LOCAL_SIZE = 4;

private static long[] global_work_size_rect = new long[]{width,height,1};
private static long[] local_work_size_rect = new long[]{LOCAL_SIZE,LOCAL_SIZE,1};

private static void transpose(cl_mem input_mem, cl_mem output_mem) {

    int[] widthPtr = new int[1];
    widthPtr[0]=width;
    int[] heightPtr = new int[1];
    heightPtr[0]=height;

    clSetKernelArg(kernel_transpose, 0, Sizeof.cl_mem, Pointer.to(input_mem));
    clSetKernelArg(kernel_transpose, 1, Sizeof.cl_mem, Pointer.to(output_mem));
    clSetKernelArg(kernel_transpose, 2, Sizeof.cl_int, Pointer.to(widthPtr));
    clSetKernelArg(kernel_transpose, 3, Sizeof.cl_int, Pointer.to(heightPtr));

    clEnqueueNDRangeKernel(queue,
            kernel_transpose, 2, null,
            global_work_size_rect,
            local_work_size_rect,
            0, null, null);
}

最後にタイルの構成について補記します。以下のコードでは、本来16要素分のローカルメモリ空間があれば十分なのに、size個の要素を余分に割り当てをしています。

__local float2 tile[size * (size+1)];

これは、バンクコンフリクトを避けるためのOpenCLにおけるデータアクセスデザインパターンの一種です。

このようにアクセスが一つずらすと、バンクがsize(コードでは4に設定)の倍数とするのであれば、1つずらすと連続するアクセスが同じバンクをヒットする可能性を低く抑えることが可能となります。

実装例では以下のようなアクセスのパターンとなり、仮に0,4,8,12がバンクだとすると、0,5,10,15とするので、第一バンク(0-0=0)、第二バンク(5-4=1)、第三バンク(10-8=2)、第四バンク(15-12=3)と巧妙に同一バンクでの重複（衝突）がおきないように回避がなされていることがわかるかと思います。

gx

gy 

index_tile

0

0

0

5

10

15

1

6

11

16

index_out

0

0

0

1

2

3

8

9

10

11

index_tile

0

1

0

5

10

15

1

6

11

16

index_out

0

1

4

5

6

7

12

13

14

15

index_tile

1

1

0

5

10

15

1

6

11

16

index_out

1

1

36

37

38

39

44

45

46

47