• トップ
  • G-DEPについて
  • ご購入ガイド
  • サポート
  • お問い合わせ

G-DEPトップ  >  第13回 NVIDIAが提供する数学ライブラリ

第13回 NVIDIAが提供する数学ライブラリ

 

第13回 「NVIDIAが提供する数学ライブラリ」

<< 第12回   |   目次に戻る   

 数学ライブラリとは、プログラム開発において頻出するルーチンをその道の専門家が作成したものです。非常に高速な上、ルーチンを呼出すだけで簡単に利用できるというメリットがあります。また、ライブラリ作成者がハードの進化に合わせてその都度最適化してくれるという点も大変大きなメリットです。GPUに対応したライブラリは既に多く開発されており、いくつかはNVIDIAのサイトで紹介されています。これらの内、CUFFTやCUBLAS、CUSPARSE等はNVIDIA自身が作成・提供しています。一方で、IMSLライブラリ等はサードパーティから提供されています。

今回は、この前者のNVIDIA自身が提供するGPU対応数学ライブラリを簡単にご紹介いたします。実は、CUDAをインストールすると、これらも同時に無償でインストールされます。これらはハードベンダーであるNVIDIAが自社のハードであるGPUに最適化したライブラリだけあって、抜群の最適化がなされています。今後のGPUの進化にも適時対応してくれる事間違いなしのライブラリです。今回は、密ベクトル・密行列の線形演算を扱うライブラリとして有名なBLASのGPU版 "CUBLAS" をご紹介します。なお、CUBLASのガイドはこちらでご覧になれますし、CUDAフォルダ内のdocにもpdfとして置かれています。一旦CUBLASに慣れてしまえば、他のライブラリも同じ要領で使いこなせます。

 

13.1 行列積のコード -Dgemmの利用例-


 一般にBLASのルーチンは、演算タイプに応じて以下の3つのレベルに分類されます。ここでα, βはスカラー、x, yはベクトル、A, B, Cは行列です。

 レベル 1: y ← αx+y

 レベル2: y ← αAx+βy

 レベル 3: C ← αAB+βC

また、簡単な指示により、行列Aは転置させたり、エルミート共役にしたりできます。

以下では、この内のレベル3に分類されるDgemmルーチンを使ってみましょう。CUBLASではcublasDgemmというルーチン名になっています。大文字のDはdouble、つまり実倍精度演算という意味です。実単精度にはcublasSgemm、複素単精度にはcublasCgemm、複素倍精度にはcublasZgemmというルーチンが用意されています。BLASを使ったことがある方は分かるかと思いますが、BLASの関数にcublas~と付けるだけになります。

なお、CUDA 3.2まではCUBLAS APIを使うにはヘッダファイル "cublas.h" をインクルードする必要がありましたが、CUDA 4.0以降は一新されたAPIが "cublas_v2.h" で提供されています。以下ではこちらを利用します。(その特徴はこちらにまとめられています。)

cublas_v2.hでのcublasDgemmは、以下で与えられます:

 

cublasDgemm( cublasHandle_t handle,

         cublasOperation_t transa, cublasOperation_t transab,

         int m, int n, int k,

         const double *alpha,

         const double *A, int lda,

         const double *B, int ldb,

         const double *beta,

         double *C, int ldc)

 

一見、引数が沢山あって複雑に見えますが、行列積で必要最小限の情報だけを渡しています。慣れてしまえば簡単です。このルーチンは次のような演算を行えます:

C ← α op(A) op(B) + βC

op(A)はtransaで、op(B)はtransbで決まる変換です。

        A         if  transa == CUBLAS_OP_N

op(A) =   Aの転置行列    if  transa == CUBLAS_OP_T

   Aのエルミート共役      if  transa == CUBLAS_OP_C

整数m ,n, kはop(A)が m x k行列、op(B)が k x n行列、Cがm x n行列であるという指示です。CUBLASでは通常のBLASと同じく、行列のデータは1次元配列に格納されます。ここで注意すべきは、CUBLASではデータレイアウトがcolumn-majorであるという点です。例えば、1次元配列 A=(1, 2, 3, .., 2048x2048) を 2048x2048行列としてルーチンに渡すと、デフォルト(CUBLAS_OP_N)では

なる行列として解釈されます。それでは早速、ライブラリをつかったコードを書いてみましょう。

C ← α op(A) op(B) + βC

について、行列や定数の定義、デバイスメモリの確保など、途中まではライブラリを使わない場合(第9回 「プログラミング演習 ~基礎編②~」をご参照ください)と同じになります。具体的には

 

 //cublas_v2.hをインクルード

#include "cublas_v2.h"

 

//Nを定義

#define N 2048

 

int main() {

 int i;

 double *A, *B, *C;

 

 //デバイス(GPU)側での名称を定義

 double *Ad, *Bd, *Cd, alpha=3.0, beta=2.0;

 cublasHandle_t handle;

 cublasCreate(&handle);

 

 //メインメモリの確保

 A = (double*)malloc(N*N*sizeof(double));

 B = (double*)malloc(N*N*sizeof(double));

 C = (double*)malloc(N*N*sizeof(double));

 

 //A,B,Cを定義(#1)

  for (i=0;i<N*N;i++) A[i] = B[i] = (double) i;

  for (i=0;i<N*N;i++) C[i] = 0.0;

 

 //デバイス側でのメモリー確保

 cudaMalloc((void**)&Ad, N*N*sizeof(double));

 cudaMalloc((void**)&Bd, N*N*sizeof(double));

 cudaMalloc((void**)&Cd, N*N*sizeof(double));

 

 //タイマーを生成して、計測開始

 unsigned int timer =0;

 CUT_SAFE_CALL( cutCreateTimer( &timer));

 CUT_SAFE_CALL( cutStartTimer( timer));

 

 //デバイスへのメモリ転送

 cublasSetMatrix(N, N, sizeof(double), A, N, Ad, N);

 cublasSetMatrix(N, N, sizeof(double),B, N, Bd, N);

 

 //C ← α op(A) op(B) + βCの計算

 cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, N, N, &alpha,

         Ad, N, Bd, N, &beta, Cd, N);

 

 //メインへのメモリ転送

 cublasGetMatrix(N,N,sizeof(double), Cd, N, C, N);

 

 //タイマーを停止し、時間を表示

 CUT_SAFE_CALL( cutStopTimer( timer));

 Printf("Procesing time: %f (msec)\n",cutGetTimerValue( timer));

 CUT_SAFE_CALL( cutDeleteTimer( timer));

 

 //デバイス側のメモリを解放

 cudaFree(Ad);

 cudaFree(Bd);

 cudaFree(Cd);

 cublasDestroy(handle);

 

 //メイン側のメモリを解放

 free(A);

 free(B);

 free(C);

 

 return 0;

}

 

#1 前述のようにデータレイアウトがcolumn-majorなので

と入力されています。

以上が、ライブラリを使用したプログラミングになります。それでは次にちゃんとGPUがつかえているのか比較してみましょう。

 

13.2 実際に試してみる


 では、コンパイルして実行してみましょう。今回は

1. CPU 1コアのみ

2. Intel MKLを使用

3. naiveなCUDA行列積コード(ライブラリ不使用)

4. CUBLASを使用

の4つの場合を比較してCUDAライブラリの有用性について考えてみます。

行うのは、先ほどと同じ

C ← α op(A) op(B) + βC

として、α=2, β=3を代入してあります。また、行列については 2048×2048 の正方行列で、各要素についてはランダムの実数を当てはめるようにしています。

☆計算環境☆
CPU:Intel Core i7-950 3.07GHz
GPU:NVIDIA Geforce GTX480
OS:Windows 7 64bit

それぞれのコードは以下の通りです。

 

1. CPU 1コアのみ

// CPU Matrix
 
#include<stdio.h>
#include<malloc.h>
#include<stdlib.h>
#include<time.h>
 
#define N 2048
 
int main( int argc, char**argv) {
 
 int col, row, scan;
 
 int*matA;
 int*matB;
 int*matC;
 
 time_t timeStart, timeStop;
 
 //行列のメインメモリを確保
 matA = (int*)malloc(sizeof(int) * N * N);
 matB = (int*)malloc(sizeof(int) * N * N);
 matC = (int*)malloc(sizeof(int) * N * N);
 
 //行列の各要素にランダムの実数値を挿入
 for(col = 0; col<N; col++) {
  for(row = 0; row<N; row++) {
   matA[col*N+row]=rand()%(N * N);
   matB[col*N+row]=rand()%(N * N);
   matC[col*N+row]=rand()%(N * N);
  }
 }
 
 //タイマーをスタート
 time(&timeStart);
 
 //C=3*A*B+2*C を計算
 for (col = 0; col<N; col++) {
    for (row =0; row<N; row++) {
      for (scan = 0; scan<N; scan++) {
         matC[col*N+row] += 3*matA[col*N+scan]*matB[scan*N+row]+2*matC[col*N+row];
   }
  }
 }
 
 //タイマーをストップ
 time(&timeStop);
 
 //結果を表示
 printf("It takes %dsec\n",timeStop - timeStart);
 
 //メモリを解放
 free(matA);
 free(matB);
 free(matC);
 
 return 0;
}

さて、これを実行してみますと、


 

CPU 1コアの場合は77秒です。

次に、非常に高速で演算できるとIntel MKL BLASを用いた演算について見てみましょう。

 

2. Intel MKLを使用

// Using MKL

#include <stdio.h>
#include <malloc.h>
#include <time.h>
#include "mkl.h"

#define N 2048 // 正方行列のサイズを指定(N×N)

int main( int argc, char** argv)
{
    int col, row, scan;
    double *matA, *matB, *matC;
    double alpha = 3.0;
    double beta  = 2.0;

    // 時間計測
    clock_t Start, Stop;

    // 行列A~Cの動的メモリ確保
    matA = (double*)malloc(sizeof(double) * N * N);
    matB = (double*)malloc(sizeof(double) * N * N);
    matC = (double*)malloc(sizeof(double) * N * N);

    // 行列A、Bの各要素に乱数を入れる
    // 行列Cは初期値0

    for (col = 0; col < N; col++) {
        for (row = 0; row < N; row++) {
            matA[col * N + row] = rand() % (N * N);
            matB[col * N + row] = rand() % (N * N);
            matC[col * N + row] = 0;
        }
    }


    // 時間計測開始
    Start=clock();

    // MKL
   cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
         N, N, N, alpha, matA, N, matB, N, beta, matC, N);

    // 時間計測終了
    Stop=clock();

    // 計算時間結果の表示
    printf("Using MKL\n It takes %dmsec\n", Stop - Start);

    // 各行列のメモリ開放
    free(matA);
    free(matB);
    free(matC);

    return 0;
}

これを実行してみますと、

 

0.7秒と、ライブラリ(MKL)を使用することで大幅に計算を短縮することができます。

 

では実際にGPU計算はどれだけ大きな効力を発揮するのでしょう。GPUを使用した2つの場合を比較してみます。

 

3. naiveなCUDA行列積コード(ライブラリ不使用)

#include <stdio.h>

#include <malloc.h>

#include <stdlib.h>

#include <cutil_inline.h>

#include "cublas_v2.h"

 

#define N 2048 // 正方行列のサイズを指定(N×N)

#define BLOCK 16// ブロックのサイズを指定

 

__global__ void

matrixMul(double* inMatA, double* inMatB, double* inMatC);

 

int main(int argc, char** argv) {

 

 // 行列のサイズをバイト単位で算出

 int matrixSize = sizeof(double) * N * N;

 

 double *hMatA, *hMatB, *hMatC; // ホスト側の行列変数設定

 double *dMatA, *dMatB, *dMatC;// デバイス側の行列変数設定

 double alpha=3.0, beta=2.0;

 

 printf("MatrixSize:%d*%d\n\n", N, N);

 

 // 行列変数のメモリ確保

 hMatA = (double*)malloc(matrixSize);

 hMatB = (double*)malloc(matrixSize);

 hMatC = (double*)malloc(matrixSize);

 

 // 初期値設定

 int col, row;

 

 for (col = 0; col < N; col++) {

    for (row = 0; row < N; row++) {

   hMatA[col * N + row] = rand() % (N * N);

   hMatB[col * N + row] = rand() % (N * N);

  }

 }

 

 // ブロックサイズとグリッドサイズの設定

 dim3 block(BLOCK, BLOCK);

 dim3 grid( N / BLOCK, N / BLOCK);

 cublasHandle_t handle;

 

 // デバイスメモリ領域の確保

 cudaMalloc((void**)&dMatA, matrixSize);

 cudaMalloc((void**)&dMatB, matrixSize);

 cudaMalloc((void**)&dMatC, matrixSize);

 

 // タイマーを作成して計測開始

 unsigned int timer = 0;

 cutCreateTimer( &timer);

 cutStartTimer( timer);

 

 // ホストからデバイスへの変数の受け渡し

 cudaMemcpy(dMatA, hMatA, matrixSize, cudaMemcpyHostToDevice);

 cudaMemcpy(dMatB, hMatB, matrixSize, cudaMemcpyHostToDevice);

 

 // カーネルの起動

 matrixMul<<<grid, block>>>(dMatA, dMatB, dMatC);

 cudaThreadSynchronize();

 

 // 結果の領域確保とデバイス側からのメモリ転送

 cudaMemcpy(hMatC, dMatC, matrixSize, cudaMemcpyDeviceToHost);

 

 // 時間計測終了

 CUT_SAFE_CALL( cutStopTimer( timer));

 

 // 計算時間結果の表示

 printf("Using kernel(non CUBLAS)\n it takes %f (msec)\n\n", cutGetTimerValue( timer));

 cutDeleteTimer( timer);

 

 cudaFree(dMatA);

 cudaFree(dMatB);

 cudaFree(dMatC);

 

__global__ void 

matrixMul(double* inMatA, double* inMatB, double* inMatC) {

 

 int col = blockIdx.x * blockDim.x + threadIdx.x;

 int row = blockIdx.y * blockDim.y + threadIdx.y;

 int scan;

 int target = 0;

 

 // 行列の演算を行う

 for (scan = 0; scan < N; scan++) {

  target += inMatA[col * N + scan] * inMatB[scan * N + row];

  __syncthreads();

 }

 

 inMatC[col * N + row] = target;

}

 

 

続けて、CUBLASを使用した場合です。

. CUBLASを使用

#include <stdio.h>

#include <malloc.h>

#include <stdlib.h>

#include <cutil_inline.h>

#include "cublas_v2.h"

 

#define N 2048 // 正方行列のサイズを指定(N×N)

#define BLOCK 16// ブロックのサイズを指定

 

__global__ void

matrixMul(double* inMatA, double* inMatB, double* inMatC);

 

int main(int argc, char** argv) {

 

 // 行列のサイズをバイト単位で算出

 int matrixSize = sizeof(double) * N * N;

 

 double *hMatA, *hMatB, *hMatC; // ホスト側の行列変数設定

 double *dMatA, *dMatB, *dMatC;// デバイス側の行列変数設定

 double alpha=3.0, beta=2.0;

 

 printf("MatrixSize:%d*%d\n\n", N, N);

 

 // 行列変数のメモリ確保

 hMatA = (double*)malloc(matrixSize);

 hMatB = (double*)malloc(matrixSize);

 hMatC = (double*)malloc(matrixSize);

 

 // 初期値設定

 int col, row;

 

 for (col = 0; col < N; col++) {

     for (row = 0; row < N; row++) {

   hMatA[col * N + row] = rand() % (N * N);

   hMatB[col * N + row] = rand() % (N * N);

     }

 }

 

 // ブロックサイズとグリッドサイズの設定

 dim3 block(BLOCK, BLOCK);

 dim3 grid( N / BLOCK, N / BLOCK);

 cublasHandle_t handle;

 

 cublasCreate(&handle);

 

 // デバイスメモリ領域の確保

 cudaMalloc((void**)&dMatA, matrixSize);

 cudaMalloc((void**)&dMatB, matrixSize);

 cudaMalloc((void**)&dMatC, matrixSize);

 

 // タイマーを作成して計測開始

 cutCreateTimer( &timer);

 cutStartTimer( timer);

 

 //デバイスへのメモリ転送

  cublasSetMatrix(N,N,sizeof(double),hMatA,N,dMatA,N);

  cublasSetMatrix(N,N,sizeof(double),hMatB,N,dMatB,N);

 

 // ライブラリ計算

 cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, N, N, &alpha,

         dMatA, N, dMatB, N, &beta, dMatC, N);

 

 // 結果の領域確保とデバイス側からのメモリ転送

 cublasGetMatrix(N,N,sizeof(double),dMatC,N,hMatC,N);

 

 // 時間計測終了

 CUT_SAFE_CALL( cutStopTimer( timer));

 

 cublasDestroy(handle);

 

 // 計算時間結果の表示

 printf("Using CUBLAS\n it takes %f (msec)\n\n", cutGetTimerValue( timer));

 cutDeleteTimer( timer);

 

 // ホスト・デバイスメモリの解放

 free(hMatA);

 free(hMatB);

 free(hMatC);

 cudaFree(dMatA);

 cudaFree(dMatB);

 cudaFree(dMatC);

 

 // 終了処理

 cudaThreadExit();

 cutilExit(argc, argv);

}

 

__global__ void

matrixMul(double* inMatA, double* inMatB, double* inMatC){

  

 int col = blockIdx.x * blockDim.x + threadIdx.x;

 int row = blockIdx.y * blockDim.y + threadIdx.y;

 int scan;

 int target = 0;

 

 // 行列の演算を行う

 for (scan = 0; scan < N; scan++) {

  target += inMatA[col * N + scan] * inMatB[scan * N + row];

  __syncthreads();

 }

 

 inMatC[col * N + row] = target;

}

 

さてお待ちかねの結果を見てみましょう。(以下の結果は上記コードをひとまとめにしたプログラムで表示しています。)

 

CUBLASを使用していない場合は1.45秒とIntel MKLの2倍ぐらいかかってしまいますが、CUBLASを使用した場合はなんと0.16秒になっており、断然CUBLASを使用したほうが早いことがわかります。

また、ライブラリを使うことでソースコードが見やすくなるといった点も良いですね。

 

以上が今回のCUBLASを使用した場合の例になります。

実際に計算してみるとよりその速さを実感できますので、ぜひともお試しくださいね!

 

<< 第12回   |   目次に戻る   


参考資料

 CUBLAS LIBRALY, User Guide (http://docs.nvidia.com/cuda/cublas/index.html)

 

(執筆: Associate Research Engineer 東京大学大学院工学系研究科 岡安優、久保田英司 共著)