ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

GPU コンピューティングへの 3 つのアプローチの説明: マンデルブロ集合

この例では、単純な、よく知られた数学上の問題であるマンデルブロ集合を、MATLAB® コードでどのように記述できるかを説明します。次に Parallel Computing Toolbox™ を使用して、GPU ハードウェアを利用するように 3 つの方法でこのコードを適合させます。

  1. GPU データを入力として、既存のアルゴリズムを使用

  2. arrayfun を使用して各要素で個別にアルゴリズムを実行

  3. MATLAB/CUDA インターフェイスを使用して既存の CUDA/C++ コードを実行

設定

下記の値は、大きなカージオイドとその左にある p/q バルブ間の谷にあたる、マンデルブロ集合を大きくズームした部分を指定します。

この範囲内に実数部 (X) と虚数部 (Y) の 1000x1000 のグリッドが作成され、マンデルブロ アルゴリズムがそれぞれのグリッドの位置で繰り返されます。この特定位置では、画像を完全にレンダリングするには 500 回の反復で十分です。

maxIterations = 500;
gridSize = 1000;
xlim = [-0.748766713922161, -0.748766707771757];
ylim = [ 0.123640844894862,  0.123640851045266];

MATLAB でのマンデルブロ集合

以下に示すのは、CPU で実行される標準 MATLAB コマンドを使用したマンデルブロ集合の実装です。これは、Cleve Moler の電子書籍『Experiments with MATLAB』に記載されているコードに基づいています。

この計算は、すべての位置が同時に更新されるようにベクトル化されています。

% Setup
t = tic();
x = linspace( xlim(1), xlim(2), gridSize );
y = linspace( ylim(1), ylim(2), gridSize );
[xGrid,yGrid] = meshgrid( x, y );
z0 = xGrid + 1i*yGrid;
count = ones( size(z0) );

% Calculate
z = z0;
for n = 0:maxIterations
    z = z.*z + z0;
    inside = abs( z )<=2;
    count = count + inside;
end
count = log( count );

% Show
cpuTime = toc( t );
fig = gcf;
fig.Position = [200 200 600 600];
imagesc( x, y, count );
colormap( [jet();flipud( jet() );0 0 0] );
axis off
title( sprintf( '%1.2fsecs (without GPU)', cpuTime ) );

gpuArray の使用方法

MATLAB が GPU でデータに遭遇した場合、そのデータの計算は GPU で実行されます。クラス gpuArray はデータ配列の作成に使用できる数多くの関数の GPU バージョンを備えており、ここで必要な関数 linspacelogspacemeshgrid も含まれています。同様に、配列 count は関数 ones を使用して GPU で直接、初期化されます。

データ初期化にこれらの変更を加えて、計算を GPU で実行します。

% Setup
t = tic();
x = gpuArray.linspace( xlim(1), xlim(2), gridSize );
y = gpuArray.linspace( ylim(1), ylim(2), gridSize );
[xGrid,yGrid] = meshgrid( x, y );
z0 = complex( xGrid, yGrid );
count = ones( size(z0), 'gpuArray' );

% Calculate
z = z0;
for n = 0:maxIterations
    z = z.*z + z0;
    inside = abs( z )<=2;
    count = count + inside;
end
count = log( count );

% Show
count = gather( count ); % Fetch the data back from the GPU
naiveGPUTime = toc( t );
imagesc( x, y, count )
axis off
title( sprintf( '%1.3fsecs (naive GPU) = %1.1fx faster', ...
    naiveGPUTime, cpuTime/naiveGPUTime ) )

要素単位の演算

アルゴリズムは入力のすべての要素で等しく動作するため、コードを補助関数に配置して arrayfun で呼び出すことができます。GPU 配列の入力では、arrayfun で使用される関数がネイティブの GPU コードにコンパイルされます。ここでは、ループを pctdemo_processMandelbrotElement.m に配置しました。

function count = (x0,y0,maxIterations)
z0 = complex(x0,y0);
z = z0;
count = 1;
while (count <= maxIterations) && (abs(z) <= 2)
    count = count + 1;
    z = z*z + z0;
end
count = log(count);

この関数は 1 つの要素だけを処理するので、早い段階で中止されていることに注意してください。マンデルブロ集合のほとんどの表示では、かなりの数の要素で処理がごく初期に停止し、これによって処理が大幅に軽減されます。また、for ループも while ループに置き換えられていますが、これは後者の効率が通常はより高いからです。この関数は GPU について記述しておらず、GPU の特定の機能も使用しない標準の MATLAB コードです。

arrayfun を使用すると、GPU 用に最適化された別々の演算に対し何千もの呼び出し (反復ごとに少なくとも 6 回) を行う代わりに、すべての計算を実行する並列化された GPU 演算を 1 回呼び出すことになります。これによってオーバーヘッドがかなり軽減されます。

% Setup
t = tic();
x = gpuArray.linspace( xlim(1), xlim(2), gridSize );
y = gpuArray.linspace( ylim(1), ylim(2), gridSize );
[xGrid,yGrid] = meshgrid( x, y );

% Calculate
count = arrayfun( @pctdemo_processMandelbrotElement, ...
                  xGrid, yGrid, maxIterations );

% Show
count = gather( count ); % Fetch the data back from the GPU
gpuArrayfunTime = toc( t );
imagesc( x, y, count )
axis off
title( sprintf( '%1.3fsecs (GPU arrayfun) = %1.1fx faster', ...
    gpuArrayfunTime, cpuTime/gpuArrayfunTime ) );

CUDA の利用

Experiments with MATLAB』では、基本のアルゴリズムを C MEX 関数に変換することでパフォーマンスの向上を実現しています。C/C++ での作業をいとわなければ、Parallel Computing Toolbox を使用して、事前に MATLAB データを使って書いた CUDA カーネルを呼び出すことができます。これは、parallel.gpu.CUDAKernel 機能を用いて実行します。

要素処理アルゴリズムの CUDA/C++ による実装は、pctdemo_processMandelbrotElement.cu に手作業で作成されています。これを nVidia の NVCC コンパイラを使用して手動でコンパイルし、アセンブリレベルの pctdemo_processMandelbrotElement.ptx を生成しなければなりません (.ptx は "並列スレッド実行言語、Parallel Thread eXecution language" を表します)。

この CUDA/C++ コードは、C++ に複素数がないため、これまで見てきた MATLAB バージョンよりやや複雑になります。しかし、アルゴリズムの本質は変わりません。

__device__
unsigned int doIterations( double const realPart0,
                           double const imagPart0,
                           unsigned int const maxIters ) {
   // Initialize: z = z0
   double realPart = realPart0;
   double imagPart = imagPart0;
   unsigned int count = 0;
   // Loop until escape
   while ( ( count <= maxIters )
          && ((realPart*realPart + imagPart*imagPart) <= 4.0) ) {
      ++count;
      // Update: z = z*z + z0;
      double const oldRealPart = realPart;
      realPart = realPart*realPart - imagPart*imagPart + realPart0;
      imagPart = 2.0*oldRealPart*imagPart + imagPart0;
   }
   return count;
}

マンデルブロ集合内の位置に対し 1 つの GPU スレッドが必要で、スレッドがグループ化されてブロックを構成します。カーネルには 1 つのスレッドブロックの大きさが示され、以下のコードでは、これを使用して必要なスレッドブロックの数を計算します。これが GridSize となります。

% Load the kernel
cudaFilename = 'pctdemo_processMandelbrotElement.cu';
ptxFilename = ['pctdemo_processMandelbrotElement.',parallel.gpu.ptxext];
kernel = parallel.gpu.CUDAKernel( ptxFilename, cudaFilename );

% Setup
t = tic();
x = gpuArray.linspace( xlim(1), xlim(2), gridSize );
y = gpuArray.linspace( ylim(1), ylim(2), gridSize );
[xGrid,yGrid] = meshgrid( x, y );

% Make sure we have sufficient blocks to cover all of the locations
numElements = numel( xGrid );
kernel.ThreadBlockSize = [kernel.MaxThreadsPerBlock,1,1];
kernel.GridSize = [ceil(numElements/kernel.MaxThreadsPerBlock),1];

% Call the kernel
count = zeros( size(xGrid), 'gpuArray' );
count = feval( kernel, count, xGrid, yGrid, maxIterations, numElements );

% Show
count = gather( count ); % Fetch the data back from the GPU
gpuCUDAKernelTime = toc( t );
imagesc( x, y, count )
axis off
title( sprintf( '%1.3fsecs (GPU CUDAKernel) = %1.1fx faster', ...
    gpuCUDAKernelTime, cpuTime/gpuCUDAKernelTime ) );

まとめ

この例では、GPU ハードウェアを利用するよう MATLAB アルゴリズムを適合させる 3 つの方法を示しました。

  1. gpuArray を使用して入力データが GPU で使われるように変換し、アルゴリズムは変更しない

  2. arrayfungpuArray 入力で使用し、アルゴリズムを入力の各要素で個別に実行する

  3. parallel.gpu.CUDAKernel を使用して、MATLAB データを使う既存の CUDA/C++ コードを実行する

title('The Mandelbrot Set on a GPU')

この情報は役に立ちましたか?