CUDAプログラミング
CUDAはNVIDIAが開発した，GPU（グラフィカルプロセッシングユニット）をプログラムするための，Cに似た一般的なプログラム言語である．CUDALink はGPUをプログラムするのに必要な多くのステップを省く使いやすいインターフェースを提供する．コンパイル，リンキング，データの転送等はすべてWolfram言語のCUDALink が行う．これによりユーザはインターフェースやコードを書くのではなく，アルゴリズムを書くことに集中できるようになる．
このセクションではWolfram言語でCUDAプログラミングを開始する方法について述べる．
| CUDAFunctionLoad | CUDA関数をWolfram言語にロードする |
このドキュメントでは，GPUアーキテクチャと，CUDAカーネルの書き方について述べる．最後には CUDALink を使って書かれた多くの適用例を示す．
はじめに
Common Unified Device Architecture (CUDA)はNVIDIAによって2007年後半に開発されたGPUをより一般的にする方法である．CUDALink はGPUプログラミングを簡単にしGPUのプログラミングは何年も前存在したが，プログラミングの難しさからその利用は限られていた．CUDALink はGPUプログラミングを簡単にし，より多く採用されることを狙いとしている．
Wolfram言語を使うときには，多くのステップについて考えなければならない．Wolfram言語では，CUDAをカーネルを書くだけである．これはNVIDIAアーキテクチャスタックの全レベルを使って行われる．
CUDAアーキテクチャ
CUDAのプログラミングはデータ並列モデルに基づいている．高レベルの観点から見ると，問題はまず計算のために何百から何千ものスレッドに分割される．以下の計算を行うとする：
OutputData = Table[fun[InputData[[i, j]]], {i, 10000}, {j, 10000}]
CUDAのプログラミングパラダイムでは，この計算は以下に等しい：
CUDALaunch[fun, InputData, OutputData, {10000, 10000}]
は計算関数である．上式は10000×10000スレッドを開始し，指標を各スレッドに渡し，InputDataに関数を適用し，OutputDataに結果を置く．CUDALink のCUDALaunchに相当するものはCUDAFunctionである．
CUDAが何千ものスレッドを開始できる理由はハードウェアアーキテクチャにある．次のセクションでは，それについてと，スレッドが実行するのにどのように分割されるかについて述べる．
CUDAのグリッドとブロック
CUDAプログラミングはメッセージの受け渡しやスレッドベースの並列計算モデルとは異なり，問題を一，二，あるいは三次元のグリッドにマップする．以下では，一般的な二次元CUDAスレッドの構成を示す．
各グリッドには複数のブロックが含まれ，各ブロックには複数のスレッドが含まれる．上の図におけるグリッド，ブロック，スレッドを以下に示す．
一，二，あるいは三次元のどのスレッド構成を選ぶかは問題による．例えば画像処理の場合は，次の図のようにして画像をスレッド上にマップし，各画素に関数を適用する．
一方，一次元セルオートマトンは一次元問題なので，次のようにマップできる：
CUDAプログラムサイクル
CUDAプログラミングの概要は，多く（通常数千）のスレッドを開始してデータをコピーし，GPUの実行が終了するのを待ち（あるいは待機中にCPU計算を行い），最後に結果をデバイスからホストにコピーするというものである．
上図はCUDAプログラムの一般的なサイクルの詳細を示している．
1. GPUのメモリを割り当てる．GPUメモリはCPUメモリとは別のもので，プログラマは割当てコピーを管理しなければならない．
3. スレッドの構成を設定する．問題に適したブロックとグリッドの次元を選ぶ．
5. CUDAスレッドを同期させると，デバイスがGPUメモリ上でその後の操作を行う前にすべてのタスクを確実に完了させるようにできる．
6. スレッドが完了したら，メモリはGPUからCPUにコピーして返される．
Wolfram言語を使うときには，多くのステップについて考えなければならない．Wolfram言語では，CUDAをカーネルを書くだけである．
メモリの階層
CUDAのメモリは複数のレベルに分割され，それぞれが独自の利点と制限を持つ．次の図はCUDAで利用できるすべてのタイプのメモリを示す．
大域メモリ
GPU上で利用できる，最も豊富だが最も遅いメモリ．このメモリは128，256，512MBのパッケージングがある．すべてのスレッドは大域メモリの要素にアクセスできるが，パフォーマンス上の理由からこれらのアクセスは最小限に抑えられ，さらに制約がある．
大域メモリにおけるパフォーマンスの制約は，最近のハードウェアでは緩められ，今後さらに緩められる予定である．一般に，大域メモリは複数回アクセスされると，パフォーマンスは劣化する．
テクスチャメモリ
テクスチャメモリは大域メモリと同じ場所にあるが，読取り専用である．テクスチャメモリは大域メモリに見られるパフォーマンスの低下はないが，サポートしている型はchar，int，floatのみである．
コンスタントメモリ
グリッド内のどのスレッドからもアクセスできる高速なコンスタントメモリ．メモリはキャッシュされるが，全部で64KBに制限される．
共有メモリ
特定のブロックに局所的な高速メモリ．現在のハードウェア上では共有メモリ量はブロックごとに16KBに制限されている．
局所メモリ
局所メモリは各スレッドに局所的であるが，音波医らが変数をレジスタに置かない限り，大域メモリ内に存在する．パフォーマンスについては，メモリアクセスに時間がかかるため，一般に局所メモリは最小限に抑えるべきであると考えられているが，このメモリについては大域メモリのような問題は生じない．
計算能力
使用中のシステムの計算能力についての情報はCUDAInformationで取り出すことができる．
システム上の全CUDAデバイスと"Compute Capabilities"を表示する．
複数のCUDAデバイス
システムハードウェアがサポートしている場合はCUDAでは計算を行うデバイスを選択することができる．デフォルトでは最も速いものが選ばれるが，ユーザはこれを変更することができる．一旦デバイスを設定したら（自動的に選ばれたものでもユーザが指定したものでも），カーネルセッション中はデバイスを変更することはできない．
CUDAカーネルを書く
CUDAカーネルは何度も呼ばれる原子関数である．これらは通常プログラムのForループの中の数行である．次の例は2つのベクトルを足し合せる．
1つ目のカーネル
CUDAカーネルは入力リストの各要素について計算を行う小さいコードである．この例の1つ目のカーネルは各要素に2を足す．
__global__ void addTwo_kernel(mint * arry, mint len) {
int index = threadIdx.x + blockIdx * blockDim.x;
if (index >= len) return;
arry[index] += 2;
}
は関数をGPU上で実行するようにコンパイラに指示する関数修飾子である．
関数はCから呼び出すことができる．もう一つの関数修飾子は
で，他の
関数や
関数から呼び出せるがCからは呼び出せない関数を表す．
CUDAカーネルは無効の出力を持たなければならない．そのためその結果を得るためにはポインタ入力を私，それらを上書きしなければならない．この例では
を渡し，要素を上書きする（
は入力・出力パラメータと考えることができる）．
これは指標を取得する．CUDAカーネルは起動グリッドとブロックサイズによって設定される次の変数を提供する．
これらのパラメータはカーネル起動パラメータ（ブロックとグリッドの次元）に基づき，CUDAによって自動的に設定される．低い次元の計算を起動するとき，より高い次元は自動的に1に設定される．そのため1Dのグリッドを起動すると，
と
は1に設定される．
スレッドは複数のブロック次元で起動されるので，入力の次元がブロック次元の倍数でない場合は境界を上書きしないようにしなければならない．これはそれを行っている．
これはリストの各要素に適用される関数である．この場合は入力に2を足す．
2つ目のカーネル
__device__ float f(float x) {
return tanf(x);
}
__global__ void secondKernel(float * diff, float h, mint listSize) {
mint index = threadIdx.x + blockIdx.x * blockDim.x;
float f_n = f(((float) index) / h);
float f_n1 = f((index + 1.0) / h);
if( index < listSize) {
diff[index] = (f_n1 - f_n) / h;
}
}
これはデバイス関数である．
関数は
関数から呼び出せるが，Wolfram言語から直接呼び出すことはできない．
浮動小数点入力にTanを適用する．
リストをdiffを返し，2つのスカラ入力h，listSizeを取る関数．
index/hで関数を評価する．indexとhは整数なので，切り捨てられないようにするために一つを浮動小数点数にしなければならない．デバイス関数をfと呼ぶ．
次のセクションでは，このCUDAプログラムをWolfram言語にロードする方法を示す．
CUDAカーネルをWolfram言語にロードする
CUDAFunctionをロードする．第1引数にカーネルコード，第2引数に実行するカーネル名("secondKernel")，第3引数に関数パラメータ，第4引数にブロックサイズを渡す．結果はsecondKernelに保管される．
CUDAFunctionは一旦ロードされたら他のWolfram言語関数と同様用に使える．実行する前に結果を保管するためのバッファが必要である．サイズ1000.のバッファを作成する．
直列的なCからCUDAプログラミングへ
以下に直列CPUバージョンからCUDAバージョンへのプログラムの工程の詳細を示す．
プログラムは半径3の移動平均を実装し，隣接画素に基づく配列の平均値を計算する．
実装は直列的なものから始めるのが一般的である．これにより，問題を正しく読み取って並列化の可能な方法を見付けることができるようになる．直列実装は参照のための実装として使えるだけでなく，CUDAがタスクにフィットするかどうかを確定するためにも使える．
void xAverage_cpu (mint * in, mint * out, mint width, mint height, mint pitch) {
int xIndex, yIndex, idx;
mint left, curr, right;
for (yIndex = 0; yIndex < height; yIndex++) {
for (xIndex = 0; xIndex < width; xIndex++) {
idx = x + y*pitch;
left = idx > 0 ? in[idx - 1] : in[yIndex*pitch];
curr = in[idx];
right = idx < width - 1 ? in[idx + 1] : in[width + yIndex*pitch];
out[x + y*pitch] = (left + curr + right)/3;
}
}
}
最適化の可能な範囲を見付けるのはよい考えである．この場合は，要素が右端であるかどうかを確認するifと，左端であるかどうかを確認するifの2つのif文を削除する．
void xAverage_cpu(mint * in, mint * out, mint width, mint height, mint pitch) {
int xIndex, yIndex, idx;
for (yIndex = 0; yIndex < height; yIndex++) {
out[yIndex*pitch] = (2*in[yIndex*pitch] + in[yIndex*pitch + 1])/3;
for (xIndex = 1; xIndex < width-1; xIndex++) {
idx = x + y*pitch;
out[idx] = (in[idx-1] + in[idx] + in[idx+1])/3;
}
out[width - 1 + yIndex*pitch] = (in[width - 2 + yIndex*pitch] +
2*in[width - 1 + yIndex*pitch])/3;
}
}
もう一つ考えるべきことは隣接画素に依存するコードを避けることである．
void xAverage_cpu(mint * in, mint * out, mint width, mint height, mint pitch) {
mint accum;
mint left, curr, right;
for (yIndex = 0; yIndex < height; yIndex++) {
accum = 2*in[yIndex*pitch] + in[yIndex*pitch + 1]
out[yIndex*pitch] = accum/3;
accum -= in[yIndex*pitch];
for (xIndex = 0; xIndex < width-1; xIndex++) {
idx = x + y*pitch;
accum += in[index+1];
out[idx] = accum/3;
accum -= in[idx-1];
}
accum += in[width - 2 + yIndex*pitch];
out[width - 1 + yIndex*pitch] = accum/3;
}
}
__global__ void xAverage_kernel(mint * in, mint * out, mint width, mint height, mint pitch) {
int xIndex = blockDim.x*blockIdx.x + threadIdx.x;
int yIndex = blockDim.y*blockIdx.y + threadIdx.y;
int index = xIndex + yIndex*pitch;
if (xIndex >= width || yIndex >= height) return;
if (xIndex > 0) {
out[index] = (2*in[index] + in[index + 1])/3;
} else if (xIndex < width) {
out[index] = (in[index - 1] + in[index] + in[index + 1])/3;
} else {
out[index] = (in[index - 1] + 2*in[index])/3;
}
}
大域メモリコピーが同じデータについて3回行われることに注目されたい．これは共有メモリに一時的にデータを置くことで避けられる．
__global__ void xAverage_improved_kernel(mint * in, mint * out, mint width, mint height, mint pitch) {
extern __shared__ smem[];
int tx = threadIdx.x, dx = blockDim.x;
int xIndex = dx*blockIdx.x + tx;
int yIndex = blockDim.y*blockIdx.y + threadIdx.y;
int index = xIndex + yIndex*pitch;
if (yIndex >= height) return;
smem[threadIdx.x+1] = index < width ? in[index] : in[yIndex*pitch + width - 1];
if (tx == 0) {
smem[0] = index > 0 ? in[index-1] : in[yIndex * pitch];
smem[dx+1] = index+bx < width ? in[index+dx] : in[yIndex*pitch + width - 1];
}
__syncthreads();
if (xIndex < width)
out[index] = (smem[tx] + smem[tx+1] + smem[tx+2])/3;
}
このCUDAプログラムはCUDAFunctionLoadを使って保存，ロードできる．
CUDAのためのコンパイル
CUDALink はCUDAコードのためのポータブルなコンパイルメカニズムを提供する．これはCCompilersでシステム上で見付かったNVCCコンパイラを登録し，DLLと実行ファイルをコンパイルするNVCCCompilerによって行われる．
CUDALink にNVCCCompilerをロードする前に，システム上で利用できるコンパイラを調べる．
上記出力は，NVCCCompilerをロードしていない場合のものである．これが当てはまる場合は CUDALink アプリケーションをロードする．
これでNVCCCompilerが使えるようになった．
PTXファイルの生成
PTXは別のCUDAカードで実行できるCUDAのバイトコードである．バイトコードはその後オンザフライで（JITメカニズムを使って）コンパイルされる．
実行ファイルの作成時に"CreatePTX"->Trueと設定すると，コンパイル後にPTXになる．
CUBINファイルの生成
CUBINファイルは特定のCUDAアーキテクチャのためにコンパイルされるCUDAバイナリファイルである．例えばCUDAFunctionLoadは入力コードをコンパイルしてCUBINファイルにする．
実行ファイルの作成時に"CreateCUBIN"->Trueと設定すると，コンパイル後にCUBINになる．
ASCII命令集合を含むPTXファイルとは違い，CUBINファイルはELFバイナリである．
デバッグ情報とレジスタ情報の表示
CUDAコードをデバッグしているときまたは最適化しているときは，コンパイラ出力が見られると便利である．例えばCUDA関数が使用するレジスタ数を求めるのは定番である．
"Debug"->Trueと"CompileOptions"->"--ptxas-options=-v"を設定すると，デバッグ情報を表示し，バイナリを生成するときに冗長にするようにコンパイラに指示する．
CUDALink 用のカーネルを書く
CUDAプログラミングに詳しいユーザは，CUDALink 用にプログラミングするときに以下の方法を考えるかもしれない．この方法はCUDAコードをよりポータブルで，正しく，デバッグしやすく，効率的にする．
実数型を使う
多くのCUDAカードで倍精度はサポートされていないので，CUDA Cプログラマは浮動小数点数を表すのにfloatを使わなければならない．CUDALink では，CUDAカードで使える最高精度を使う浮動小数点数を定義するのに
型が使える．
__device __ Real_t f(Real_t x) {
return tanf(x);
}
__global __ void secondKernel(Real_t * diff, Real_t h, mint listSize) {
int index = threadIdx.x + blockIdx.x * blockDim.x;
Real_t f_n = f(index) / h);
Real_t f_n1 = f((index + 1.0) / h);
if( index < listSize) {
diff[index] = (f_n1 - f_n) / h;
}
}
__device __ float f(float x) {
return tanf(x);
}
__global __ void secondKernel(float * diff, float h, mint listSize) {
int index = threadIdx.x + blockIdx.x * blockDim.x;
float f_n = f(((double) index) / h);
float f_n1 = f((index + 1.0) / h);
if( index < listSize) {
diff[index] = (f_n1 - f_n) / h;
}
}
これにより能力を最大限に引き出し，利用できる最高精度が使えるようになる．
最適化
ループパラメータを定数にすると，コンパイラはループを最適化できることがある．次の例は，コンパイラがさらにコード最適化を行うことができるため，チャンネルを渡すものよりも少し最適化される．
__global __ void colorNegate(unsigned char * img, mint width, mint height) {
int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
int yIndex = threadIdx.y + blockIdx.y * blockDim.y;
int index = (xIndex + yIndex*width)*CHANNELS;
if (xIndex < width && yIndex < height) {
#pragma unroll
for (mint ii = 0; ii < CHANNELS; ii++)
img[index + ii] = 255 - img[index + ii];
}
}
コードの中でを定数にするのではなく，CUDAFunctionLoadを使ってそれを"Defines"として渡すことができる．
CUDAFunctionLoad[...,"Defines"{"CHANNELS"3}]
オフラインコンパイル
ソースコードをCUDAFunctionLoadに渡すと，コードは裏でコンパイルされる．毎回のパフォーマンスヒットを避けるために，ユーザはNVCCCompilerでCreateExecutableを使ってCUDAコードをコンパイルしてバイナリにし，CUDAバイナリからCUDAFunctionをロードすることができる．
OpenCLをCUDAにポートする
CUDALink と OpenCLLink は両方ともCUDAとOpenCLの潜在的なメモリ管理をすべて扱うので，ユーザはカーネルコードのみに集中ればよい．このセクションではCUDAからOpenCLへの変換について述べる．ユーザは関数の開発に OpenCLLink を使ったことが想定されている．
OpenCLプログラムをCUDAにポートするときに必要なカーネルコードの一般的な変化を以下の表にまとめる．
OpenCL関数を開始するには，次の式を実行しなければならない．
多くの場合，Wolfram言語におけるOpenCLからCUDAへ４のポートの際に必要な変更は上記のものだけである．OpenCLまたはCUDAのプログラミングガイドに他の一対一の関数変換がいくつか記載されている．
次はOpenCL用に行われた，Wolfram言語関数そのままであるが，OpenCLFunctionLoadのロードをCUDAFunctionLoadで置き換える．
上で示したメカニズムを検索・置換する方法の代替として，カーネルコードを記号的に生成し，それをWolfram言語で操作するというものがある．以下のセクションでその方法を示す．
記号的なコードの生成
Wolfram言語は，Cコードの階層的な表示をWolfram言語独自の言語として見ることを可能にするSymbolicCをサポートする．これにより，Wolfram言語はCコードの作成，操作，最適化に適したものとなる．
この機能に加え，CUDAカーネルコードをいくつかの別々の目的のために生成することができるため，よりよいポータブル性，より少ないプラットフォーム依存性，よりよいコード最適化等が可能となる．
SymbolicCを使うことの利点としては，Cコードが確実に文法エラーを含んでいないようにできることや，SymbolicCを使ってCUDAプログラムを生成し，CUDAFunctionLoadを使ってそれらを実行できること，そしてCUDAからOpenCLへの（またはその反対の）ポートが簡単になるということが挙げられる．
このセクションは CUDALink の記号コード生成機能を使って「RGB」から「HSV」，および「RGB」から「HSB」への色変換を書く．この変換は画像・映像処理において非常に一般的である．
このチュートリアルでは，まず最初に CUDALink アプリケーションをロードしなければならない．
はじめに
CUDALink には，SymbolicCを使ってカーネルを書く手助けとなる記号ツールが含まれている．
| SymbolicCUDAFunction | CUDA関数の記号表現 |
| SymbolicCUDABlockIndex | CUDAのブロックインデックス呼出しの記号表現 |
| SymbolicCUDAThreadIndex | CUDAのスレッドインデックス呼出しの記号表現 |
| SymbolicCUDABlockDimension | CUDAのブロック次元呼出しの記号表現 |
| SymbolicCUDACalculateKernelIndex | CUDAのインデックス計算の記号表現 |
| SymbolicCUDADeclareIndexBlock | CUDAのインデックス宣言の記号表現 |
ToCCodeStringで生成されたCコードを表示する．
コードを記号的に書くもう一つの理由は，"CUDA"を"OpenCL"に変更するだけで，上記のカーネルをOpenCL形式に書き換えることができるというものである．
ToCCodeStringを使って再びCコードを表示する．
簡単なカーネル
CUDALink には上述のようなよくあるGPUカーネルプログラミングタスクを行うユーティリティ関数がある．それらを使って最初のカーネル関数を書いてみる．
最初の関数は整数配列を取り，その全要素を0に設定する．まず関数のプロトタイプを定義する必要がある．関数は"firstKernel"で，入力リストとリストサイズを取る．以下のようにユーティリティ関数を使う．
ここではこの例を1徐々に構築するため，また，上記を毎回書くことを避けるため，この例のためのヘルパー関数を定義する．
インデックスの取得
SymbolicCUDADeclareIndexBlockを使い，スレッドインデックス，ブロックインデックス，ブロック次元に基づいてリスト要素のインデックスを求める．
本体の宣言
要素の数より多いスレッドを起動することができるので，関数本体をif文の中に置く必要がある．これを行わないと，プログラムが予測不可能となる可能性がある．
これでカーネル本体が定義できる．本体は簡単なもので，リスト内の各要素を0に設定する．
最初のカーネルの検証
CUDALink をまだインポートしていない場合はインポートする．
高度なカーネル
このセクションではRGBからHSV，RGBからHSBへの変換器を書く．これは画像・映像処理において一般的な変換である．
ヘルパー関数
SymbolicCはCの構文のみを含んでいる．コードをよりWolfram言語風にするために，以下を定義する．
これを使ってFileNameJoin[{$CUDALinkPath,"SupportFiles","colorConvert.cu"}]のコードを再実装する．
色の変換
問題を部分要素に分割する．まずRGBからHueに変換する関数を定義する．この関数はRGBからHSV，RGBからHSBに共通である．
カーネル関数を生成するためにWolfram言語の関数型機能を利用する．ヘルパー関数を定義する．
上記関数を使い，RGBからHSBへの操作が単精度演算を使うように定義する．
プログラムの実行
上記コードはCUDAFunctionLoadを使って実行できる．まず単精度でソースを生成する．
GPU上でコードを実行する前に，必要なパラメータを定義する必要がある．以下は入出力画像を定義する．HSVとHSBは実数値なので，出力もReal型であることに注意．