cuFFT: フーリエ変換で雑音を消す
フーリエ変換(Fourier Transform)は信号処理/解析のド定番。
周期を持ったあらゆる波は異なる周波数のサイン波の重ね合わせで作り出すことができ、与えられた波形から、その成分(サイン波)を割り出すのがフーリエ変換、CUDAにはフーリエ変換ライブラリ: cuFFT が入ってます。軽く使ってみましょうね。
まず入力波を作ります。350,400,450Hzのサイン波を重ね合わせ、さらに一様乱数で生成したノイズを乗せましょう。
const size_t N = 44100U; // データ数(44.1kHzサンプリングでの1秒分) vector<float> signal(N); vector<float> h_in(N); // 振幅 ±2 のホワイト・ノイズ mt19937 mt; uniform_real_distribution<float> rnd(-2.0f, 2.0f); // 350,400,450Hzのサイン波にノイズを乗せる float omega = 2.0f * 3.1416f / N; for ( unsigned int i = 0; i < N; ++i ) { signal[i] = sinf(omega * 350.0f * (float)i) * 1.0f + sinf(omega * 400.0f * (float)i) * 0.8f + sinf(omega * 450.0f * (float)i) * 0.6f ; h_in[i] = signal[i] + rnd(mt); }
こんなのができました。'赤'はサイン波の合成、それにノイズを乗せたのが'青'です。
この(ノイズまみれの)信号にフーリエ変換を施します。
// device-memoryの確保(入/出力兼用) float* d_real = nullptr; cudaMalloc(&d_real, N*sizeof(float)); float2* d_cplx = reinterpret_cast<float2*>(d_real); // フーリエ変換 cudaMemcpy(d_real, h_in.data(), N*sizeof(float), cudaMemcpyHostToDevice); cufftHandle plan_f; cufftPlan1d(&plan_f, N, CUFFT_R2C, 1); // Real to Complex (forward) cufftExecR2C(plan_f, d_real, d_cplx); vector<float2> h_mid(N/2); // スペクトル(フーリエ変換の結果) cudaMemcpy(h_mid.data(), d_cplx, N*sizeof(float), cudaMemcpyDeviceToHost);
変換結果がコレ。
350,400,450に大きなピークが見られますね。ノイズは様々な周波数の波がちょっとずつ重なったものなのでグラフの底に貼りつく'モジョモジョ'した部分に現れます。
で、このデータから300Hz以下と500Hz以上の部分をばっさり削ってしまいます。帯域フィルタ(band-pass filter)ってやつです。
// band-pass filter // 300Hz以下/500Hz以上の信号をカットする cudaMemset(d_cplx , 0, 300U * sizeof(float2)); cudaMemset(d_cplx+500U, 0, (N/2-500U) * sizeof(float2));
よーするに邪魔なノイズ成分の多くを削り取ったことになります。
しかるのち逆フーリエ変換をかけて、周波数軸から時間軸に戻します。
// 逆フーリエ変換 cufftHandle plan_i; cufftPlan1d(&plan_i, N, CUFFT_C2R, 1); // Complex to Real (inverse) cufftExecC2R(plan_i, d_cplx, d_real); // 結果の出力 vector<float> h_out(N); cudaMemcpy(h_out.data(), d_real, N*sizeof(float), cudaMemcpyDeviceToHost);
結果がコレ。ノイズが消えました♪
コチラ↓が全コード:
/* * Noise Reduction with cuFFT */ #include <cuda_runtime.h> #include <cufft.h> #include <iostream> #include <random> #include <vector> #include <cmath> using namespace std; int main() { const size_t N = 44100U; // データ数(44.1kHzサンプリングでの1秒分) vector<float> signal(N); vector<float> h_in(N); // 振幅 ±2 のホワイト・ノイズ mt19937 mt; uniform_real_distribution<float> rnd(-2.0f, 2.0f); // 350,400,450Hzのサイン波にノイズを乗せる float omega = 2.0f * 3.1416f / N; for ( unsigned int i = 0; i < N; ++i ) { signal[i] = sinf(omega * 350.0f * (float)i) * 1.0f + sinf(omega * 400.0f * (float)i) * 0.8f + sinf(omega * 450.0f * (float)i) * 0.6f ; h_in[i] = signal[i] + rnd(mt); } // device-memoryの確保(入/出力兼用) float* d_real = nullptr; cudaMalloc(&d_real, N*sizeof(float)); float2* d_cplx = reinterpret_cast<float2*>(d_real); // フーリエ変換 cudaMemcpy(d_real, h_in.data(), N*sizeof(float), cudaMemcpyHostToDevice); cufftHandle plan_f; cufftPlan1d(&plan_f, N, CUFFT_R2C, 1); // Real to Complex (forward) cufftExecR2C(plan_f, d_real, d_cplx); vector<float2> h_mid(N/2); // スペクトル(フーリエ変換の結果) cudaMemcpy(h_mid.data(), d_cplx, N*sizeof(float), cudaMemcpyDeviceToHost); // band-pass filter // 300Hz以下/500Hz以上の信号をカットする cudaMemset(d_cplx , 0, 300U * sizeof(float2)); cudaMemset(d_cplx+500U, 0, (N/2-500U) * sizeof(float2)); // 逆フーリエ変換 cufftHandle plan_i; cufftPlan1d(&plan_i, N, CUFFT_C2R, 1); // Complex to Real (inverse) cufftExecC2R(plan_i, d_cplx, d_real); // 結果の出力 vector<float> h_out(N); cudaMemcpy(h_out.data(), d_real, N*sizeof(float), cudaMemcpyDeviceToHost); cout << "signal, noised, processed, spectrum" << endl; for ( unsigned int i = 0; i < 500; ++i ) { cout << signal[i] << ',' << h_in[i] << ',' << h_out[i]/N << ',' << cuCabsf(h_mid[i]) << endl; } cudaFree(d_real); cufftDestroy(plan_f); cufftDestroy(plan_i); }
(original: 2015-05-21 #36 #37)
cuSOLVER: 鶴亀算を解く
「つる と かめ があわせて3匹います。
足の数はあわせて10本でした。
問1: つる と かめ はそれぞれ何匹ですか?
問2: 足の数があわせて8本なら、それぞれ何匹ですか?」
小学校の算数でやりましたよね、鶴亀算。
コ難しくいえば 連立一次方程式の解を求める問題です。
行列で表現するとこんな。
与えられた A と B から X を求めよ、ってわけ。
CUDA 7.0からこんな問題をCUDAで解いてくれるライブラリ: cuSOLVER がついてます。早速使ってみましたよ。
#include <cuda_runtime.h> #include <cusolverDn.h> // dense LAPACK #include <cassert> #include <iostream> using namespace std; template<typename T> inline size_t bytesof(unsigned int s) { return s * sizeof(T); } template<typename T> T* allocate(unsigned int size) { T* result = nullptr; cudaError_t status = cudaMalloc(&result, bytesof<T>(size)); assert( status == cudaSuccess ); return result; } int main() { cusolverStatus_t status; // dense LAPACK cusolverDnHandle_t handle; status = cusolverDnCreate(&handle); assert( status == CUSOLVER_STATUS_SUCCESS ); int n = 2; // 2x2 正方行列 float A[] = { 1.0f, 1.0f, 2.0f, 4.0f }; float* dA = allocate<float>(n*n); cudaMemcpy(dA, A, bytesof<float>(n*n), cudaMemcpyHostToDevice); int lda = 2; // 必要なバッファ量を求め、確保する int worksize; status = cusolverDnSgetrf_bufferSize( handle, n, // 行 n, // 列 dA, // A lda, // Aのヨコハバ &worksize); assert( status == CUSOLVER_STATUS_SUCCESS ); #ifdef _DEBUG cout << "worksize = " << worksize << endl; #endif float* workspace = allocate<float>(worksize); // 計算結果に関する情報 int* devInfo = allocate<int>(1); // ピボット int* pivot = allocate<int>(n); // LU分解 : dAに結果が求まる(それとpivot) status = cusolverDnSgetrf( handle, n, // 行 n, // 列 dA, // A lda, // Aのヨコハバ workspace, pivot, devInfo); #ifdef _DEBUG int info; cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost); cout << "info = " << info << endl; #endif assert( status == CUSOLVER_STATUS_SUCCESS ); // 鶴と亀の総数, 足の総数 float B[] = { 3.0f, 10.0f , 3.0f, 8.0f }; int nrhs = 2; // 問題数 float* dB = allocate<float>(n*nrhs); cudaMemcpy(dB, B, bytesof<float>(n*nrhs), cudaMemcpyHostToDevice); int ldb = 2; // AX = B を解く (解XはBをoverrideする) status = cusolverDnSgetrs( handle, CUBLAS_OP_T, n, // 行(=列) nrhs, // 問題数 dA, // A lda, // Aのヨコハバ pivot, // LU分解で得られたピボット dB, // B ldb, // Bのヨコハバ devInfo); #ifdef _DEBUG cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost); cout << "info = " << info << endl; #endif assert( status == CUSOLVER_STATUS_SUCCESS ); // 結果を取得し、出力する float X[16]; cudaMemcpy(X, dB, bytesof<float>(n*nrhs), cudaMemcpyDeviceToHost); for ( int i = 0; i < nrhs; ++i ) { float* q = B + i*2; float* a = X + i*2; cout << "総数 = " << q[0] << ", 足の数 = " << q[1] << "\t 解 : " <<"鶴 = " << a[0] << " , 亀 = " << a[1] <<endl; } cudaFree(workspace); cudaFree(dA); cudaFree(dB); cudaFree(devInfo); cudaFree(pivot); cusolverDnDestroy(handle); }
実行結果:
(original: 2016-06-03 #40)
streamで速くする(3) ~からくり
CPUとGPUはPCI-busに隔てられてそれぞれが勝手に動くことができます。
CPUはGPUの仕事の完了を待って次の仕事を依頼する必要はないんですわ(そうでないとGPUが重たい仕事してる間CPUがぼーっと待ってにゃならんですから)。仕事の完了を待たずに次の仕事を叩き込むことができるってことは、仕事の待ち行列が用意されてるってことで、それがstreamです。
streamに溜まった仕事(メモリ・コピー と kernel実行)はふたつのエンジンが捌きます。kernel実行をCUDA-core群に割り当てるGigaThread Engineと、メモリ・コピーの依頼を受けてPCI-bus経由でCPU/GPU間のデータ転送を司るCopy Engine、この二つは独立して動けます。
んだから、kernel実行中には次のkernel実行は待たされ、メモリ・コピー中は次のメモリ・コピーは待たされるけど、kernel実行とメモリ・コピーとはそれぞれ異なるEngineで処理されるために同時にやれるってわけですわ。
GeForce GTX9xx, GTX10xx になるとCPU→GPU用 と GPU→CPU用のふたつのCopy Engineを持ってます。PCI-busは上り/下りを同時に転送できるのでより速くなります。
GTX980でのtimelineはこんなカンジ:
streamで速くする(2)
streamはその名の通り"一連の(処理の)流れ"です。 streamに乗せられたjob(仕事)は並んだ順に実行されます。 複数のstreamにjobを乗せておけば、GPUはその中から実行可能なものを一つ選んで処理します。 ホスト・メモリ間のデータ転送にGPUは関与しないので、データ転送の間でも他のstreamに乗ったjobを処理できます。
ここまでにいくつかサンプル・コードを紹介してきましたが、それらの中には stream を明示的に使ったものはありません。 けども実は暗黙の default stream が一本だけあって、その default stream に cudaMemcpy
やら kernel呼び出しやらのjobを乗っけていたんです。
3本のstream(ssin, scos, stan)を用意し、sin,cos,tanの表作成をそれぞれの stream に乗せて処理しましょう。
#include <cuda_runtime.h> #include <device_launch_parameters.h> #include <math_functions.h> // 時間稼ぎ __device__ float time_waste() { float result = 0.0f; const int TIMES = 10U; for ( int i = 0; i < TIMES; ++i ) { result += 1.0f / TIMES; } return result; } // 度→ラジアン __device__ __forceinline__ float deg2rad(float angle) { return angle / 180.0f * 3.01416f * time_waste(); } // sine __global__ void kernel_sin(float* out, const float* angle, unsigned int size) { unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { out[i] = sinf(deg2rad(angle[i])); } } // cosine __global__ void kernel_cos(float* out, const float* angle, unsigned int size) { unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { out[i] = cosf(deg2rad(angle[i])); } } // tangent __global__ void kernel_tan(float* out, const float* angle, unsigned int size) { unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { out[i] = tanf(deg2rad(angle[i])); } } // helper funcs template<typename T> inline T* allocate_host(size_t size) { T* p = nullptr; cudaMallocHost(&p, size*sizeof(T)); return p; } template<typename T> inline void free_host(T* host) { cudaFreeHost(host); } template<typename T> inline T* allocate_device(size_t size) { T* p = nullptr; cudaMalloc(&p, size*sizeof(T)); return p; } template<typename T> inline void free_device(T* device) { cudaFree(device); } inline cudaEvent_t create_event() { cudaEvent_t e = nullptr; cudaEventCreate(&e); return e; }; inline void destroy_event(cudaEvent_t event) { cudaEventDestroy(event); } inline cudaStream_t create_stream() { cudaStream_t s = nullptr; cudaStreamCreate(&s); return s; }; inline void destroy_stream(cudaStream_t stream) { cudaStreamDestroy(stream); } #include <iostream> using namespace std; const unsigned int N = 1000U * 1000U; const unsigned int BN = N*sizeof(float); void multiple_stream() { cout << "--- multiple stream ---" << endl; float* angle = allocate_host<float>(N); // in float* sin = allocate_host<float>(N); // out float* cos = allocate_host<float>(N); // out float* tan = allocate_host<float>(N); // out for ( unsigned int i = 0; i < N; ++i ) { angle[i] = 360.0f * i / N; } float* dangle = allocate_device<float>(N); float* dsin = allocate_device<float>(N); float* dcos = allocate_device<float>(N); float* dtan = allocate_device<float>(N); unsigned int size = N; unsigned int block = 128U; unsigned int grid = (size + block -1U) / block; cudaEvent_t start = create_event(); cudaEvent_t mid = create_event(); cudaEvent_t stop = create_event(); cudaStream_t ssin = create_stream(); cudaStream_t scos = create_stream(); cudaStream_t stan = create_stream(); cudaEventRecord(start); cudaMemcpy(dangle, angle, BN, cudaMemcpyHostToDevice); cudaEventRecord(mid); kernel_sin<<<grid,block,0,ssin>>>(dsin, dangle, size); cudaMemcpyAsync(sin, dsin, BN, cudaMemcpyDeviceToHost, ssin); kernel_cos<<<grid,block,0,scos>>>(dcos, dangle, size); cudaMemcpyAsync(cos, dcos, BN, cudaMemcpyDeviceToHost, scos); kernel_tan<<<grid,block,0,stan>>>(dtan, dangle, size); cudaMemcpyAsync(tan, dtan, BN, cudaMemcpyDeviceToHost, stan); cudaEventRecord(stop); cudaDeviceSynchronize(); float elapsed; cudaEventElapsedTime(&elapsed, start, mid); cout << elapsed << " + "; cudaEventElapsedTime(&elapsed, mid, stop); cout << elapsed << " = "; cudaEventElapsedTime(&elapsed, start, stop); cout << elapsed << " [ms]" << endl; for ( unsigned int i = 0; i < 5; ++i ) { cout << angle[i] << " :" << " sin=" << sin[i] << " cos=" << cos[i] << " tan=" << tan[i] << endl; } cout << "..." << endl; for ( auto event : { start, mid, stop }) destroy_event(event); for ( auto stream : { ssin, scos, stan }) destroy_stream(stream); for ( auto host : { angle, sin, cos, tan }) free_host(host); for ( auto device : { dangle, dsin, dcos, dtan }) free_device(device); } int main() { multiple_stream(); }
default stream 1本で処理した時のタイムラインがコレ。
対して3本(+default)のstreamでやったのがコチラ。
ね、データ転送とkernel実行が同時に行われてるでしょ。重なった分だけ速くなるんです。
実装上のキモはふたつ:
cudaMemcpy
じゃなくcudaMemcpyAsync
を使うべし。 後者はコピー完了を待たずに(job投入後ただちに)返って来るので、すぐさま次のjobを投入できます。ホスト側メモリを
malloc
/free
やnew
/delete
じゃなくcudaMallocHost
/cudaFreeHost
で確保/解放すべし。 そうしないとデータ転送中に kernelが動いてくれないのですよ。
... streamのおハナシはまだまだ続きます。
(original: 2015-04-16)
streamで速くする(1) ~イントロ
三角関数表を作りましょう。 ひと回り360度を百万等分し、360マイクロ度刻みの sin, cos, tan を求めます。
#include <cuda_runtime.h> #include <device_launch_parameters.h> #include <math_functions.h> // 時間稼ぎ __device__ float time_waste() { float result = 0.0f; const int TIMES = 10U; for ( int i = 0; i < TIMES; ++i ) { result += 1.0f / TIMES; } return result; } // 度→ラジアン __device__ __forceinline__ float deg2rad(float angle) { return angle / 180.0f * 3.01416f * time_waste(); } // sine __global__ void kernel_sin(float* out, const float* angle, unsigned int size) { unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { out[i] = sinf(deg2rad(angle[i])); } } // cosine __global__ void kernel_cos(float* out, const float* angle, unsigned int size) { unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { out[i] = cosf(deg2rad(angle[i])); } } // tangent __global__ void kernel_tan(float* out, const float* angle, unsigned int size) { unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { out[i] = tanf(deg2rad(angle[i])); } } // helper funcs template<typename T> inline T* allocate_host(size_t size) { T* p = nullptr; cudaMallocHost(&p, size*sizeof(T)); return p; } template<typename T> inline void free_host(T* host) { cudaFreeHost(host); } template<typename T> inline T* allocate_device(size_t size) { T* p = nullptr; cudaMalloc(&p, size*sizeof(T)); return p; } template<typename T> inline void free_device(T* device) { cudaFree(device); } inline cudaEvent_t create_event() { cudaEvent_t e = nullptr; cudaEventCreate(&e); return e; }; inline void destroy_event(cudaEvent_t event) { cudaEventDestroy(event); } #include <iostream> using namespace std; const unsigned int N = 1000U * 1000U; const unsigned int BN = N*sizeof(float); void single_stream() { cout << "--- single stream ---" << endl; // host memories float* angle = allocate_host<float>(N); // in float* sin = allocate_host<float>(N); // out float* cos = allocate_host<float>(N); // out float* tan = allocate_host<float>(N); // out for ( unsigned int i = 0; i < N; ++i ) { angle[i] = 360.0f * i / N; } // device memories float* dangle = allocate_device<float>(N); float* dsin = allocate_device<float>(N); float* dcos = allocate_device<float>(N); float* dtan = allocate_device<float>(N); unsigned int size = N; unsigned int block = 128U; unsigned int grid = (size + block -1U) / block; // events cudaEvent_t start = create_event(); cudaEvent_t mid = create_event(); cudaEvent_t stop = create_event(); cudaEventRecord(start); // copy angle[] from host to device cudaMemcpy(dangle, angle, BN, cudaMemcpyHostToDevice); cudaEventRecord(mid); // compute sine kernel_sin<<<grid,block>>>(dsin, dangle, size); // write-back to host cudaMemcpy(sin, dsin, BN, cudaMemcpyDeviceToHost); // compute cosine kernel_cos<<<grid,block>>>(dcos, dangle, size); // write-back to host cudaMemcpy(cos, dcos, BN, cudaMemcpyDeviceToHost); // compute tangent kernel_tan<<<grid,block>>>(dtan, dangle, size); // write-back to host cudaMemcpy(tan, dtan, BN, cudaMemcpyDeviceToHost); cudaEventRecord(stop); cudaDeviceSynchronize(); float elapsed; cudaEventElapsedTime(&elapsed, start, mid); cout << elapsed << " + "; cudaEventElapsedTime(&elapsed, mid, stop); cout << elapsed << " = "; cudaEventElapsedTime(&elapsed, start, stop); cout << elapsed << " [ms]" << endl; for ( unsigned int i = 0; i < 5; ++i ) { cout << angle[i] << " :" << " sin=" << sin[i] << " cos=" << cos[i] << " tan=" << tan[i] << endl; } cout << "..." << endl; for ( auto event : { start, mid, stop }) destroy_event(event); for ( auto host : { angle, sin, cos, tan }) free_host(host); for ( auto device : { dangle, dsin, dcos, dtan }) free_device(device); } int main() { single_stream(); }
なんともオーソドックスなコードです。入力(angle)をホストからデバイスに送り、"kernel_うんちゃら で計算して、デバイスからホストに書き戻す"を三回やってるだけ。
CUDAを使う限り、ホスト・デバイス間のデータ転送は避けられません。そのことが大きなボトルネックになるのは事実ですが、だからってこれ以上速くはならんかというと、そうでもない。
上記コードでは、ホスト・デバイス間のデータ転送をやってる間、GPUは仕事せずデータ転送が終わるのを待ってます。
ピザが焼け、バイトの兄ちゃんがバイク飛ばして配達に行って、帰って来るまでピザ窯空けたままぼーっと待ってるあほゥなピザ職人はいないでしょ? 窯が空いたら次の注文受けて焼いとけば、兄ちゃんが配達から帰って来次第すぐさま次の配達できるやないですか。
...てなオハナシを次のエントリで。
(original: 2015-04-15)
Zero Copy
CUDA C Best Practices Guide 9.1.3 に Zero Copy てのが出てきます...なにコレ? ってんで調べてみました。
サンプル: float列の各要素に対し、その平方根をを求める処理を書きました:
#include <cuda_runtime.h> #include <device_launch_parameters.h> // out[i] = √in[i] where i = 0..size-1 __global__ void kernel_square_root(const float* in, float* out, unsigned int size) { unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { out[i] = sqrtf(in[i]); } } #include <vector> #include <numeric> #include <iostream> void normal() { using namespace std; const unsigned int N = 10000U; // host領域確保 float* in = new float[N]; float* out = new float[N]; // device領域確保 size_t before, after, total; cudaMemGetInfo(&before, &total); cout << "normal : device free = " << before << " -> "; float* din; cudaMalloc(&din, 2*N*sizeof(float)); float* dout; cudaMalloc(&dout, N*sizeof(float)); cudaMemGetInfo(&after, &total); cout << after << " (" << (before-after) << " consumed)" << endl; // in[] = { 0.0, 1.0, 2.0, 3.0 ... } iota(in, in+N, 0.0f); // host -> device cudaMemcpy(din, in, N*sizeof(float), cudaMemcpyHostToDevice); // kernel-call unsigned int block = 256; unsigned int grid = (N + block -1U) / block; kernel_square_root<<<grid,block>>>(din, dout, N); // device -> host cudaMemcpy(out, dout, N*sizeof(float), cudaMemcpyDeviceToHost); for ( unsigned int i : { 2, 3, 200, 300} ) { cout << out[i] << endl; } // device領域解放 cudaFree(din); cudaFree(dout); // host領域解放 delete[] in; delete[] out; }
いつものダンドリです:
一方 zero copy ではデバイス側に領域を確保しません。 そのかわり、デバイスから読み書きできる領域をホスト側に確保します。 zero copyによる「float列の平方根」はこんなコードになります:
void zero_copy() { using namespace std; const unsigned int N = 10000U; // host領域確保 size_t before, after, total; cudaMemGetInfo(&before, &total); cout << "zero-copy: device free = " << before << " -> "; float* in; cudaHostAlloc(&in, N*sizeof(float), cudaHostAllocMapped); float* out; cudaHostAlloc(&out, N*sizeof(float), cudaHostAllocMapped); // hostにマップされたdeviceポインタを取得 float* din; cudaHostGetDevicePointer(&din, (void*)in, 0); float* dout; cudaHostGetDevicePointer(&dout, (void*)out, 0); cudaMemGetInfo(&after, &total); cout << after << " (" << (before-after) << " consumed)" << endl; // in[] = { 0.0, 1.0, 2.0, 3.0 ... } iota(in, in+N, 0.0f); // kernel-call unsigned int block = 256; unsigned int grid = (N + block -1U) / block; kernel_square_root<<<grid,block>>>(din, dout, N); // kernel処理完了を待つ cudaStreamSynchronize(nullptr); for ( unsigned int i : { 2, 3, 200, 300} ) { cout << out[i] << endl; } // host領域解放 cudaFreeHost(in); cudaFreeHost(out); }
違いは3つ:
- ホスト側の領域確保には
cudaHostAlloc( ... cudaHostAllocMapped)
を用いる。 - デバイス側の領域確保は行わず、
cudaHostGetDevicePointer
でホスト側領域に対応するデバイス・ポインタを手に入れ、kernelに渡す。 - kernel呼び出しのあと、kernelの処理が完了するのを待つ(
cudaMemcpy
は処理完了を暗黙裡に待つが、それが行われないため)。
両者を実行しました:
int main() { // zero copy できるか確認 cudaDeviceProp prop; cudaGetDeviceProperties(&prop, 0); if ( !prop.canMapHostMemory ) { std::cerr << "sorry, not supported." << std::endl; return 0; } // zero copy を有効にする cudaSetDeviceFlags(cudaDeviceMapHost); normal(); zero_copy(); }
ご覧のとおり、zero copy版はデバイス・メモリを消費していません。
cudaHostAlloc( ... cudaHostAllocMapped)
で確保されたホスト側領域は物理アドレスが固定(ピン留め)されていて、kernelはこの領域に対し、PCI-bus越しに読み書きするってからくりです。
読み書きが発生するたびにPCI-busを行き来しますから、kernelがこの領域に対して何度も読み書きするとか、(連続領域をリニアに、ではなく)飛び飛びのアドレスにアクセスすると、kernel-callをcudaMemcpy
で挟むいつもの手順に比べてパフォーマンスを著しく損なうことになります。
Tegra K1/X1 なんてな integrated-GPU だとCPUとGPUはひとつのメモリ空間を共有してるんで cudaMemcpy
は無駄であり、ZeroCopyが効果的に機能しますですよきっと。
(original: 2015-08-03)
後付けピン留め
cudaMallocHost
を使ってピン留めされたHost-memoryを確保すればcudaMemcpy
のスキマを潰すことができ、その分速くなる(てかそうしないと複数streamでのoverlapができん)のですが、あらかじめ(malloc
やnew
で)確保された領域を後付けでピン留め...できるんですわこれが。
cudaHostRegister
ってAPIがありまして、領域の先頭アドレスと大きさ(バイト数)を与えることでピン留めしてくれます。ピンを外すのはcudaHostUnregister
で。
#include <cuda_runtime.h> #include <device_launch_parameters.h> #include <math_functions.h> __global__ void kernel_sine(const float* angle, float* sine, unsigned int size) { int i = blockDim.x * blockIdx.x + threadIdx.x; if ( i < size ) { sine[i] = sinf(angle[i]); } } #include <cstdlib> #include <iostream> #define HOSTREGISTER int main() { const unsigned int N = 100000; // Host-memoryを確保 float* h_angle; float* h_sine; h_angle = (float*)malloc(N*sizeof(float)); h_sine = (float*)malloc(N*sizeof(float)); #if defined(HOSTREGISTER) // 後付けピン留め cudaHostRegister(h_angle, N*sizeof(float), cudaHostRegisterDefault); cudaHostRegister(h_sine , N*sizeof(float), cudaHostRegisterDefault); #endif // Device-memoryを確保 float* d_angle; float* d_sine; cudaMalloc(&d_angle, N*sizeof(float)); cudaMalloc(&d_sine , N*sizeof(float)); // HtoD, kernel, DtoH を二回実行 cudaMemcpyAsync(d_angle, h_angle, N*sizeof(float), cudaMemcpyDefault); kernel_sine<<<(N+255)/256,256>>>(d_angle, d_sine, N); cudaMemcpyAsync(h_sine , d_sine , N*sizeof(float), cudaMemcpyDefault); cudaMemcpyAsync(d_angle, h_angle, N*sizeof(float), cudaMemcpyDefault); kernel_sine<<<(N+255)/256,256>>>(d_angle, d_sine, N); cudaMemcpyAsync(h_sine , d_sine , N*sizeof(float), cudaMemcpyDefault); cudaStreamSynchronize(nullptr); // defaul-stream上の処理が完了するまで待機 std::cout << "done." << std::endl; // あとしまつ cudaFree(d_angle); cudaFree(d_sine ); #if defined(HOSTREGISTER) // ピン留め解除 cudaHostUnregister(h_angle); cudaHostUnregister(h_sine ); #endif free(h_angle); free(h_sine ); cudaDeviceReset(); }
timelineはこんな↓、スキマ詰まってますね。
ただね、このcudaHostRegister
使った後付けピン留めってかなり時間がかかるみたい。なので頻繁に留め/外しを繰り返すもんじゃなさそうです。