Help us understand the problem. What is going on with this article?

arrayfunを利用したステンシル計算

More than 3 years have passed since last update.

はじめに

この記事はarrayfunを利用したステンシル計算の方法のメモになります.
バージョンはR2015bを使用しています.
下記のGPUを利用する場合はParallel Computing Toolboxが必要です.

arrayfunとは

[B1, B2, ..., Bm] = arrayfun(fun, A1, A2, ..., An)はMATLABで提供されている関数の一つで, n個の引数からm個の値を返す関数funを, A1, A2, ... , Anの各要素に適用してB1, B2, ..., Bmに返す関数です.
反復法的に書くと,

for i = 1:n
    [B1(i), B2(i), ..., Bm(i)] = fun(A1(i), A2(i), ..., An(i))
end

という処理をまとめてやってくれるものになります.
この反復法的な記法で言うところのiに対してすべて独立な処理を行う場合はarrayfunで簡単に実装できるのですが, 配列の中心からその近傍を見るステンシル計算をやる場合には若干面倒になります.

実装例

今回は入力引数に行列のサイズを与え, 全要素が1のN x Nサイズの行列を作り, 配列Aの各要素に対して上:A(i-1, j), 下:A(i+1, j), 左:A(i, j-1), 右:A(i, j+1)の値と自分自身A(i, j)の値の値を足すという処理を行います.
また, 領域外は0となるようにします.

ステンシル計算を行う場合は, arrayfunの引数である関数funに与える引数として, 行番号と列番号を与えるようにします.
この場合, 対象となる配列Aをfunが読み込めるような状況にする必要があります.
この解決策としては, グローバル変数を使う, 入れ子関数を利用するなどが挙げられます.

理由は後述しますが, 今回は入れ子関数を利用して実現します.
そのサンプルコードが以下のとおりです.

stencil.m
function X = stencil(N)
%% 初期化と準備
A = ones(N);
rows = [1:N];
cols = [1:N];
[cols, rows] = meshgrid(cols, rows);

%% ステンシル計算の中身の定義
    function out = mysum(row, col)
        % 上下左右の行番号と列番号
        up = row-1;
        down = row+1;
        left = col-1;
        right = col+1;

        % 上下の境界処理
        if(row == 1)
            rowsum = A(down, col);
        elseif(row==N)
            rowsum = A(up, col);
        else
            rowsum = A(up, col) + A(down, col);
        end

        % 左右の境界処理
        if(col == 1)
            colsum = A(row, right);
        elseif(col==N)
            colsum = A(row, left);
        else
            colsum = A(row, left) + A(row, right);
        end

        out = A(row, col) + rowsum + colsum;
    end

%% arrayfunの実行
X = arrayfun(@mysum, rows, cols);
end

GPUでの処理

GPUで同様の処理を行う場合は, 通常の行列の代わりにgpuArrayを利用すればよいです.
ただし, GPUでarrayfun利用する場合にはグローバル変数機能を利用できません.
これが入れ子関数を利用する理由となります.

また, GPUでarrayfunを利用する場合, gpuArrayのメソッドの一つであるcolonを利用して行番号と列番号を並べたベクトルを作ってそれをそのまま引数として与えると, 自動的にmeshgridを行ってくれる機能があるようです.

上記処理をGPUでも対応できるようにしたバージョンがこちらになります.

mystencil_GPU.m
function stencil_GPU(Usegpu, N)
%% 初期化と準備
if(Usegpu==1)
    gpu = gpuDevice;
    A = gpuArray(ones(N));
    rows = gpuArray.colon(1, N)';
    cols = gpuArray.colon(1, N);
else
    A = ones(N);
    rows = [1:N];
    cols = [1:N];
    [cols, rows] = meshgrid(cols, rows);
end

%% ステンシル計算の中身の定義
    function out = mysum(row, col)
        up = row-1;
        down = row+1;
        left = col-1;
        right = col+1;
        if(row == 1)
            rowsum = A(down, col);
        elseif(row==N)
            rowsum = A(up, col);
        else
            rowsum = A(up, col) + A(down, col);
        end

        if(col == 1)
            colsum = A(row, right);
        elseif(col == N)
            colsum = A(row, left);
        else
            colsum = A(row, left) + A(row, right);
        end                
        out = A(row, col) + rowsum + colsum;
    end

%% arrayfunの実行
A = arrayfun(@mysum, rows, cols);

%% GPU利用の場合はgatherでローカルワークスペースに移す
if(Usegpu==1)
    wait(gpu);
    X = gather(A);
else
    X = A;
end

end

CPUとGPUでの処理速度の比較

CPUとGPUでの処理速度を比較してみたいと思います.
私のGPUの環境は以下のとおりです.

>> gpuDevice

ans = 

  CUDADevice のプロパティ:

                      Name: 'GeForce GTX 760 Ti OEM'
                     Index: 1
         ComputeCapability: '3.0'
            SupportsDouble: 1
             DriverVersion: 7.5000
            ToolkitVersion: 7
        MaxThreadsPerBlock: 1024
          MaxShmemPerBlock: 49152
        MaxThreadBlockSize: [1024 1024 64]
               MaxGridSize: [2.1475e+09 65535 65535]
                 SIMDWidth: 32
               TotalMemory: 2.1475e+09
           AvailableMemory: 1.7663e+09
       MultiprocessorCount: 7
              ClockRateKHz: 980000
               ComputeMode: 'Default'
      GPUOverlapsTransfers: 1
    KernelExecutionTimeout: 1
          CanMapHostMemory: 1
           DeviceSupported: 1
            DeviceSelected: 1

この状況で, 時間計測のため先ほどのstencil_GPUを以下のように変更します.

stencil_GPU_time.m
function X = stencil_GPU_time(Usegpu, N)
%% 初期化と準備
if(Usegpu==1)
    gpu = gpuDevice;
    A = gpuArray(ones(N));
    rows = gpuArray.colon(1, N)';
    cols = gpuArray.colon(1, N);
else
    A = ones(N);
    rows = [1:N];
    cols = [1:N];
    [cols, rows] = meshgrid(cols, rows);
end
%% ステンシル計算の中身の定義
    function out = mysum(row, col)
        up = row-1;
        down = row+1;
        left = col-1;
        right = col+1;
        if(row == 1)
            rowsum = A(down, col);
        elseif(row==N)
            rowsum = A(up, col);
        else
            rowsum = A(up, col) + A(down, col);
        end

        if(col == 1)
            colsum = A(row, right);
        elseif(col == N)
            colsum = A(row, left);
        else
            colsum = A(row, left) + A(row, right);
        end                
        out = A(row, col) + rowsum + colsum;
    end

%% 計測開始
timer = tic();
A = arrayfun(@mysum, rows, cols);
if(Usegpu==1)
    wait(gpu);
    X = gather(A);
else
    X = A;
end
%% 計測終了
toc(timer)

end

それでは実行してみましょう.

>> stencil_GPU_time(0, 100);
経過時間は 0.030878 秒です。
>> stencil_GPU_time(1, 100);
経過時間は 0.168936 秒です。
>>
>> stencil_GPU_time(0, 1000);
経過時間は 2.847233 秒です。
>> stencil_GPU_time(1, 1000);
経過時間は 0.009505 秒です。
>>
>> stencil_GPU_time(0, 10000);
経過時間は 287.057677 秒です。
>> stencil_GPU_time(1, 10000);
経過時間は 0.849652 秒です。

以上のように, 要素数が少ない場合は高速化されませんでしたが, 要素数が非常に多くなると圧倒的にスピードが違うことがわかります.

まとめ

arrayfunでステンシル計算をさせる場合には, 適用させる関数が配列にアクセスできるようにするために, 入れ子関数などで適用させる関数を定義する必要があります.
GPUを利用する場合は配列の代わりにgpuArrayを利用すれば簡単に使えます.

参考

MATLABドキュメンテーション : arrayfun
http://jp.mathworks.com/help/matlab/ref/arrayfun.html
MATLABドキュメンテーション : GPUでのステンシル演算
http://jp.mathworks.com/help/distcomp/examples/stencil-operations-on-a-gpu_ja_JP.html

g_rej55
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした