特集記事

CUDA：「超並列コンピューティング」はじめました

επιστημη[著]

更新日: 2023/09/25
公開日: 2015/02/12

ポスト

ダウンロードサンプルファイル (184.0 KB)

行列のかけ算

　CUDAが得意な大量のデータ処理を行うサンプルとして行列の積を求めてみましょう。2つの行列：A（K行L列）とB（L行M列）の積C（K行M列）を求めます。

　float行列の積：C = A * Bにおいて、Cのy行x列：C[y][x]は Aのy行とBのx列との内積、すなわち

list04

float sum = 0.0f;
for ( unsigned int t = 0; t < l; ++t ) {
  sum += A[y][t] * B[t][x];
}
C[y][x] = sum;

　これをy=0,1,……K-1, x=0,1,……M-1の組み合わせ（K*M通り）について求めるわけで、行われるかけ算の総数はK*L*M回、K,L,Mが1000のとき入力要素数は2メガ個、1ギガ回のかけ算を行うことになります。

　やってみましょう、プロジェクト・テンプレートは既存アプリケーションのCUDA化を想定し"Win32 コンソールアプリケーション"を使います。

list05: mtxMul.cpp

#include <iostream>
#include <iomanip>
#include <chrono>

using namespace std;

/*
 * 行列の積:
 * C[k][m] = A[k][l] * B[l][m] 
 */
template<typename T>
void host_mtxMul(T* C, const T* A, const T* B, 
	             unsigned int k, unsigned int l, unsigned int m) {

  for ( unsigned int y = 0U; y < k; ++y ) {
    for ( unsigned int x = 0U; x < m; ++x ) {
      T sum = T(0);
      for ( unsigned int t = 0U; t < l; ++t ) {
        sum += A[y*l + t] * B[t*m + x];
      }
      C[y*m + x] = sum;
    }
  }
}

int main() {

  typedef float item;

  const unsigned int k = 1000U;
  const unsigned int l = 1000U;
  const unsigned int m = 1000U;

  item* A  = new item[k*l];
  item* B  = new item[l*m];
  item* hC = new item[k*m];

  for ( unsigned int i = 0; i < k*l; ++i ) A[i] = 0.001f * i;
  for ( unsigned int i = 0; i < l*m; ++i ) B[i] = 0.001f * i;

  chrono::high_resolution_clock::time_point start, stop;

  cout << "-----  host  ----- ";
  start = chrono::high_resolution_clock::now();
  host_mtxMul((item*)hC, (const item*)A, (const item*)B, k, l, m); // run on host
  stop = chrono::high_resolution_clock::now();
  chrono::milliseconds host_time = chrono::duration_cast<chrono::milliseconds>(stop - start);
  cout << setw(5) << host_time.count() << " [ms]" << endl;

  delete[] A;
  delete[] B;
  delete[] hC;

}

　僕の愛機：i7-2600K（Win8.1、x64）で14秒ほどかかりました。

　このプロジェクトにCUDA版行列積：device_mtxMulを実装します。メニュー：プロジェクト→新しい項目の追加でdevice_mtxMul.cuを追加し、メニュー：プロジェクト→ビルドのカスタマイズでCUDA 6.5にチェックを入れます。

　さらにdevice_mtxMul.cuのプロパティ：項目の種類を"ビルドに含めない"から"CUDA C/C++"に変更。

　もう一息、プロジェクトプロパティ：CUDA C++／Common／CUDA Toolkit Custon Dirを$(CUDA_PATH)とし、CUDAランタイム：cudart.libをリンクします。

　これで準備完了。device_mtxMulの実装に取り掛かりましょう。戦術は配列のたし算と同様、要素1つ分について計算するkernel_mtxMulをデバイス側に配置し、そいつをK*M個着火、つまりkernel_mtxMul<<<1,dim3(M,K)>>>(……)すればいい。dim3(x,y,z)はスレッド/ブロックを幅x：高さy：奥行zに配置します（省略時は1と見なされ、dim3(M,K)はdim3(M,K,1)と同じ）。

　……いや、残念ながらオトナの事情でそう単純ではありません。ブロックあたりのスレッド数には上限がありまして、大抵のGPUは1024個までなんです。CUDA SamplesにあるdeviceQueryを実行するとこの値を確認できます。

deviceQuery.exe

　なので1000*1000の配列ならば、例えばkernel_mtxMul<<<dim3(100,100),dim3(10,10)>>>(……);とします。10*10スレッドからなるブロックを100*100ブロック実行する、と。

　聞くところによりますとブロックあたり256スレッド程度がパフォーマンスが良いそうな。そんなこんなでデバイス側：kernel_mtxMulおよびデータの転送とkernel_mtxMulの呼び出しでホスト-デバイス間を仲介するdevice_mtxMulはこんなコードになりました。

list06 device_mtxMul.cu

#include <cuda_runtime.h>

#include <iostream>
#include <chrono>
#include <stdexcept>

using namespace std;

#pragma region both for device and host
/*
 * dim3: 要素ごとの 加減乗除 および 端数切り上げの除算
 */
__host__ __device__ inline dim3 operator+(const dim3& a, const dim3& b) 
  { return dim3(a.x + b.x, a.y + b.y, a.z + b.z); }

__host__ __device__ inline dim3 operator-(const dim3& a, const dim3& b) 
  { return dim3(a.x - b.x, a.y - b.y, a.z - b.z); }

__host__ __device__ inline dim3 operator*(const dim3& a, const dim3& b) 
  { return dim3(a.x * b.x, a.y * b.y, a.z * b.z); }

__host__ __device__ inline dim3 operator/(const dim3& a, const dim3& b) 
  { return dim3(a.x / b.x, a.y / b.y, a.z / b.z); }

__host__ __device__ inline dim3 roundup_div(const dim3& a, const dim3& b)
  { return dim3((a.x + b.x - 1)/b.x, (a.y + b.y - 1)/b.y, (a.z + b.z - 1)/b.z); }

/*
 * 配列を行列とみなすための wrapper
 */
template<typename T>
class ary2mtx {
  char*  base_;  // 先頭
  size_t pitch_; // 一行のバイト幅
public:
  __device__ ary2mtx(T* base, size_t pitch) : base_((char*)base), pitch_(pitch) {}
  __device__ T& operator()(unsigned int row, unsigned int column)
               { return *(T*)(base_ + row*pitch_) + column; }
  __device__ T* operator[](unsigned int row)
               { return (T*)(base_ + row*pitch_); }
};
#pragma endregion

#pragma region device
/*
 * C[y][x] = Σ(t:0..l-1) { A[y][t] * B[t][x] } 
 *  を device 側で計算する
 */
template<typename T>
__global__ void kernel_mtxMul(T* C, const T* A, const T* B, unsigned int k, unsigned int l, unsigned int m) {
  dim3 pos = blockDim * blockIdx + threadIdx;
  if ( pos.y < k && pos.x < m ) {
    ary2mtx<const T> a(A, sizeof(T)*l);
    ary2mtx<const T> b(B, sizeof(T)*m);
    ary2mtx<T>       c(C, sizeof(T)*m);
    T sum = T(0);
    for ( unsigned int t = 0; t < l; ++t ) {
      sum += a[pos.y][t] * b[t][pos.x];
    }
    c[pos.y][pos.x] = sum;
  }
}
#pragma endregion

#pragma region host
// dim3の体積
inline unsigned int volume(const dim3& d)
  { return d.x * d.y * d.z; }

// dim3の最も大きい要素を半減
inline void half(dim3& d) { 
  if ( d.x < d.y ) {
    if ( d.y < d.z ) d.z /= 2; else d.y /= 2;
  } else {
    if ( d.x < d.z ) d.x /= 2; else d.x /= 2;
  }
}

// cuda-API 失敗時に runtime_error を throw する
inline void assert_by_exception(cudaError_t result) {
  if ( result !=  cudaSuccess ) throw runtime_error(cudaGetErrorString(result));
}

/*
 * 行列の積:
 * C[k][m] = A[k][l] * B[l][m] 
 *  を device 側で行う
 */
template<typename T>
chrono::milliseconds device_mtxMul(T* C, const T* A, const T* B, 
                                   unsigned int k, unsigned int l, unsigned int m) {
  
  const T* dA = nullptr;
  const T* dB = nullptr;
  T*       dC = nullptr;

  chrono::high_resolution_clock::time_point start, stop;
  try {
    // device 側に領域確保
    assert_by_exception( cudaMalloc((void**)&dA, k*l*sizeof(T)) );
    assert_by_exception( cudaMalloc((void**)&dB, l*m*sizeof(T)) );
    assert_by_exception( cudaMalloc((void**)&dC, k*m*sizeof(T)) );

    // 入力データを host から device へ
    assert_by_exception( cudaMemcpy((void*)dA, A, k*l*sizeof(T), cudaMemcpyHostToDevice) );
    assert_by_exception( cudaMemcpy((void*)dB, B, l*m*sizeof(T), cudaMemcpyHostToDevice) );

    // thread と block の大きさを決定する
    dim3 desired(m,k);
    dim3 thr = desired;
    while ( volume(thr) > 256 ) half(thr);
    dim3 blk = roundup_div(desired,thr);

    // 実行し、完了を待つ
    start = chrono::high_resolution_clock::now();
    kernel_mtxMul<<<blk,thr>>>(dC, dA, dB, k, l, m);
    assert_by_exception( cudaGetLastError() );

    assert_by_exception( cudaDeviceSynchronize() );
    stop = chrono::high_resolution_clock::now();

    // 結果を device から host へ
    assert_by_exception( cudaMemcpy(C, dC, k*m*sizeof(T), cudaMemcpyDeviceToHost) );
  } catch ( const runtime_error& err ) {
    cerr << err.what() << endl;
  }

  // device 側領域の解放
  cudaFree((void*)dA);
  cudaFree((void*)dB);
  cudaFree(dC);
  return chrono::duration_cast<chrono::milliseconds>(stop - start);
}

chrono::milliseconds device_mtxMul_float(float* C, const float* A, const float* B, 
                                   unsigned int k, unsigned int l, unsigned int m) {
  return device_mtxMul(C,A,B,k,l,m);
}

chrono::milliseconds device_mtxMul_double(double* C, const double* A, const double* B, 
                                   unsigned int k, unsigned int l, unsigned int m) {
  return device_mtxMul(C,A,B,k,l,m);
}
#pragma endregion

　コード中に現れる__host__と__device__はそれぞれホスト側／デバイス側に展開されるコードであることを意味します（__global__,__host__,__device__のいずれも指定しないと__host__として扱われます）。デバイス側コードもC++で書けるのでC++屋の僕には快適です。デバイス側コードのコンパイラ、すなわちnvccはv6.5ではC++03に対応していますが、開発者登録すると手に入るv7.0RCだとC++11に対応しています。

　main()を少し修正し、ホストとデバイスとでスピードを比較してみました。

　ほほー、デバイスはホストの80倍、データ転送時間を除けば300倍のスピードです。floatのかけ算を1ギガ回やってますから、計算部分だけだと20GFLOPSを超えてます。CUDAすげー。

この記事は参考になりましたか？

連載通知を行うには会員登録(無料)が必要です。
既に会員の方はを行ってください。

印刷用を表示

ポスト

特集記事連載記事一覧: 【随時更新・まとめ読み】開発者のための「AWS re:Invent 2024」関連記事

開発者のための「WWDC 2024」関連発表まとめ

LeSSやSAFeなど、大規模アジャイルにおけるソフトウェア品質の課題とは？

もっと読む

この記事の著者: επιστημη（エピステーメー）

C++に首まで浸かったプログラマ。Microsoft MVP, Visual C++ (2004.01～2018.06) "だった"りわんくま同盟でたまにセッションスピーカやったり中国茶淹れてにわか茶...

※プロフィールは、執筆時点、または直近の記事の寄稿時点での内容です

この著者の最近の執筆記事

この記事は参考になりましたか？

この記事をシェア

ポスト

CUDA：「超並列コンピューティング」はじめました

目次

行列のかけ算

イベント