CUDA&#12503;&#12525;&#12464;&#12521;&#12511;&#12531;&#12464;

システム上の全CUDAデバイスとを表示する．

Out[2]=

複数の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カーネルは起動グリッドとブロックサイズによって設定される次の変数を提供する．

— 現在のスレッドの指標．スレッド指標は0との間である．

— 現在のブロックの指標．ブロック指標は0との間である．

— ブロックサイズ次元．

— グリッドサイズ次元．

これらのパラメータはカーネル起動パラメータ（ブロックとグリッドの次元）に基づき，CUDAによって自動的に設定される．低い次元の計算を起動するとき，より高い次元は自動的に1に設定される．そのため1Dのグリッドを起動すると，とは1に設定される．

1Dの最も一般的な場合では，大域オフセットを求めるのにはを使い，2Dの大域的位置を求めるのにを使う．

スレッドは複数のブロック次元で起動されるので，入力の次元がブロック次元の倍数でない場合は境界を上書きしないようにしなければならない．これはそれを行っている．

これはリストの各要素に適用される関数である．この場合は入力に2を足す．

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

これはデバイス関数である．関数は関数から呼び出せるが，Mathematica から直接呼び出すことはできない．

浮動小数点入力にTanを適用する．

リストをを返し，2つのスカラ入力，を取る関数．

大域スレッドオフセットを得る．

で関数を評価する．とは整数なので，切り捨てられないようにするために一つを浮動小数点数にしなければならない．デバイス関数をと呼ぶ．

で関数を評価する．

出力バッファを上書きしないようにする．

前進差分を計算し，結果をに保管する．

次のセクションでは，このCUDAプログラムを Mathematica にロードする方法を示す．

CUDAカーネルを Mathematica にロードする

まだCUDALink をロードしていない場合はロードする．

全セクションのコードを文字列に置く．

CUDAFunctionをロードする．第1引数にカーネルコード，第2引数に実行するカーネル名()，第3引数に関数パラメータ，第4引数にブロックサイズを渡す．結果はに保管される．

In[3]:=

Out[3]=

裏では以下のことが行われている．

システムのCUDAが利用できるかどうか確認される．

CUDAがコンパイルされて選択されたGPUに最適化されたバイナリファイルにされる．

コンパイルされたコードはさらにコンパイルされるのを防ぐためにキャッシュされる．

コンパイルされたコードにカーネルが存在するかどうか確認される．

CUDAFunctionは一旦ロードされたら他の Mathematica 関数と同様用に使える．実行する前に結果を保管するためのバッファが必要である．サイズ1000.のバッファを作成する．

In[4]:=

Out[4]=

を実行する．

In[5]:=

Out[5]=

出力バッファの最初の30要素を取得する．

Out[6]=

出力バッファをプロットする．

In[7]:=

Out[7]=

バッファはユーザが手動で解放しなければならない．

In[28]:=

直列的な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;
    }
}

最後に最も最適化されたCPU実装に基づくプログラムを書く．

__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をロードする前に，システム上で利用できるコンパイラを調べる．

Out[1]=

上記出力は，NVCCCompilerをロードしていない場合のものである．これが当てはまる場合は CUDALink アプリケーションをロードする．

もう一度コンパイラを確認する．

In[3]:=

Out[3]=

これでNVCCCompilerが使えるようになった．

テストファイルの例をロードする．

In[4]:=

Out[4]=

ファイルの内容が見られる．

In[5]:=

コードをコンパイルしてライブラリにする．

Out[6]=

ライブラリは上記の場所に置かれる．

PTXファイルの生成

PTXは別のCUDAカードで実行できるCUDAのバイトコードである．バイトコードはその後オンザフライで（JITメカニズムを使って）コンパイルされる．

簡単なCUDAカーネルを定義する．

実行ファイルの作成時に"CreatePTX"->Trueと設定すると，コンパイル後にPTXになる．

Out[2]=

コンパイルされたコードの最初の500文字を表示する．

In[3]:=

Out[3]=

CUBINファイルの生成

CUBINファイルは特定のCUDAアーキテクチャのためにコンパイルされるCUDAバイナリファイルである．例えばCUDAFunctionLoadは入力コードをコンパイルしてCUBINファイルにする．

簡単なCUDA関数を定義する．

In[4]:=

実行ファイルの作成時に"CreateCUBIN"->Trueと設定すると，コンパイル後にCUBINになる．

In[5]:=

Out[5]=

コンパイルされたコードの最初の100文字を表示する．

Out[6]=

ASCII命令集合を含むPTXファイルとは違い，CUBINファイルはELFバイナリである．

デバッグ情報とレジスタ情報の表示

CUDAコードをデバッグしているときまたは最適化しているときは，コンパイラ出力が見られると便利である．例えばCUDA関数が使用するレジスタ数を求めるのは定番である．

簡単なCUDA関数を定義する．

In[7]:=

"Debug"->Trueとを設定すると，デバッグ情報を表示し，バイナリを生成するときに冗長にするようにコンパイラに指示する．

In[8]:=

Out[8]=

この簡単なカーネルは上で示したように2つのレジスタを使う．

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を使ってそれをとして渡すことができる．

CUDAFunctionLoad[...,"Defines"→{"CHANNELS"→3}]

Cコードを生成してライブラリとしてロードする

プロトタイプ化が完了し，コードが頻繁に使われるなら，コードのCコード表現はCUDACodeStringGenerateを使ってエキスポートすることができる．結果のCコードはコンパイルしてCUDAFunctionLoadにロードできる．これによりCUDAFunctionLoadはエラーチェックの頻度が少ない別の方法で実行されることになる．

オフラインコンパイル

ソースコードをCUDAFunctionLoadに渡すと，コードは裏でコンパイルされる．毎回のパフォーマンスヒットを避けるために，ユーザはNVCCCompilerでCreateExecutableを使ってCUDAコードをコンパイルしてバイナリにし，CUDAバイナリからCUDAFunctionをロードすることができる．

OpenCLをCUDAにポートする

CUDALink と OpenCLLink は両方ともCUDAとOpenCLの潜在的なメモリ管理をすべて扱うので，ユーザはカーネルコードのみに集中ればよい．このセクションではCUDAからOpenCLへの変換について述べる．ユーザは関数の開発に OpenCLLink を使ったことが想定されている．

OpenCLプログラムをCUDAにポートするときに必要なカーネルコードの一般的な変化を以下の表にまとめる．

入力画像を取り，色を反転する簡単なOpenCL関数である．

OpenCL関数を開始するには，次の式を実行しなければならない．

Out[5]=

CUDAにポートするために，次のものが変更される．

はに変更された．

はに変更された．

はに変更された．についても同じである．

多くの場合，Mathematica におけるOpenCLからCUDAへ４のポートの際に必要な変更は上記のものだけである．OpenCLまたはCUDAのプログラミングガイドに他の一対一の関数変換がいくつか記載されている．

次はポートされたCUDAコードである．

次はOpenCL用に行われた，Mathematica 関数そのままであるが，OpenCLFunctionLoadのロードをCUDAFunctionLoadで置き換える．

In[7]:=

Out[10]=

上で示したメカニズムを検索・置換する方法の代替として，カーネルコードを記号的に生成し，それを Mathematica で操作するというものがある．以下のセクションでその方法を示す．

記号的なコードの生成

Mathematica は，Cコードの階層的な表示を Mathematica 独自の言語として見ることを可能にするSymbolicCをサポートする．これにより，Mathematica はCコードの作成，操作，最適化に適したものとなる．

この機能に加え，CUDAカーネルコードをいくつかの別々の目的のために生成することができるため，よりよいポータブル性，より少ないプラットフォーム依存性，よりよいコード最適化等が可能となる．

SymbolicCを使うことの利点としては，Cコードが確実に文法エラーを含んでいないようにできることや，SymbolicCを使ってCUDAプログラムを生成し，CUDAFunctionLoadを使ってそれらを実行できること，そしてCUDAからOpenCLへの（またはその反対の）ポートが簡単になるということが挙げられる．

このセクションは CUDALink の記号コード生成機能を使って「RGB」から「HSV」，および「RGB」から「HSB」への色変換を書く．この変換は画像・映像処理において非常に一般的である．

このチュートリアルでは，まず最初に CUDALink アプリケーションをロードしなければならない．