適用例—Wolfram Documentation

適用例

画像処理	乱数生成
線形代数とリスト処理	コードの生成

GPUはSIMDマシンなので，CUDAの潜在力を利用するにはSIMDの方法で問題を提示しなければならない．各スレッドが1つの要素を独立して計算できる方法で分割できる計算はGPUに理想的である．

アルゴリズムによっては並列で書けなかったり，CUDA上で使えなかったり（アーキテクチャの制約のため）する．そのような場合に，それらの計算を行うためにGPUを使用する別の方法を導入するための研究が行われている．

このセクションでは，Wolfram言語内におけるCUDAプログラミングの使い方を示す．以下のすべての例において，CUDAソース，バイナリまたはライブラリをWolfram言語にロードすることができるようになるCUDAFunctionLoadが使われる．

CUDAFunctionLoad

CUDA関数をWolfram言語にロードする

CUDAFunctionLoadを使うと，CUDAのソース，バイナリ，ライブラリをWolfram言語にロードすることが可能になる

画像処理

このセクションでは画像処理操作を行うCUDAの適用例を示す．CUDALink にはCUDAImageConvolve，CUDABoxFilter，CUDAErosion，CUDADilation，CUDAOpening，CUDAClosing，CUDAImageAdd等，画像処理を行う関数がいくつか組み込まれている．

画像の二値化

Binarizeは入力画像を取り，画素が閾値より上の場合は白，それ以外は黒である二値画像を出力する．

CUDALink アプリケーションをまだインポートしていない場合はインポートする．

Wolfram Language code: Needs["CUDALink`"]

入力画像を定義する．GPU上の必要なメモリ量を減らすために，画像を表すのにを使う．

Wolfram Language code:

binarizeImageCode = "__global__ void binarize(unsigned char * in, unsigned char * out, mint threshold, mint width, mint height, mint channels) {
   	mint xIndex = threadIdx.x + blockIdx.x*blockDim.x;
   	mint yIndex = threadIdx.y + blockIdx.y*blockDim.y;
   	mint index = channels*(xIndex + yIndex*width);
   	if (xIndex < width && yIndex < height) {
  		 mint accum = 0;
   		for (mint ii = 0; ii < channels; ii++)
   			accum += in[index+ii];
  		 out[xIndex + yIndex*width] = accum > channels*threshold ? 255 : 0;
   	}
}";

CUDAFunctionをロードする．型はパラメータリスト中で"UnsignedByte"で表される．

Wolfram Language code:

binarizeImage = CUDAFunctionLoad[binarizeImageCode, "binarize", {{"UnsignedByte", "Input"}, {"UnsignedByte", "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}]

入力画像を定義し，出力に"UnsignedByte"メモリを割り当てる．

Wolfram Language code:

input = [image];
{height, width} = ImageDimensions[input];
channels = ImageChannels[input];
output = CUDAMemoryAllocate["UnsignedByte", {width, height}];

閾値を150として二値化関数を呼び出す．

Wolfram Language code: binarizeImage[input, output, 150, width, height, channels, {width, height}]

出力画像を表示する．

Wolfram Language code: Image[output, ImageSize -> Medium]

結果はWolfram言語のものと一致する．

Wolfram Language code: Binarize[[image]]

割り当てられたメモリをアンロードする．

Wolfram Language code: CUDAMemoryUnload[output]

ボックスフィルタ

ぼっくするフィルタはカーネルがBoxMatrixの時に最適化された畳込みである．ここでこれを実装する．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "boxFilter.cu"}]

ボックスフィルタを適用するのに必要なCUDA関数をロードする．

Wolfram Language code:

boxFilterHorizontal = CUDAFunctionLoad[{srcf}, "d_boxfilter_rgba_x", {{"UnsignedByte", "Input"}, {"UnsignedByte", "Output"}, _Integer, _Integer, _Integer}, {16, 1}, "IncludeDirectories" -> FileNameJoin[{$CUDALinkPath, "SupportFiles"}], "UnmangleCode" -> False];
boxFilterVertical = CUDAFunctionLoad[{srcf}, "d_boxfilter_rgba_y", {{"UnsignedByte", "Input"}, {"UnsignedByte", "Output"}, _Integer, _Integer, _Integer}, {16, 1}, "IncludeDirectories" -> FileNameJoin[{$CUDALinkPath, "SupportFiles"}], "UnmangleCode" -> False];

入力パラメータを設定する．

Wolfram Language code:

input = CUDAMemoryLoad[[image], "UnsignedByte"];
tmp = CUDAMemoryAllocate["UnsignedByte", {512, 512, 4}];
output = CUDAMemoryAllocate["UnsignedByte", {512, 512, 4}];
width = height = 512;

半径は5に設定される．

Wolfram Language code: radius = 5;

関数を呼び出す．

Wolfram Language code:

boxFilterHorizontal[input, tmp, width, height, radius];
boxFilterVertical[tmp, output, width, height, radius];

画像を得る．

Wolfram Language code: Image[CUDAMemoryGet[tmp], "Byte", ImageSize -> Small]

割り当てられたメモリを解放する．

Wolfram Language code: CUDAMemoryUnload[input, output, tmp]

画像の調整

これはCUDAにおけるImageAdjustの実装である．

Wolfram Language code:

src = "
__device__ mint xclamp(mint val, mint low, mint high) {
	return val <= low ? low : (val >= high ? high : val);
}

__device__ mint adjust(mint pixel, float lowIn, float highIn, float lowOut, float highOut, float gamma) {

	float res, val;
	val = xclamp(pixel, lowIn, highIn);
	
	res = pow((val - lowIn) / (highIn - lowIn), gamma);
	res = res * (highOut - lowOut) - lowOut;
	
	return res + 0.5f;
}

__global__ void imageAdjust(mint * img, mint width, mint height, mint channels, float lowIn, float highIn, float lowOut, float highOut, float gamma) {
	
	int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
	int yIndex = threadIdx.y + blockIdx.y * blockDim.y;
	if (xIndex >= width || yIndex >= height)
		return ;
	int pos = channels * (yIndex * width + xIndex);
	for (mint ii = 0; ii < channels; ii++) {
		img[pos + ii] = adjust(img[pos + ii], 255*lowIn, 255*highIn, 255*lowOut, 255*highOut, gamma);
	}
}";

ソース文字列からCUDAFunctionをロードし，定数値にFloatを使用する．

Wolfram Language code:

cCUDAImageAdjust = CUDAFunctionLoad[src, "imageAdjust", {{_Integer}, _Integer, _Integer, _Integer, "Float", "Float", "Float", "Float", "Float"}, {16, 16}]

CUDAFunctionをラップし，CUDAImageAdjustをImageAdjustと似たシンタックスで作る．

Wolfram Language code:

CUDAImageAdjust[img_Image, {lowIn_Real, highIn_Real}, gamma_ : 1.0] /; Head[gamma] == Real := CUDAImageAdjust[img, {lowIn, highIn}, {0.0, 1.0}, gamma]
CUDAImageAdjust[img_Image, {lowIn_Real, highIn_Real}, {lowOut_Real, highOut_Real}, gamma_ : 1.0] := 
	Module[{width, height, channels}, 
	{height, width, channels} = Flatten[{ImageDimensions[img], ImageChannels[img]}];
	cCUDAImageAdjust[img, width, height, channels, lowIn, highIn, lowOut, highOut, gamma, {width, height}]//First
	]

関数はImageAdjustと同じように使える．

Wolfram Language code: CUDAImageAdjust[[image], {0.3, 0.8}]

キャニー(Canny)エッジ検出

キャニーエッジ検出は10種類程度のフィルタを組み合せ，画像のエッジを見付ける．Wolfram言語のEdgeDetectは同様の機能を提供する．以下がその実装である．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "canny.cu"}]

エッジ検出を行うのに必要なCUDA関数をロードする．

Wolfram Language code:

gaussianVertical = CUDAFunctionLoad[{srcf}, "gaussianVert_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
gaussianHorizontal = CUDAFunctionLoad[{srcf}, "gaussianHoriz_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
sobelGXVertical = CUDAFunctionLoad[{srcf}, "sobelGxVert_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
sobelGXHorizontal = CUDAFunctionLoad[{srcf}, "sobelGxHoriz_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
sobelGYVertical = CUDAFunctionLoad[{srcf}, "sobelGyVert_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
sobelGYHorizontal = CUDAFunctionLoad[{srcf}, "sobelGyHoriz_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
magnitude = CUDAFunctionLoad[{srcf}, "magnitude_kernel", {{_Integer, "Input"}, {_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
atan = CUDAFunctionLoad[{srcf}, "atan_kernel", {{_Integer, "Input"}, {_Integer, "Input"}, {"Float", "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
zeroCrossing = CUDAFunctionLoad[{srcf}, "zeroCrossing_kernel", {{"Float", "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
nonMaximalSupression = CUDAFunctionLoad[{srcf}, "nonMaximalSuppression_kernel", {{_Integer, "Input"}, {_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
hysteresis = CUDAFunctionLoad[{srcf}, "hysteresis_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];
binarize = CUDAFunctionLoad[{srcf}, "binarize_kernel", {{_Integer, "Input"}, {_Integer, "Output"}, _Integer, _Integer, _Integer, _Integer}, {16, 16}];

入力画像とパラメータを設定する．

Wolfram Language code:

img = [image];
{{width, height}, channels} = {ImageDimensions[img], ImageChannels[img]};
pitch = width * channels;

ここでCUDAマネージャにホストとデバイスを加える．これらは入力と出力の両方を保持する．

Wolfram Language code:

inhost = CUDAMemoryLoad[img, Integer];
outhost = CUDAMemoryAllocate[Integer, {width, height, channels}];

次に計算に使うためのデバイスのみの一時メモリを定義する．キャニーエッジ検出には多くのフィルタを使うので，多くのメモリが必要である．

Wolfram Language code:

tmpdev = CUDAMemoryAllocate[Integer, {width, height, channels}];
gxdev = CUDAMemoryAllocate[Integer, {width, height, channels}];
gydev = CUDAMemoryAllocate[Integer, {width, height, channels}];
cdev = CUDAMemoryAllocate[Integer, {width, height, channels}];
magdev = CUDAMemoryAllocate[Integer, {width, height, channels}];
tandev = CUDAMemoryAllocate["Float", {width, height, channels}];

キャニーエッジ検出関数を呼び出す．

Wolfram Language code:

gaussianVertical[inhost, tmpdev, width, height, channels, pitch];
gaussianHorizontal[tmpdev, outhost, width, height, channels, pitch];
sobelGXHorizontal[outhost, tmpdev, width, height, channels, pitch];
sobelGXVertical[tmpdev, gxdev, width, height, channels, pitch];
sobelGYHorizontal[outhost, tmpdev, width, height, channels, pitch];
sobelGYVertical[tmpdev, gydev, width, height, channels, pitch];
magnitude[gxdev, gydev, magdev, width, height, channels, pitch];
atan[gxdev, gydev, tandev, width, height, channels, pitch];
zeroCrossing[tandev, cdev, width, height, channels, pitch];
nonMaximalSupression[magdev, cdev, tmpdev, width, height, channels, pitch];
hysteresis[tmpdev, outhost, width, height, channels, pitch];
binarize[outhost, outhost, width, height, channels, pitch];

出力を画像として表示する．

Wolfram Language code: Image[outhost]

割り当てられたメモリをアンロードする．

Wolfram Language code: CUDAMemoryUnload[inhost, outhost, tmpdev, gxdev, gydev, cdev, magdev, tandev]

線形代数とリスト処理

このセクションでは線形代数操作を行うCUDA適用例を示す．ここで取り上げる関数の多くはCUDATransposeかCUDADotを使って実行できる．

行列の転置

行列の転置は多くのアルゴリズムにおいて不可欠である．CUDALink はCUDATransposeという形ですぐに使える実装を提供する．しかし，ユーザは独自のものを実装することもできる．

行列の転置のためにCUDAFunctionをロードし，実数行列を取ってその転置を出力する新しい関数newCUDATransposeを定義する．

Wolfram Language code:

newCUDATranspose[matrix_] := Module[
	{dInputMatrix, dOutputMatrix, hOutputMatrix, blockDim, transposeFun}, 
	
	blockDim = {16, 16};
	
	dInputMatrix = CUDAMemoryLoad[matrix];
	dOutputMatrix = CUDAMemoryAllocate[Real, Reverse@ Dimensions[matrix]];
	
	transposeFun = CUDAFunctionLoad[{FileNameJoin[{$CUDALinkPath, "SupportFiles", "transpose.cu"}]}, "transpose_kernel", {{_Real, "Input"}, {_Real, "Output"}, _Integer, _Integer}, blockDim];
	
	transposeFun[dInputMatrix, dOutputMatrix, Sequence@@Reverse[Dimensions[matrix]]];
	
	hOutputMatrix = CUDAMemoryGet[dOutputMatrix];
	
	CUDAMemoryUnload[dInputMatrix];
	CUDAMemoryUnload[dOutputMatrix];
	
	Return[hOutputMatrix];
	];

行列を設定する．

Wolfram Language code: MatrixForm[A = Table[1.0 * i - j, {i, 0, 3}, {j, 0, 8}]]

行列Aを転置する．

Wolfram Language code: MatrixForm[newCUDATranspose[A]]

Wolfram言語の結果と一致する．

Wolfram Language code: MatrixForm[Transpose[A]]

行列とベクトルの乗算

行列とベクトルの乗算は線形代数，有限要素解析等において一般的な操作である．CUDAFunctionをロードし，行列とベクトルの乗算を実装する．

Wolfram Language code:

CUDAMatrixVectorMultiply[matrix_, vector_] := Module[
	{pmv, blockDim, gridDim, dInputMatrix, dInputVector, dOutputVector, hOutputVector}, 
	
	blockDim = 16;
	gridDim = First@Dimensions[matrix];
	
	dInputMatrix = CUDAMemoryLoad[matrix];
	dInputVector = CUDAMemoryLoad[vector];
	dOutputVector = CUDAMemoryAllocate[Integer, Length@vector];
	
	pmv = CUDAFunctionLoad[{FileNameJoin[{$CUDALinkPath, "SupportFiles", "matrixVectorMul.cu"}]}, "matrixVecMul_kernel", {{_Integer, _, "InputOutput"}, {_Integer, _, "Input"}, {_Integer, _, "Input"}, _Integer, _Integer}, blockDim];
	
	pmv[dOutputVector, dInputMatrix, dInputVector, Sequence@@Reverse[Dimensions[matrix]], gridDim];
	
	hOutputVector = CUDAMemoryGet[dOutputVector];
	
	CUDAMemoryUnload /@ {dOutputVector, dInputMatrix, dInputVector};
	
	Return[hOutputVector];
	];

入力行列，ベクトルを設定する．

Wolfram Language code:

A = Table[i + j, {i, 0, 15}, {j, 0, 15}];
B   = Table[i, {i, 0, 15}];

上で定義した関数を呼び出し，MatrixFormを使って結果を表示する．

Wolfram Language code: MatrixForm[CUDAMatrixVectorMultiply[A, B]]

Wolfram言語の結果と一致する．

Wolfram Language code: MatrixForm[A.B]

行列と行列の乗算

行列と行列の乗算は多くのアルゴリズムにおいて重要な関数である．ソースファイルからCUDAFunctionをロードし，ブロック次元を4に設定する．

Wolfram Language code:

blockSize = 4;
MatrixMultiply = CUDAFunctionLoad[{FileNameJoin[{$CUDALinkPath, "SupportFiles", "matrixMul.cu"}]}, "matrixMul", {{_Real, 2, "Output"}, {_Real, 2, "Input"}, {_Real, 2, "Input"}, _Integer, _Integer}, {blockSize, blockSize}, "Defines" -> {"BLOCK_SIZE" -> blockSize}];

入力値を設定し，出力のためのメモリを割り当てる．

Wolfram Language code:

A = RandomReal[1.0, {8, 8}];
B = RandomReal[1.0, {8, 8}];
out = CUDAMemoryAllocate[Real, {8, 8}];

計算を行う．

Wolfram Language code: MatrixMultiply[out, A, B, Sequence@@Dimensions[A]];

MatrixFormを使って結果を表示する．

Wolfram Language code: CUDAMemoryGet[out]//MatrixForm

Wolfram言語の結果と一致する．

Wolfram Language code: Dot[A, B]//MatrixForm

出力バッファをアンロードする．

Wolfram Language code: CUDAMemoryUnload[out]

ドット積

2つのベクトルのドット積は線形代数において一般的な操作である．ベクトルの集合を取り，それぞれのドット積を返す関数を実装する．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "dotProduct.cu"}]

関数をWolfram言語にロードする．

Wolfram Language code:

dotProduct = CUDAFunctionLoad[{srcf}, "scalarProdGPU", {{_Real, "Output"}, {_Real, "Input"}, {_Real, "Input"}, _Integer, _Integer}, 1024]

50のランダムベクトルを生成する．

Wolfram Language code:

A = RandomReal[1.0, {50, 50}];
B = RandomReal[1.0, {50, 50}];
out = ConstantArray[0.0, First[Dimensions[A]]];

関数を呼び出し，結果を返す．

Wolfram Language code: dotProduct[out, A, B, Sequence@@Dimensions[A], Length[out]]//First

結果はWolfram言語の出力と一致する．

Wolfram Language code: MapThread[Dot, {A, B}]

凸包

凸包はGPU上で並列化が難しいため，この例ではGPUプログラミングへのハイブリッドアプローチを取る．つまり，GPUに適した計算をGPU上で行い，残りをCPUで行う．

この凸包の実装はAndrewのアルゴリズムの教科書通りの実装で，簡単になるように設計されており，効率的ではない．

Wolfram言語内でLeftOf述語を定義する．

Wolfram Language code:

pdirection[pa_List, pb_List, px_List] := ((pa[[1]] - px[[1]]) * (pb[[2]] - px[[2]])) - ((pb[[1]] - px[[1]]) * (pa[[2]] - px[[2]]))

端点を結ぶ与えられた線の上か下かどちらかに点集合を分割する関数を定義する．

Wolfram Language code:

code = "
__global__ void partitionPts(float *inpts, mint *idx, mint length) {
	int index = threadIdx.x + blockIdx.x*blockDim.x;
	int length1 = length-1;
	int index2, tindex;
	Real_t  leftx, lefty;
	Real_t  rightx, righty;
	Real_t  xx, yy, det;

	if(index == 0 || index == length1) {
		idx[index] = -1;
		return;
	}
	if(index >= length)
		return;

	tindex = 2*length1;
	index2 = 2*index;
	leftx = inpts[0];
	lefty = inpts[1];
	rightx = inpts[tindex];
	righty = inpts[tindex+1];

	xx = inpts[index2];
	yy = inpts[index2+1];

	det = ((leftx - xx)*(righty - yy)) - ((rightx - xx)*(lefty - yy));

	if(det >=0)
		idx[index] = 1;
	else
		idx[index] = 0;
}";

上記関数をCUDA関数としてロードする．

Wolfram Language code: partitionPts = CUDAFunctionLoad[code, "partitionPts", {{"Float[2]", _, "Input"}, {_Integer}, _Integer}, 32];

2つの分割点集合を取り，それらの凸包を求める．

Wolfram Language code:

makeHalfHull[lcoord_List, ppart_List, factor_] := 
	Module[{hull, length, ii, ihull, idx, pa, pb, px, ptop}, 
		length = Length[ppart];
		hull = ConstantArray[0, length];
		hull[[1]] = 1;
	
		ii = 2; ihull = 1;
		While[ii ≤ length, 
			ihull++;
			hull[[ihull]] = ppart[[ii]];
			While[ihull ≥ 3, 
				idx = {hull[[ihull - 2]], hull[[ihull]], hull[[ihull - 1]]};
				{pa, pb, px} = Part[lcoord, idx];
				If[factor * pdirection[pa, pb, px] < 0.0, 
					hull[[ihull - 1]] = hull[[ihull]];
					hull[[ihull]] = 0;
					ihull--, 
					Break[]
				]
			];
			ii++
		];
		ptop = Position[hull, 0];
		If[ptop  == {}, Return[hull]];
		Take[hull, First[ptop][[1]] - 1]
	]

上記関数を呼び出す，リストはpartitionPtsが処理する前にCUDASortでソートされる．

Wolfram Language code:

CUDAConvexHull[pts_List] := 
	Module[{mem, slst, ilst, imem, ilength, ppts, upts, lpts, lhull, uhull, ihull}, 
		ilength = Length[pts];
	
		(* CUDA Sort points *)
		mem = CUDAMemoryLoad[pts, "Float[2]"];
		CUDASort[mem];
		slst = CUDAMemoryGet[mem];
	
		(* CUDA Classify points as bottom/top, upts and lpts *)
		ilst = ConstantArray[0, ilength];
		imem = CUDAMemoryLoad[ilst, Integer];
		partitionPts[mem, imem, ilength, ilength];
		ppts = CUDAMemoryGet[imem];
		ppts[[1]] = 1;ppts[[ilength]] = 1;
		upts = Flatten[Position[ppts, 1]];
		ppts[[1]] = 0;ppts[[ilength]] = 0;
		lpts = Flatten[Position[ppts, 0]];
	
		(* Get lower and upper half hulls *)
		lhull  = makeHalfHull[slst, lpts, -1.0];
		uhull  = makeHalfHull[slst, upts, 1.0];
		
		(* Construct Final ConvexHull *)
		If[Length[lhull] == 0, 
			ihull = uhull, 
			If[Length[uhull] == 0, 
			ihull = lhull, 
			ihull = Join[lhull, Reverse[Drop[uhull, -1]]]
		]
		];
	
		(* Delete CUDA Memory *)
		CUDAMemoryUnload[mem, imem];
	
		(* Get result *)
		Part[slst, ihull]
		]

検証するために，一様に分布したランダムな20,000の点を作成する．

Wolfram Language code:

npts = 20000;
lst = RandomReal[10, {npts, 2}];

凸包を計算する．

Wolfram Language code: hullpts = CUDAConvexHull[lst];

結果を可視化する．凸包の点の間に線が引かれる．

Wolfram Language code: Graphics[{Green, Point[lst], Magenta, Line[hullpts]}]

上記アルゴリズムはすべて，またはほとんどの点が凸包に含まれる極端な場合も扱う．単位円板じょうに一様に分布した点を生成する．

Wolfram Language code:

npts = 500;
lst = Table[pt = RandomReal[{-10, 10}, 2];pt / Norm[pt], {i, npts}];

凸包の点を計算する．

Wolfram Language code: hullpts = CUDAConvexHull[lst];

凸包を可視化する．

Wolfram Language code: Graphics[{Green, Point[lst], Magenta, Line[hullpts]}]

これは直列コードでしか書けないようなアルゴリズムを，Wolfram言語とCUDAプログラミングを組み合せて部分的に並列にする重要な例である．

乱数生成

レイトレーシングから偏微分方程式の解放に至るまでの多くのアルゴリズムで，入力に乱数が必要となる．しかしこれは各コアが他のコアと同じ状態である複数のコアのシステムの多くでは難しい問題である．それを避けるために，通常は並列乱数生成器は一日の中の時間等のエントロピー値を乱数生成器のシードに利用するが，これらの呼出しはCUDAでは利用できない．

次のセクションでは乱数生成の3クラスのアルゴリズムを取り上げる．最初は一様乱数生成器である（擬似乱数生成器と準乱数生成器を例示する）．2つ目は乱数生成にハッシング関数の一様分布を利用する乱数生成器，そして最後は正規乱数生成器である．

擬似乱数生成器

擬似乱数生成器はランダムであるように見える数を生成する決定論アルゴリズムである．これはシード値を利用して乱数を生成する．

このセクションでは，簡単な線形合同乱数生成器（Park–Millerのアルゴリズム）と，より複雑なメルセンヌツイスタ(Mersenne Twister)を取り上げる．

Park–Miller

Park–Miller乱数生成器は次の再帰方程式で定義される．

これはWolfram言語で簡単に実装できる．

Wolfram Language code: ParkMiller[x_, a_ : 16807, n_ : 2147483647] := Mod[x * a, n]

とに共通の値を使う．NestListを使って1000個の数のリストを生成し，プロットすることができる．

Wolfram Language code: ListPlot[NestList[ParkMiller, 1, 1000]]

百万個の数の生成にかかる時間である．

Wolfram Language code: AbsoluteTiming[NestList[ParkMiller, 1, 10000000];]

代りにSeedRandomでMethodオプションを使うこともできる．これも前回と同じように使える．

Wolfram Language code:

ListPlot[BlockRandom[SeedRandom[1, Method -> {"Congruential", "Multiplier" -> 16807, "Increment" -> 0, "Modulus" -> 2147483647}];RandomInteger[2 ^ 31 - 1, 1000]]]

Wolfram言語での実装と比較すると，これは300倍も速い．

Wolfram Language code:

AbsoluteTiming[BlockRandom[SeedRandom[1, Method -> {"Congruential", "Multiplier" -> 16807, "Increment" -> 0, "Modulus" -> 2147483647}];RandomInteger[2 ^ 31 - 1, 10000000]];]

CUDAでの実装はWolfram言語で書かれたものと似ている．実装は CUDALink で分配され，以下に場所を示す．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

CUDA関数をWolfram言語にロードする．

Wolfram Language code:

CUDAParkMiller = CUDAFunctionLoad[{srcf}, "ParkMiller", {{_Integer, _, "Input"}, {_Integer, _, "Output"}, _Integer}, 256]

出力メモリを割り当てる．

Wolfram Language code: mem = CUDAMemoryAllocate[Integer, {2048}]

CUDAFunctionを呼び出す．

Wolfram Language code: CUDAParkMiller[RandomInteger[Developer`$MaxMachineInteger, 256], mem, 2048]

結果はランダムである．

Wolfram Language code: ListPlot[CUDAMemoryGet[mem]]

時間を測定すると，Wolfram言語の組込みメソッドより2倍，純粋な Wolfram言語の実装より600倍速いことが分かる．

Wolfram Language code: mem = CUDAMemoryAllocate[Integer, {10000000}]

Wolfram Language code: AbsoluteTiming[CUDAParkMiller[RandomInteger[Developer`$MaxMachineInteger, 256], mem, 10000000];]

Compileに対する計算時間もCUDAのものと同程度である．Compile文から整数オーバーフロー検出のないCコードを生成する．

Wolfram Language code:

CompiledParkMiller = Compile[{{n, _Integer}}, NestList[# * 16807&, 1, n], RuntimeOptions -> {"CatchMachineIntegerOverflow" -> False}, CompilationTarget -> "C"];

計算時間を求める．少し違いがある．

Wolfram Language code: AbsoluteTiming[CompiledParkMiller[50000000];]

メモリの割当てに必要な時間を無視すると（生成された乱数はGPU上で再利用されるので，これは正当である），速度は10倍改善されたことが分かる．

Wolfram Language code:

mem = CUDAMemoryAllocate[Integer, {10000000}];
CUDAMemoryCopyToDevice[mem];
AbsoluteTiming[CUDAParkMiller[RandomInteger[Developer`$MaxMachineInteger, 512], mem, 10000000];]

メルセンヌツイスタ

メルセンヌツイスタは乱数の生成にシフトレジスタを利用する．実装は簡単なので，GPUにうまくマップされる．以下のファイルに実装が含まれている．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

ファイルから"MersenneTwister"関数をロードする．

Wolfram Language code:

mersenneTwister = CUDAFunctionLoad[{srcf}, "MersenneTwister", {{_Real, _, "Output"}, {_Integer, _, "Input"}, {_Integer, _, "Input"}, {_Integer, _, "Input"}, {_Integer, _, "Input"}, _Integer}, 32];

メルセンヌツイスタ入力パラメータを定義する．つい巣他派シード値を必要とする．これらの値はオフラインで計算してファイルに保管したり，Wolfram言語で生成したりすることができる．ここでは後者を示す．

Wolfram Language code:

MTRNGCount = 4096;
PATHN = 2 ^ 14;
NPerRNG = Ceiling[PATHN / MTRNGCount];
NPerRNG = If[EvenQ[NPerRNG], NPerRNG, NPerRNG + 1];
RANDN = MTRNGCount * NPerRNG;
{hsMatrixA, hsMaskB, hsMaskC} = RandomInteger[{-Developer`$MaxMachineInteger, Developer`$MaxMachineInteger}, {3, MTRNGCount}];
hsSeed = RandomInteger[{-Developer`$MaxMachineInteger, Developer`$MaxMachineInteger}, MTRNGCount];

出力メモリを割り当てる．出力は上書きされるので，Wolfram言語からGPU上にメモリをロードする必要はない．

Wolfram Language code: output = CUDAMemoryAllocate[Real, RANDN]

CUDAFunctionをパラメータを与えて呼び出す．

Wolfram Language code: mersenneTwister[output, hsMatrixA, hsMaskB, hsMaskC, hsSeed, NPerRNG, MTRNGCount]

出力をプロットするとランダムであることが示される．

Wolfram Language code: ListPlot[CUDAMemoryGet[output]]

生成する乱数を入力として取り，パラメータを設定して必要な割当てを行い，ランダムな出力メモリを返すWolfram言語関数を書くことができる．

Wolfram Language code:

MersenneTwister[n_] := Module[{MTRNGCount, PATHN, NPerRNG, RANDN, hsMatrixA, hsMaskB, hsMaskC, hsSeed, output}, 
	MTRNGCount = 4096;
	PATHN = n;
	NPerRNG = Ceiling[PATHN / MTRNGCount];
	NPerRNG = If[EvenQ[NPerRNG], NPerRNG, NPerRNG + 1];
	RANDN = MTRNGCount * NPerRNG;
	{hsMatrixA, hsMaskB, hsMaskC} = RandomInteger[{-Developer`$MaxMachineInteger, Developer`$MaxMachineInteger}, {3, MTRNGCount}];
	hsSeed = RandomInteger[{-Developer`$MaxMachineInteger, Developer`$MaxMachineInteger}, MTRNGCount];
	output = CUDAMemoryAllocate[Real, RANDN];
	First@mersenneTwister[output, hsMatrixA, hsMaskB, hsMaskC, hsSeed, NPerRNG, 128]
	]

最初の10,000要素のプロットを生成する．

Wolfram Language code: ListPlot[CUDAMemoryGet[MersenneTwister[10000]][[ ;; 10000]]]

一億個の数を精製するのに乱数生成器が必要とする時間を測定する．

Wolfram Language code: MersenneTwister[100000000];//AbsoluteTiming

これはWolfram言語の乱数生成器が必要とする時間と同程度である．

Wolfram Language code: BlockRandom[SeedRandom[1, Method -> "MersenneTwister"];RandomReal[1, 100000000]];//AbsoluteTiming

乱数が他の問題のシードであると考えると，Wolfram言語の所要時間がCUDA実装のものより優れていても，アルゴリズム全体としてパフォーマンスが改善される．

準乱数生成器

このセクションでは準乱数生成器について述べる．擬似乱数生成器と異なり，これらの列は一様ではなく，数値メソッドにおいて便利であることがある潜在的な構造を持つ．例えば，これらの列は通常多次元モンテカルロ(Monte Carlo)積分でより素早く収束する．

ハルトン(Halton)列

ハルトン列は単位区間で一様な準乱数を精製する．コードは任意の次元で動作するが，ここでは一次元空間でのみ動作するvan der Corput列について述べる．これは比較の意味で適切である．

ハルトン（またはvan der Corput）列は決定論的であるが，単位区間において相違が少ない．これは空間を一様に埋めるので，モンテカルロ積分など，適用分野によっては擬似乱数生成器より好ましいことがある．

指定のと基底．

基底における一次元ハルトン（またはvan der Corput）の値．

長さの列は次のように書ける．

van der Corput列は，基底の表現に数が与えられると，数を小数点の反対側にミラーし，列の値はとなる．

Wolfram言語では，IntegerDigitsを使って列が求められる．

Wolfram Language code:

VanDerCorput[base_][len_] := Table[
	With[{digits = Reverse@IntegerDigits[n, base]}, 
	Sum[2 ^ (-ii) * digits[[ii]], {ii, Length[digits]}]
	], {n, len}]

基底を2に設定すると，列の最初の1000要素が計算できる．

Wolfram Language code: x = VanDerCorput[2][1000];

結果をプロットすると，空間が一様に埋められることが分かる．

Wolfram Language code: ListPlot[%]

相違の少ない列の特性は，列の中の隣の要素は，前の要素がどこに位置するのかを知っているというものである．これはManipulateで見ることができる．

Wolfram Language code:

Manipulate[
	ListPlot[VanDerCorput[2][n]], {n, 10, 10000, 10}]

CUDA実装については，IntegerDigitsに相当する独自のものを実装しなければならないが，これは難しくはない．まず実装のソースコードをロードする．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

CUDAFunctionをロードする．

Wolfram Language code: CUDAHaltonSequence = CUDAFunctionLoad[{srcf}, "Halton", {{_Real, _, "Output"}, _Integer, _Integer}, 256]

出力用にメモリを割り当てる．乱数は1024個だけ生成する．

Wolfram Language code: mem = CUDAMemoryAllocate[Real, {1024}]

次元1で関数を実行する．

Wolfram Language code: CUDAHaltonSequence[mem, 1, 1024]

結果をプロットする．

Wolfram Language code: ListPlot[CUDAMemoryGet[mem]]

Sobol列

Sobol列も相違の少ない列である．これは次のCUDAファイルに実装されている．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

ブロック次元を{64,1}として関数をロードする．

Wolfram Language code:

sobolFun = CUDAFunctionLoad[{srcf}, "Sobol", {_Integer, _Integer, {_Integer, _, "Input"}, {_Real, _, "Output"}}, {64, 1}]

入力パラメータをロードする．Sobol列が必要とする方向ベクトルは前もって計算され，ファイルに保管される．

Wolfram Language code:

vectorCount = 100000;
dimension = 100;
directions = CUDAMemoryLoad[Flatten[Import[FileNameJoin[{$CUDALinkExampleDataPath, "soboldirection.txt"}], "Data"]]];
output = CUDAMemoryAllocate[Real, vectorCount * dimension];
gridDim = {64, dimension};

Sobol関数にパラメータを渡して実行する．

Wolfram Language code: sobolFun[vectorCount, dimension, directions, output, gridDim];

列の中の最初の10000個の値をプロットする．空間が一様に埋められることが分かる（準乱数生成器）．

Wolfram Language code: ListPlot[Take[CUDAMemoryGet[output], {1, 10000}]]

完了したら，メモリはアンロードしなければならない．

Wolfram Language code: CUDAMemoryUnload[output, directions]

Niederreiter列

Niederreiter列もよく使われる，相違の少ない列である．これは以下のCUDAファイルに実装されている．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

乱数ジェネレータで使うシードデータをロードする．

Wolfram Language code: hTable = Flatten[Import[FileNameJoin[{$CUDALinkExampleDataPath, "niederreitaer.txt"}], "Data"]];

CUDAFunctionをロードする．

Wolfram Language code:

niederreiterFun = CUDAFunctionLoad[{srcf}, "Niederreiter", {{_Integer, _, "Input"}, {_Real, _, "Output"}, _Integer, _Integer}, {128, 3}]

計算に必要なメモリを割り当て，適切なパラメータを設定する．次元3の列を生成し，解像度は31とする．

Wolfram Language code:

qrngDimensions = 3;
qrngResolution = 31;
numRNG = 100000;
dTable = CUDAMemoryLoad[hTable];
dOutput = CUDAMemoryAllocate[Real, qrngDimensions * numRNG];
gridDim = {128 * 128, 3};

関数を実行する．

Wolfram Language code: niederreiterFun[dTable, dOutput, 1234, numRNG, gridDim];

GPUから出力を取得する．

Wolfram Language code: hOutput = CUDAMemoryGet[dOutput];

列の最初の10,000個の値をプロットする．Sobol列のときのように空間が一様に埋められているのが分かる．

Wolfram Language code: ListPlot[Take[hOutput, 10000]]

完了したら，メモリはアンロードしなければならない．

Wolfram Language code: CUDAMemoryUnload[dTable, dOutput]

ハッシング乱数生成器

ハッシングを利用する乱数生成器は質のより低い乱数を生成するが，高速に生成する．これは多くの適用分野において十分である．

小型暗号化アルゴリズムハッシング

小型暗号化アルゴリズム(TEA)は非常に簡単なハッシングアルゴリズムで，以下のファイルに実装されている．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

CUDAFunctionをロードする．

Wolfram Language code: CUDATeaEncryption = CUDAFunctionLoad[{srcf}, "Tea", {{_Integer, _, "Output"}, _Integer}, 256]

出力用のメモリを割り当てる．

Wolfram Language code: mem = CUDAMemoryAllocate[Integer, {2048}]

CUDAFunction関数を呼び出す．

Wolfram Language code: CUDATeaEncryption[mem, 2048]

結果をプロットする．

Wolfram Language code: ListPlot[CUDAMemoryGet[mem]]

割り当てられたメモリを削除する．

Wolfram Language code: CUDAMemoryUnload[mem]

MD5ハッシング

乱数生成器にはその他の一般的なハッシング法が使える．以下はよく知られたハッシングアルゴリズムである，MD5アルゴリズムの実装である．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "md5_rand.cu"}]

ソースからCUDAFunctionをロードする．

Wolfram Language code: md5 = CUDAFunctionLoad[{srcf}, "gen_randMD5", {{"Integer32[4]", _, "Output"}, _Integer, _Integer}, 32]

出力メモリをロードする．

Wolfram Language code: mem = CUDAMemoryLoad[ConstantArray[0, {1024 * 4}], "Integer32[4]"]

CUDAFunctionを呼び出す．

Wolfram Language code: md5[mem, 1024, 7, 1024]

結果をプロットする．

Wolfram Language code: ListPlot[CUDAMemoryGet[mem]]

割り当てられたメモリを削除する．

Wolfram Language code: CUDAMemoryUnload[mem]

正規乱数

以下のアルゴリズムは正規分布の乱数を生成する．

逆累積正規分布

以下は正規分布の乱数を生成する方法を実装する．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

CUDAFunctionをロードする．

Wolfram Language code: CUDAInverseCND = CUDAFunctionLoad[{srcf}, "InverseCND", {{_Real, _, "InputOutput"}, _Integer, _Integer}, 256]

100,000個の乱数のためのメモリを割り当てる．

Wolfram Language code:

sampleCount = 100000;
mem = CUDAMemoryAllocate[Real, sampleCount];

CUDAFunctionを呼び出す．

Wolfram Language code: CUDAInverseCND[mem, sampleCount, 0]

メモリをWolfram言語内に取得する．

Wolfram Language code: samples = CUDAMemoryGet[mem];

Histogramを使って結果をプロットする．

Wolfram Language code: Histogram[samples, Automatic, "ProbabilityDensity"]

メモリをアンロードする．

Wolfram Language code: CUDAMemoryUnload[mem]

Box–Muller

Box–Mullerは一式の一様分布の乱数を与えられると，正規分布の数を生成する方法である．CUDAの実装は以下のファイルにある．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}]

CUDAFunctionをロードする．

Wolfram Language code: CUDABoxMuller = CUDAFunctionLoad[{srcf}, "BoxMuller", {{_Real, _, "InputOutput"}, _Integer}, 128]

入力引数を設定する．

Wolfram Language code:

MTRNGCount = 4096;
PATHN = 2 ^ 14;
NPerRNG = Ceiling[PATHN / MTRNGCount];
NPerRNG = If[EvenQ[NPerRNG], NPerRNG, NPerRNG + 1];
RANDN = MTRNGCount * NPerRNG;

メルセンヌツイスタ（2セクション前に定義したもの）を使って一様に分布した乱数を生成する．

Wolfram Language code: mem = MersenneTwister[RANDN]

一様な乱数のリストを正規分布乱数に変換する．

Wolfram Language code: CUDABoxMuller[mem, NPerRNG, MTRNGCount]

Histogramを使うとベル状の曲線が示される．

Wolfram Language code: Histogram[CUDAMemoryGet[mem], Automatic, "ProbabilityDensity"]

T割り当てられたメモリを削除する．

Wolfram Language code: CUDAMemoryUnload[mem];

乱数ジェネレータの適用例

乱数は多くの分野に適用できる．ここでは，モンテカルロ積分（と任意の関数を近似して）とブラウン(Brownian)運動のシミュレーションという2つの重要な適用例を示す．

π の近似

の値はモンテカルロ積分を使って近似できる．まず単位正方形で一様乱数を生成する．次に単位円の最初の四分円内の点の数を数える．それから結果を点の数で割る．これがとなる．

単位円内の点の数を数え，リダクションする実装である．

Wolfram Language code: srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "reduceInCircle.cu"}]

CUDAFunctionをロードする．

Wolfram Language code: CountInCircle = CUDAFunctionLoad[{srcf}, "countInCircle", {{_Real, _, "Input"}, {_Integer, _, "Output"}, _Integer}, 256]

1,000,000個の点を使う．

Wolfram Language code:

size = 1000000;
numThreads = 256;
numBlocks = Floor[(size + (numThreads * 2 - 1)) / (numThreads * 2)];

前述のメルセンヌツイスタアルゴリズムを使って乱数を生成する．

Wolfram Language code: randomNumbers = MersenneTwister[2 * size]

出力メモリを割り当てる．

Wolfram Language code: output = CUDAMemoryAllocate[Integer, numBlocks]

計算を行う．

Wolfram Language code: CountInCircle[randomNumbers, output, size]

出力メモリを取得する．

Wolfram Language code: CUDAMemoryGet[output]//Total

結果はWolfram言語のものと一致する．

Wolfram Language code: Select[Partition[CUDAMemoryGet[randomNumbers], 2][[ ;; size]], #[[1]] ^ 2 + #[[2]] ^ 2 ≤ 1&]//Length

計算速度は格段に速い．

Wolfram Language code: CountInCircle[randomNumbers, output, size];//AbsoluteTiming

Wolfram言語と比較する．

Wolfram Language code: Select[Partition[CUDAMemoryGet[randomNumbers], 2][[ ;; size]], #[[1]] ^ 2 + #[[2]] ^ 2 ≤ 1&];//AbsoluteTiming

モンテカルロ積分

モンテカルロ積分は多くの分野で使われる．ここではSqrt[x]を0から1まで積分する．

Wolfram Language code:

src = "
__device__ Real_t integratedFunction(Real_t x) {
	return sqrt(x);
}
__global__ void monteCarlo(Real_t * evals, Real_t * randList, mint length) {
   	int index = threadIdx.x + blockIdx.x*blockDim.x;
   	if (index < length) {
			evals[index] = integratedFunction(randList[index]);
   	}
}
";

関数をロードする．

Wolfram Language code: monteCarlo = CUDAFunctionLoad[src, "monteCarlo", {{_Real, "Output"}, {_Real, "Input"}, _Integer}, 256]

乱数にSobel準乱数生成器を使う．

Wolfram Language code:

srcf = FileNameJoin[{$CUDALinkPath, "SupportFiles", "random.cu"}];
sobolFun = CUDAFunctionLoad[{srcf}, "Sobol", {_Integer, _Integer, {_Integer, _, "Input"}, {_Real, _, "Output"}}, {64, 1}];
vectorCount = 100000;
dimension = 100;
directions = CUDAMemoryLoad[Flatten[Import[FileNameJoin[{$CUDALinkExampleDataPath, "soboldirection.txt"}], "Data"]]];
random = CUDAMemoryAllocate[Real, vectorCount * dimension];
gridDim = {64, dimension};
sobolFun[vectorCount, dimension, directions, random, gridDim];

長さに乱数生成器の数を使う．

Wolfram Language code: len = First["Dimensions" /. CUDAMemoryInformation[random]]

出力用にメモリを割り当てる．

Wolfram Language code: output = CUDAMemoryAllocate[Real, len]

関数を呼び出す．

Wolfram Language code: monteCarlo[output, random, len]

最初のいくつかの要素が適切であるかどうか確認する．

Wolfram Language code: CUDAMemoryGet[output][[ ;; 20]]

出力を足し合わせる．これはCUDAFoldを使って行う．

Wolfram Language code: CUDAFold[Plus, 0, output] / len

結果はNIntegrateのものと一致する．

Wolfram Language code: NIntegrate[Sqrt[x], {x, 0, 1}]

割り当てられたメモリをアンロードする．

Wolfram Language code: CUDAMemoryUnload[output, random]

ブラウン運動

シミュレーションのためのメモリを割り当てる．

Wolfram Language code: mem = CUDAMemoryLoad[RandomReal[{0.0, 1.0}, sampleCount]]

Wolfram Language code: CUDAInverseCND[mem, sampleCount, 1];

Wolfram Language code: samples = CUDAMemoryGet[mem];

擬似乱数列の値は正規分布している．

Wolfram Language code: brownianMotion = FoldList[Plus, 0, samples];

Wolfram Language code: ListPlot[brownianMotion]

コードの生成

CUDALink はWolfram言語に統合されているので，CUDAカーネルコードの生成にSymbolicC等のWolfram言語機能を利用することができる．CUDALink をまだインポートしていない場合はインポートする．

Wolfram Language code: Needs["CUDALink`"]

個の例ではSymbolicCパッケージが必要である．

Wolfram Language code: Needs["SymbolicC`"]

いくつかの一般的なWolfram言語構文をSymbolicC表現に変換して定義する．

Wolfram Language code:

ClearAll[toSymbolicC]
ClearAll[x]
ClearAll[xx]
SetAttributes[toSymbolicC, {HoldAll}]
toSymbolicC[x_List] := toSymbolicC /@ x
toSymbolicC[Times[-1, x_]] := "-" <> ToCCodeString[toSymbolicC[x]]
toSymbolicC[(op : (Plus | Times))[args___]] := COperator[op, toSymbolicC[{args}]]
toSymbolicC[(op : (Minus | BitNot | Not | Decrement | Increment | PreDecrement | PreIncrement))[x_]] := COperator[op, toSymbolicC[x]]
toSymbolicC[(op : (Mod | Divide | Subtract | BitShiftRight | BitShiftLeft))[x_, y_]] := COperator[op, {toSymbolicC[x], toSymbolicC[y]}]
toSymbolicC[(op : (ArcCos | ArcSin | Ceiling | Cos | Cosh | Exp | Abs | Floor | Sin | Sinh | Sqrt | Tan | Tanh | Log))[x_]] := CStandardMathOperator[op, toSymbolicC[x]]
toSymbolicC[Power[x_, r : Rational[_, _]]] := CStandardMathOperator[Power, {toSymbolicC[x], toSymbolicC[r]}]
toSymbolicC[Power[x_, 2]] := COperator[Times, {toSymbolicC[x], toSymbolicC[x]}]
toSymbolicC[Power[x_, y_]] := CStandardMathOperator[Power, {toSymbolicC[x], toSymbolicC[y]}]
toSymbolicC[CompoundExpression[stmts__]] := toSymbolicC /@ stmts
toSymbolicC[If[cond_, trueStmt_]] := CIf[toSymbolicC[cond], toSymbolicC[trueStmt]]
toSymbolicC[If[cond_, trueStmt_, falseStmt_]] := CIf[toSymbolicC[cond], toSymbolicC[trueStmt], toSymbolicC[falseStmt]]
toSymbolicC[x_Rational] := N[x]
toSymbolicC[x_] := x

検証するためにWolfram言語文を渡すと，SymbolicC出力が得られる．

Wolfram Language code: toSymbolicC[Sin[x] ^ 3 + x ^ 8 + 3]

C文字列に変換するにはToCCodeStringメソッドを使う．

Wolfram Language code: ToCCodeString[%]

以上により，Wolfram言語関数（純粋であってもそうでなくても）を取り，適切なCUDAカーネルソースを生成する関数が書けるようになる．

Wolfram Language code:

SetAttributes[CUDAMapSource, {HoldAll}];
ClearAll[CUDAMapSource];
CUDAMapSource[f_] := ToCCodeString[With[{fun  = f[xx] /. xx -> CArray["lst", "index"]}, 
	SymbolicCUDAFunction["map", {{CPointerType[{"Real_t"}], "lst"}, {"mint", "length"}}, 
	CBlock[{
	SymbolicCUDADeclareIndexBlock[1], 
	CIf[COperator[Less, {"index", "length"}], 
	CAssign[CArray["lst", "index"], toSymbolicC[fun]]
	]
	}]
	]]]

純粋関数をCUDAMapSourceに渡すと，カーネルコードが返される．

Wolfram Language code: CUDAMapSource[# + Sin[#]&]

Wolfram言語関数と入力リストを与えられると，CUDAカーネルコードを生成し，コードをCUDAFunctionとしてロードし，CUDAFunctionを実行し，結果を返す関数を定義する．

Wolfram Language code:

SetAttributes[myCUDAMap, HoldFirst];
myCUDAMap[fun_, input_List] := 
	Module[{len = Length[input], res, oclFun, output}, 
	oclFun = CUDAFunctionLoad[CUDAMapSource[fun], "map", {{_Real, "InputOutput"}, _Integer}, 256];
	First[oclFun[input, len]]
	]

入力リストの各要素に2を足す純粋関数でmyCUDAMapを検証する．

Wolfram Language code: myCUDAMap[# + 2&, ConstantArray[1.0, 100]]

toSymbolicCが変換したあらゆる構文がmyCUDAMapによってサポートされる．以下では各要素を2乗する．

Wolfram Language code: myCUDAMap[# ^ 2&, Range[100]]

モンテカルロ積分を行う．

Wolfram Language code: CUDAFold[Plus, 0, myCUDAMap[# ^ 2&, RandomReal[1, 1000000]]] / 1000000

関数は定義してmyCUDAMapに渡すことができる．次は色反転関数を定義する．

Wolfram Language code: colorNegate[x_] := 1.0 - x

色反転関数を呼び出す．

Wolfram Language code: myCUDAMap[colorNegate, Range[10]]

以上は簡単な例であるが，より複雑なものの基盤として使える．

Top

適用例

画像処理

画像の二値化

ボックスフィルタ

画像の調整

キャニー(Canny)エッジ検出

線形代数とリスト処理

行列の転置

行列とベクトルの乗算

行列と行列の乗算

ドット積

凸包

乱数生成

擬似乱数生成器

Park–Miller

メルセンヌツイスタ

準乱数生成器

ハルトン(Halton)列

Sobol列

Niederreiter列

ハッシング乱数生成器

小型暗号化アルゴリズムハッシング

MD5ハッシング

正規乱数

逆累積正規分布

Box–Muller

乱数ジェネレータの適用例

π の近似

モンテカルロ積分

ブラウン運動

コードの生成

関連するガイド

関連するテクニカルノート