描画が3桁速くなると、マンデルブロ集合をインタラクティブに楽しめるようになります。
はじめに
以前、Mathematicaでマンデルブロ集合を描く方法を紹介しました（「追悼、マンデルブロ博士」）。
Clear@f;
f[c_] := f[c, 0., 0]
f[_, _, 100] = -1;
f[c_, z_, i_] := With[{x = z^2 + c}, If[2 < Abs@x, i, f[c, x, i + 1]]]
DensityPlot[f[x + y I], {x, -2, 1}, {y, -1.5, 1.5}, PlotPoints -> 100]
描くのは確かに簡単だったのですが、「ここをもうちょっと拡大したい」という要望に応えるためのインタラクティブな仕掛けは、遅すぎて使えませんでした。「追悼記事なのに・・・」という思いがあったので、ちょっと改善することにしました。
目的
計算を高速化して、インタラクティブなマンデルブロ集合描画システムを作ります。
手法
（1）計算自体（上の例のf）高速化と（2）並列化で高速化します。（1）のために、再帰を反復に書き換えます。（2）のために、CUDAを使います。マンデルブロ集合の計算は、いわゆる「embarrassingly parallel」というやつなので、CUDAの練習に向いていると思います。
この記事のようにインタラクティブなものではありませんが、Mathematicaのマニュアルにもマンデルブロ集合の例が掲載されています。
条件
時間を計る条件は次のようにします。「有理数だと遅い」と思われるかもしれませんが、浮動小数点数を使うと、TableとParallelTableの結果が違ってしまう危険があることへの対応です。ちなみに、これはバグだと思うのですが、Wolfram社の回答は「バグでは無い」でした。
xMin = -2; xMax = 1; yMin = -3/2; yMax = 3/2; dx = (xMax - xMin)/300; dy = (yMax - yMin)/300;
計算機環境
- CPU：Core i7 930（HTオフ）
- 主記憶：6GB
- GPU：NVIDIA GeForce GTX 465
- OS：Windows 7 64 bit
- Mathematica 8.0.1
結果
比較の基準
最初のバージョンをCUDAを使うバージョンで、計算時間を比較してみましょう。
最初のバージョンの実行時間を計ります。
time0 = AbsoluteTiming[
result = Table[ f[x + y I],
{y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]]
23.5279876
CUDA
CUDAの準備をします（MathematicaからCUDAを使う時、書かなければいけないCのコードはこの部分だけなので簡単です）。
Needs["CUDALink`"]
kernelCode = "
__device__ int check(Real_t x, Real_t y, int maxSteps) {
Real_t re = 0, im = 0;
int i = 0;
while (i++ < maxSteps) {
Real_t newRe = re * re - im * im + x;
Real_t newIm = 2 * re * im + y;
if (4 < newRe * newRe + newIm * newIm) return i;
re = newRe;
im = newIm;
}
return -1;
}
__global__ void mandelbrot(int* result, int maxSteps,
int width, int height,
Real_t xStart, Real_t xEnd, Real_t yStart, Real_t yEnd) {
int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
int yIndex = threadIdx.y + blockIdx.y * blockDim.y;
if (xIndex < width && yIndex < height) {
Real_t x = xStart + (xEnd - xStart) * xIndex / width;
Real_t y = yStart + (yEnd - yStart) * yIndex / height;
result[xIndex + yIndex * width] = check(x, y, maxSteps);
}
}
";
kernel = CUDAFunctionLoad[kernelCode, "mandelbrot",
{{_Integer, _, "Output"}, _Integer, _Integer, _Integer,
_Real, _Real, _Real, _Real}, {16, 16}];
これがうまくいかないときは、CUDAのための環境がちゃんとあることを確認してください。私は、Visual Studio 2008 Professionalのデフォルトではx64のコンパイラがインストールされないために少しはまりました。本稿執筆時点でMathematicaが使うデフォルトのCUDA SDKは3.1で、これがサポートするのはVisual Studio 9 (2008)までだということにも注意が必要です（カスタマイズすればいいのでしょうが）。追記：2011年8月にCUDA SDK 4に対応したので、Visual Studio 2010がそのまま使えるようになりました。
実行してみると、4000倍くらい速くなることがわかります。
time = AbsoluteTiming[
width = 300; height = 300;
buffer = CUDAMemoryAllocate[Integer, {height, width}];
kernel[buffer, 100, width, height, N@xMin, N@xMax, N@yMin, N@yMax];
result = CUDAMemoryGet[buffer];][[1]];
CUDAMemoryUnload[buffer];
time0/time
3920.8
まあ、GPU無しの比較対象は最初に思いつくものですから、「まずGPU無しでどこまで速くなるのか追求してから」と思う方は、下の「おまけ」を読んでください。GPU無しでも最初のものより80倍くらいは速くできるので、それと比べれば、「GPUの効果は50倍くらい」というのが妥当なところでしょう。
リアルタイム描画
描画には、ArrayPlotを使うことにします。ListDensityPlotのほうが高度ですが、少し遅いのと、座標の扱いに問題があったため使えませんでした。Imageは速そうですが、やはり座標の扱いに問題がありました。
ユーザインターフェースはMathematicaのManipulateで構築します。
width = 300; height = 300;
buffer = CUDAMemoryAllocate[Integer, {height, width}];
Manipulate[With[{
xMin = center[[1]] - 10^scale, xMax = center[[1]] + 10^scale,
yMin = center[[2]] - 10^scale, yMax = center[[2]] + 10^scale},
kernel[buffer, IntegerPart[10^maxSteps], width, height, xMin, xMax,
yMin, yMax];
ArrayPlot[CUDAMemoryGet[buffer],
DataRange -> {{xMin, xMax}, {yMin, yMax}},
DataReversed -> True, ColorFunction -> "Rainbow"]],
{{center, {-0.5, 0}}, Locator},
{{scale, 0, "Scale"}, -13, 1},
{{maxSteps, 2, "Log[maxSteps]"}, 1, 3}]
画像上をクリックすると、その場所を中心にして再描画します。Scaleで描画領域の広さを、Log[maxSteps]で反復の上限を変更できます。
ぬるぬる動くというわけにはいきませんが、「ここをもうちょっと拡大したい」という要望には十分応えられると思います。もっと速くしたい場合は、ArrayPlotをどうにかするということになるのでしょうが、もう少しコードを書かなければならないでしょう。
ところで、さすがにこれってCDFにしてもだめなんですよねえ。
おまけ（細かいことなので、これ以降は読まなくていいです）
追悼記事がダメだった理由
追悼記事がダメだった最大の要因は、DensityPlotにあったようです。DensityPlotは、下に示すように場所ごとにメッシュの細かさを自動調整してくれる便利な関数なのですが、この内部が並列化されていないため、あまり速くならないのです。
CUDAなしでどこまで行くか
コンパイル＋並列化
Mathematicaでは、機械精度の数を使うことに限定すれば、ユーザ定義関数をコンパイルすることができます。最初のバージョンをコンパイルし、並列化して実行すると30倍くらい速くなります。
fc = Compile[{{c, _Complex}, {z, _Complex}, {i, _Integer}},
If[i == 100, -1,
With[{x = z^2 + c}, If[2 < Abs@x, i, fc[c, x, i + 1]]]],
CompilationTarget -> "C"];
DistributeDefinitions[fc];
time = AbsoluteTiming[
result = ParallelTable[ fc[x + y I, 0, 0],
{y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]];
time0/time
28.44619
反復
再帰を反復に書き直すともっと速く、80倍くらい速くなります。
hc = Compile[{{z, _Complex}},
Module[{c = 0. + 0. I, i = 0},
While[i++ < 100,
c = c^2 + z;
If[2 < Abs@c, Return@i]];
Return@-1]];
DistributeDefinitions[hc];
time = AbsoluteTiming[
result = ParallelTable[ hc[N@x + N@y I],
{y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]];
time0/time
77.13116
ちなみに、この例ではコンパイル時に「CompilationTarget -> “C”」を付けてもあまり効果は無いようです。
CUDA無しでのリアルタイム描画
この方法を使ってリアルタイムに描画してみます。スケールや反復の上限もリアルタイムに変更できるように、少し修正します。
hc2 = Compile[{{xIndex, _Integer}, {yIndex, _Integer},
{maxSteps, _Integer},
{width, _Integer}, { height, _Integer},
{xMin, _Real}, {xMax, _Real},
{yMin, _Real}, {yMax, _Real}},
With[{z = xMin + (xMax - xMin) xIndex/width
+ (yMin + (yMax - yMin) yIndex/height) I},
Module[{c = 0. + 0. I, i = 0},
While[i++ < maxSteps,
c = c^2 + z;
If[4 < Abs@c, Return@i]];
Return@-1]],
CompilationTarget -> "C"];
DistributeDefinitions[hc2];
width = 300; height = 300;
Manipulate[
With[{
xMin = center[[1]] - 10^scale, xMax = center[[1]] + 10^scale,
yMin = center[[2]] - 10^scale, yMax = center[[2]] + 10^scale,
maxSteps = IntegerPart[10^logMaxSteps]},
result = ParallelTable[
hc2[xIndex, yIndex, maxSteps, width, height, xMin, xMax, yMin, yMax],
{yIndex, 1, height}, {xIndex, 1, width}];
ArrayPlot[result, DataRange -> {{xMin, xMax}, {yMin, yMax}},
DataReversed -> True, ColorFunction -> "Rainbow"]],
{{center, {-0.5, 0}}, Locator},
{{scale, 0, "Scale"}, -13, 1},
{{logMaxSteps, 2, "Log[maxSteps]"}, 1, 3}]
実行結果は上に掲載したものと同じです。速いCPUがある場合はこれで十分なのかもしれませんが、私の環境では、「一応動くがストレスフル」という感じでした。
