投稿日: 投稿者: guest

MATLABでCUDAコードを動かそう！

はじめに

こんにちは。rej55です。
この記事はrogy Advent Calender 2017の3日目の記事です。

私はロボット技術研究会の中ではMATLABやSimulinkを用いていろいろなものを作る活動をしています。

過去記事

工大祭展示「Raspberry Pi アンプシミュレータ」
MATLAB FBX LoaderとMATLAB 3D Model Rendererの開発
MATLABでGPUコンピューティング

今回の記事について

今回の記事のテーマは「MATLABでCUDAコードを動かそう！」というものです。
内容としては過去記事の「MATLABでGPUコンピューティング」に近いですが、今回の内容はMATLABのGPUコンピューティング機能を直接利用するのではなく、自作のCUDAのプログラムをMATLAB上で利用しようというものになります。

MEXについて

MATLABは、C/C++やFortranといった言語で構築された関数・サブルーチンを呼び出すことができます。
C/C++およびFortranコードをMATLABから呼び出すことができるようにビルドしたファイルをMEXファイルと呼びます。
今回のテーマであるCUDAコードにおいてもこのMEX機能を利用して実現します。

まず、C/C++をベースとしたMEXの利用法について簡単に説明します。
MEXによってC/C++の関数を利用する場合、main関数に代わって次のようなゲートウェイ関数を定義します。

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
// nlhs : 出力の数
// plhs[idx] : idx番目の出力が格納されているメモリのアドレス
// nrhs : 入力の数
// prhs[idx] : idx番目の入力が格納されているメモリのアドレス
/* ここに処理を書く */
}
view raw mex_example.cpp hosted with ❤ by GitHub

MATLABからMEXを呼び出すと、まず上記のmexFunctionが呼び出されます。
mexFunction関数の引数のうちnlhs, nrhsは実際にMATLAB上で実行されたMEX関数の戻り値の個数、引数の個数に対応しています。
*plhs[], *prhs[]は各戻り値および引数が確保されているメモリの先頭アドレスとなります。

MEXにおけるメモリ管理

MEXを実際に使う上で一番難しいのは、メモリ管理だと思っています。
何が難しいのかというと、そのメモリをMATLABが所有権を持つのか、MEXが所有権を持つのかという点がわかりにくいところだと思います。

基本的には、MEXの中で動的に確保したメモリはMEXの中で解放すべきですが、ここで確保したメモリがplhsに渡される場合、すなわち確保したメモリを最終的にMATLAB上で利用する場合はMATLABがそのメモリを解放します。
逆に言えば、plhsに渡すメモリ部分をMEXの内部で解放してしまうと、MATLABが終了時に同じ部分を解放しようとするため、2重Freeが生じる可能性があります。
要するに、MEX実行後にMATLABで使う部分のメモリは解放しないでね、ということです。また、同様の理由でMATLABから渡された部分のメモリもMEX内部で解放してはいけません。

そのほかにMEXを利用するうえではまりやすいポイントとしては、plhsへデータを渡そうとしたとき、適切に行わないとメモリリークが生じる場合があるという点があります。
具体的な例を考えてみましょう。次のようなmexFunction関数を作ったものとします。

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double* ptr;
ptr = mxCalloc(5*5, sizeof(double));
/* ここでptrに数値データを入れていく */
plhs[0] = mxCreateDoubleMatrix(5, 5, mxREAL);
mxSetPr(plhs[0], ptr);
}
view raw mex_incorrect.cpp hosted with ❤ by GitHub

このコード内ではmxSetPr(plhs[0], ptr)でplhs[0]がptrの位置を指すようにポインタの上書きを行っており、一見これで問題ないように見えますが、実際は5*5*sizeof(double)のメモリがリークしています。
なぜならば、mxCreateDoubleMatrixで確保したメモリが置き去りになっているからです。
これを回避するためには、mxCreateDoubleMatrixで先にplhs[0]に対応する部分のメモリを確保し、その部分に対応するポインタを利用することで直接plhs[0]内のデータをいじることが考えられます。
これを実現したのが次の例です。

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double* ptr;
plhs[0] = mxCreateDoubleMatrix(5, 5, mxREAL);
ptr = mxGetPr(plhs[0]);
/* ここでptrに数値データを入れていく */
}
view raw mex_correct.cpp hosted with ❤ by GitHub

このようにすることで潜在的なメモリリークを回避することができます。

MEXでCUDAを使う

ここで本題のCUDAを使う方法ですが、初めに述べたようにMEXを介してCUDAを使うことになります。
GPUのメモリを使う場合は次のような流れでmexFunctionを記述することになります。

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
//変数の定義
const mxGPUArray *src;
mxGPUArray *dst;
const double *d_src;
double *d_dst;
//初期化
mxInitGPU();
//prhs[0]から読み取り専用のmxGPUArrayを生成
src = mxGPUCreateFromMxArray(prhs[0]);
//srcのアドレスを取得
d_src = (const double*)(mxGPUGetDataReadOnly(src));
//GPU上に出力用の領域を確保
dst = mxGPUCreateGPUArray(mxGPUGetNumberOfDimensions(src),
mxGPUGetDimensions(src),
mxGPUGetClassID(src),
mxGPUGetComplexity(src),
MX_GPU_DO_NOT_INITIALIZE);
//dstのアドレスを取得
d_dst = (double *)(mxGPUGetData(dst));
/* カーネル関数の実行など、dstにデータを書き込む */
//dstをMATLABに返すためのmxArrayに変換
plhs[0] = mxGPUCreateMxArrayOnGPU(dst);
//一次的なGPUArrayの解放
mxGPUDestroyGPUArray(src);
}
view raw mexcuda_example.cu hosted with ❤ by GitHub

上記のMEXの使い方がわかっていれば、それと同様に使うだけでOKです。
基本的に注意する点としては、ホスト(CPU)側とデバイス(GPU)側のメモリの扱いに注意することくらいでしょうか。
CUDAで注意すべき点と大体同じだと思います。

例：2つのベクトルの線形和

ここでは2つのベクトルの線形和を計算するMEXを紹介したいと思います。
src1, src2という2つのベクトルをそれぞれk1倍、k2倍して足し合わせるという処理です。
この処理を行うCUDAソースファイルを次のように作成します。

#include "mex.h"
#include "gpu/mxGPUArray.h"
__global__
void vec_add
(
const double* src1,
const double* src2,
const double k1,
const double k2,
double* dst,
const int N
)
{
// Calculate index
const int tid = blockDim.x * blockIdx.x + threadIdx.x;
if(tid < N)
{
dst[tid] = k1*src1[tid] + k2*src2[tid];
}
}
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
// Define variables
const mxGPUArray *src1;
const mxGPUArray *src2;
double k1;
double k2;
mxGPUArray *dst;
const double *d_src1;
const double *d_src2;
double *d_dst;
int N1, N2;
// Check the number of arguments
if ( nrhs != 4 ) {
mexErrMsgIdAndTxt("MATLAB:vec_add","The number of input arguments must be 4.");
}
if ( nlhs != 1 ) {
mexErrMsgIdAndTxt("MATLAB:vec_add","The number of output arguments must be 1.");
}
// Initialization
mxInitGPU();
// Get data from *prhs[]
src1 = mxGPUCreateFromMxArray(prhs[0]);
src2 = mxGPUCreateFromMxArray(prhs[1]);
k1 = mxGetScalar(prhs[2]);
k2 = mxGetScalar(prhs[3]);
// Check the dimension of src vectors
N1 = (int)(mxGPUGetNumberOfElements(src1));
N2 = (int)(mxGPUGetNumberOfElements(src2));
if ( N1 != N2 ) {
mxGPUDestroyGPUArray(src1);
mxGPUDestroyGPUArray(src2);
mexErrMsgIdAndTxt("MATLAB:vec_add","The dimension of input vectors must be same.");
}
// Get address of src1 and src2
d_src1 = (const double*)(mxGPUGetDataReadOnly(src1));
d_src2 = (const double*)(mxGPUGetDataReadOnly(src2));
// Allocate memory of the destination variable on device memory
dst = mxGPUCreateGPUArray(mxGPUGetNumberOfDimensions(src1),
mxGPUGetDimensions(src1),
mxGPUGetClassID(src1),
mxGPUGetComplexity(src1),
MX_GPU_DO_NOT_INITIALIZE);
d_dst = (double *)(mxGPUGetData(dst));
// Call kernel function
dim3 block(N1);
dim3 grid((N1 + block.x - 1) / block.x);
vec_add<<<grid, block>>>(d_src1, d_src2, k1, k2, d_dst, N1);
// Pass dst to plhs[0]
plhs[0] = mxGPUCreateMxArrayOnGPU(dst);
// Release memory
mxGPUDestroyGPUArray(src1);
mxGPUDestroyGPUArray(src2);
}
view raw vec_add.cu hosted with ❤ by GitHub

次に、MATLAB上でこのソースファイルが置いてあるディレクトリ移動し、mexcuda vec_add.cuを実行することでコンパイルが行われます。
すると、vec_addに対応するMEXファイルが生成されるので、これを実行してあげることでvec_add.cuで記述したカーネル関数を呼び出すことができます。

実際に次のテストコードを実行して性能の計測を行ってみました。

%% Initialize
clear
%% Set dimesion
N = 10^6;
%% Set source vectors and scalar values
src1 = ones(N, 1);
src2 = 2*ones(N, 1);
k1 = 2;
k2 = 3;
%% Measure elapsed time with MATLAB
tic;
dst1 = k1*src1 + k2*src2;
Elapsed1 = toc;
%% Measure elapsed time with MEX CUDA
tic;
d_src1 = gpuArray(src1);
d_src2 = gpuArray(src2);
d_dst2 = vec_add(d_src1, d_src2, k1, k2);
dst2 = gather(d_dst2);
Elapsed2 = toc;
%% Display result
disp(['N = ' num2str(N)]);
disp(['Elapsed time with MATLAB : ' num2str(Elapsed1) ' sec.']);
disp(['Elapsed time with MEX CUDA : ' num2str(Elapsed2) ' sec.']);
view raw run_test.m hosted with ❤ by GitHub

今回はせっかくなので新しく使えるようになった東工大のスーパーコンピュータTSUBAME 3.0を利用してみました。
gpuDeviceを実行してGPUデバイスの情報を表示してみると、
>> gpuDevice

ans =

CUDADevice with properties:

Name: ‘Tesla P100-SXM2-16GB’
Index: 1
ComputeCapability: ‘6.0’
SupportsDouble: 1
DriverVersion: 8
ToolkitVersion: 8
MaxThreadsPerBlock: 1024
MaxShmemPerBlock: 49152
MaxThreadBlockSize: [1024 1024 64]
MaxGridSize: [2.1475e+09 65535 65535]
SIMDWidth: 32
TotalMemory: 1.7067e+10
AvailableMemory: 1.6693e+10
MultiprocessorCount: 56
ClockRateKHz: 1480500
ComputeMode: ‘Default’
GPUOverlapsTransfers: 1
KernelExecutionTimeout: 0
CanMapHostMemory: 1
DeviceSupported: 1
DeviceSelected: 1
となりました。現状販売されている最強のGPGPU用GPUと言ってもいい？Tesla P100が搭載されていますね。

run_test.mの実行結果が次のようになります。
>> run_test
N = 1000000
Elapsed time with MATLAB : 0.001871 sec.
Elapsed time with MEX CUDA : 0.006699 sec.
やはり、ベクトルの線形和程度の簡単な演算ではメモリアロケーションなどの余計な処理部分がネックとなってMATLABビルトインの処理より速くなることはないみたいです。
もう少し複雑な関数で勝負すれば勝てるかもしれないですね。

まとめ

本記事ではMATLABで自作CUDAコードを利用するためにMEXを介したCUDAの利用方法について解説しました。
MEXはメモリ管理部分が少し厄介ですが、正しく理解して使えば有効なツールになると思います。
今回紹介した内容はgithub上で公開しております。こちらに公開しておりますので、興味のある方はぜひご利用ください。
この記事が皆様の参考になれば幸いです。

明日はdango_botくんの「Cheeseの活動報告」の記事となります。お楽しみに！

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です