適用例

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

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

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

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

画像処理

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

画像の二値化

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

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

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

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

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

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

出力画像を表示する．

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

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

ボックスフィルタ

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

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

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

半径は5に設定される．

関数を呼び出す．

画像を得る．

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

画像の調整

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

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

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

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

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

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

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

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

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

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

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

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

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

線形代数とリスト処理

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

行列の転置

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

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

行列を設定する．

行列Aを転置する．

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

行列とベクトルの乗算

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

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

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

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

行列と行列の乗算

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

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

計算を行う．

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

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

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

ドット積

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

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

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

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

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

凸包

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

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

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

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

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

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

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

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

凸包を計算する．

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

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

凸包の点を計算する．

凸包を可視化する．

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

乱数生成

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

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

擬似乱数生成器

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

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

Park–Miller

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

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

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

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

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

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

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

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

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

CUDAFunctionを呼び出す．

結果はランダムである．

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

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

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

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

メルセンヌツイスタ

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

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

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

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

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

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

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

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

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

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

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

準乱数生成器

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

ハルトン(Halton)列

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

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

指定の と基底 ．

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

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

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

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

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

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

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

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

CUDAFunctionをロードする．

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

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

結果をプロットする．

Sobol列

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

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

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

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

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

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

Niederreiter列

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

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

CUDAFunctionをロードする．

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

関数を実行する．

GPUから出力を取得する．

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

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

ハッシング乱数生成器

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

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

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

CUDAFunctionをロードする．

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

CUDAFunction関数を呼び出す．

結果をプロットする．

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

MD5ハッシング

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

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

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

CUDAFunctionを呼び出す．

結果をプロットする．

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

正規乱数

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

逆累積正規分布

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

CUDAFunctionをロードする．

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

CUDAFunctionを呼び出す．

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

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

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

Box–Muller

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

CUDAFunctionをロードする．

入力引数を設定する．

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

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

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

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

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

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

π の近似

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

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

CUDAFunctionをロードする．

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

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

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

計算を行う．

出力メモリを取得する．

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

計算速度は格段に速い．

Wolfram言語と比較する．

モンテカルロ積分

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

関数をロードする．

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

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

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

関数を呼び出す．

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

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

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

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

ブラウン運動

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

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

コードの生成

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

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

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

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

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

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

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

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

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

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

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

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

色反転関数を呼び出す．

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