Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

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

この例では、GPU を使用してマンデルブロ集合を計算するように MATLAB® コードを適合させる方法を示します。

この例では、既存のアルゴリズムから開始し、GPU ハードウェアを利用するように Parallel Computing Toolbox™ を使用してコードを適合させる 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 に記載されているコードに基づいています。tictocを使用して、CPU での実行時間を測定します。

複素数値の 2 次元グリッドを設定します。

tic;
x = linspace(xlim(1),xlim(2),gridSize);
y = linspace(ylim(1),ylim(2),gridSize);

[xGrid,yGrid] = meshgrid(x,y);
z0 = xGrid + 1i*yGrid;

500 回の反復について、前の値を 2 乗して初期値 z0 を加算することで、複素数グリッド z 上の点の次の値を計算します。z の大きさが 2 以下である反復の回数をカウントします。この計算は、すべての位置が同時に更新されるようにベクトル化されています。

cpuCount = ones(size(z0));
z = z0;
for n = 0:maxIterations
    z = z.*z + z0;
    inside = abs(z)<=2;
    cpuCount = cpuCount + inside;
end
cpuCount = log(cpuCount);
cpuTime = toc
cpuTime = 4.4007

カウントの自然対数をプロットします。

figure
imagesc(x,y,cpuCount);
c = colormap([jet;flipud(jet);0 0 0]);
axis off
title(sprintf("CPU Execution: %1.3f s",cpuTime));

gpuArray の使用

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

目的の GPU が使用可能で選択されていることを確認します。

gpu = gpuDevice;
disp(gpu.Name + " GPU selected.")
NVIDIA RTX A5000 GPU selected.

関数 naiveGPUMandelbrot を呼び出します。サポート関数 naiveGPUMandelbrot は GPU でグリッド上の点ごとにマンデルブロ アルゴリズムを適用します。このサポート関数はこの例の最後に記載されています。

[x,y,naiveGPUCount] = naiveGPUMandelbrot(xlim,ylim,gridSize,maxIterations);

gputimeitを使用して、GPU での関数の実行時間を測定します。GPU を使用する関数には、tic および toc、または timeit よりも gputimeit の方が適しています。必ず GPU でのすべての演算が完了してから経過時間を記録するためです。

naiveGPUTime = gputimeit(@() naiveGPUMandelbrot(xlim,ylim,gridSize,maxIterations))
naiveGPUTime = 0.2181

要素単位の演算

アルゴリズムは入力のすべての要素を同様に処理するため、コードを関数に配置してarrayfunで呼び出すことができます。この例の最後に、関数 processMandelbrotElement をサポート関数として示しています。gpuArray の入力では、arrayfun で使用される関数がネイティブの GPU コードにコンパイルされます。

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

arrayfun を使用すると、GPU 用に最適化された個別の演算を数千回 (反復ごとに少なくとも 6 回) 呼び出す代わりに、すべての計算を実行する並列化された GPU 演算を 1 回 MATLAB が呼び出すことになります。GPU で特定の関数を呼び出すために最初に arrayfun を呼び出す時点では、GPU 実行のための関数の設定に若干のオーバーヘッド時間が発生します。同じ関数を指定したその後の arrayfun の呼び出しでは、処理速度が向上します。

複素数値の 2 次元グリッドを設定します。

xGrid = gpuArray(xGrid);
yGrid = gpuArray(yGrid);

arrayfun を使用して、グリッド上の点ごとにマンデルブロ アルゴリズムを適用します。

gpuArrayfunCount = arrayfun(@processMandelbrotElement, ...
                  xGrid,yGrid,maxIterations);

gputimeit を使用して、arrayfun の実行時間を測定します。

gpuArrayfunTime = gputimeit(@() arrayfun(@processMandelbrotElement, ...
                  xGrid,yGrid,maxIterations))
gpuArrayfunTime = 0.0308

CUDA の利用

Experiments with MATLAB では、基本のアルゴリズムを C MEX 関数に変換することでパフォーマンスを向上させています。C/C++ での作業をいとわなければ、Parallel Computing Toolbox を使用して、事前に MATLAB データを使って書いた CUDA カーネルを呼び出すことができます。MATLAB での CUDA カーネルの使用の詳細については、GPU での CUDA または PTX コードの実行を参照してください。

要素処理アルゴリズムの CUDA/C++ による実装は、例 pctdemo_processMandelbrotElement.cu で示されています。単一の場所でマンデルブロ アルゴリズムを実行する CUDA/C++ コードの部分を以下に示します。

__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;
}

mexcudaを使用して、このファイルを並列スレッド実行 (PTX) ファイルにコンパイルします。

mexcuda -ptx pctdemo_processMandelbrotElement.cu
Building with 'NVIDIA CUDA Compiler'.
MEX completed successfully.

CUDA ファイルおよび PTX ファイルを関数 parallel.gpu.CUDAKernel に渡して、parallel.gpu.CUDAKernelオブジェクトを作成します。

cudaFilename = "pctdemo_processMandelbrotElement.cu";
ptxFilename = "pctdemo_processMandelbrotElement.ptx";
kernel = parallel.gpu.CUDAKernel(ptxFilename,cudaFilename);

マンデルブロ集合内の位置ごとに 1 つの GPU スレッドが必要で、スレッドがグループ化されてブロックを構成します。カーネルは 1 つのスレッドブロックの大きさを示します。必要なスレッドブロックの数を計算し、カーネルの GridSize プロパティ (実際上は、GPU によって個別に起動されるスレッド ブロックの数) をそれに従って設定します。

numElements = numel(xGrid);
kernel.ThreadBlockSize = [kernel.MaxThreadsPerBlock,1,1];
kernel.GridSize = [ceil(numElements/kernel.MaxThreadsPerBlock),1];

fevalを使用してカーネルを評価します。

count = zeros(size(xGrid),"gpuArray");
gpuCUDAKernelCount = feval(kernel,count,xGrid,yGrid,maxIterations,numElements);

gputimeit を使用して、カーネルの実行時間を測定します。

gpuCUDAKernelTime = gputimeit(@() feval(kernel,count,xGrid,yGrid,maxIterations,numElements));

まとめ

異なるメソッドの結果をプロットし、実行時間を比較します。

method = ["Naive GPU Execution" "GPU Execution Using arrayfun" "CUDAKernel Execution"];
count = cat(3,naiveGPUCount,gpuArrayfunCount,gpuCUDAKernelCount);
time = [naiveGPUTime gpuArrayfunTime gpuCUDAKernelTime];

figure
colormap(c)
tiledlayout("flow")
nexttile
imagesc(x,y,cpuCount);
axis off
title(sprintf("CPU Execution: %1.3f s",cpuTime));

for idx = 1:3
nexttile
imagesc(x,y,count(:,:,idx))
axis off
title(sprintf("%s: %1.3f s \n  (%1.1fx faster)", ...
    method(idx),time(idx),cpuTime/time(idx)))
end

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

  1. gpuArrayを使用して、GPU 用に入力データを変換し、アルゴリズムは変更しない。

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

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

この例のコードの時間測定は、NVIDIA® RTX A5000 GPU を搭載した Windows® 10、Intel® Xeon® W-2133 3.60 GHz テスト システムで行いました。

サポート関数

サポート関数 naiveGPUMandelbrot は、複素数値の 2 次元グリッドを作成し、複素数値の数 (x0,y0) が複素平面上の半径 2 の円の外側に出るまでの反復回数をカウントします。各反復では z = z^2 + z0 のマッピングが行われます。ここで、z0 = x0 + i*y0 です。この関数は、グリッド座標のベクトルと反復カウントの対数 (外側に出た場合)、または maxIterations (点が外側に出なかった場合) を返します。GPU 上のデータを初期化し、このデータを処理することにより、関数 naiveGPUMandelbrot は GPU 上でマンデルブロ アルゴリズムを実行します。

function [x,y,count] = naiveGPUMandelbrot(xlim,ylim,gridSize,maxIterations)

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");
z = z0;

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

count = log(count);

end

サポート関数 processMandelbrotElement は、複素数値の 2 次元グリッドを作成し、複素数値の数 (x0,y0) が複素平面上の半径 2 の円の外側に出るまでの反復回数をカウントします。各反復では z = z^2 + z0 のマッピングが行われます。ここで、z0 = x0 + i*y0 です。この関数は、反復カウントの対数 (外側に出た場合) または maxIterations (点が外側に出なかった場合) を返します。

function count = processMandelbrotElement(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);

end

参考

|

関連するトピック