CUDAによるIPCA
前にC++で動くIPCAのソースを載せましたが、今度はCUDAで動くIPCAのソースを載せます。
実行速度はというと、まあ早い。
CPUにやらせるより、10倍よりも早くはならないといった感じでしょうか。
/* * param[in] frame フレーム数 * param[in] arrIn 入力画像 HEIGHT*STRIDE * param[in] tableV 主成分画像の配列を格納 HEIGHT*STRIDE*DIMENSION * param[in] tableU 計算用画像の配列を格納 HEIGHT*STRIDE*DIMENSION * HEIGHT = THREAD としています */ __global__ void ipca_kernel( int frame, float* arrIn, float* tableV, float* tableU ) { extern __shared__ float sharedmem[]; float* amtV = (float*)sharedmem; float* amtUV = (float*)sharedmem + THREAD; ///// スレッドIDを取得 const unsigned int tid = threadIdx.x; ///// arrIn -> Sub[0] for(int e=0; e<STRIDE; e++) { tableU[STRIDE*tid + e] = arrIn[STRIDE*tid + e]; } __syncthreads(); int iterate = DIMENSION - 1; if( DIMENSION > frame ) iterate = frame; for(int d=0; d<=iterate; d++) { float* arrVd = tableV + THREAD*STRIDE*d;// Main float* arrUd = tableU + THREAD*STRIDE*d;// Sub //float* arrVd1 = tableV + THREAD*STRIDE*(d+1); float* arrUd1 = tableU + THREAD*STRIDE*(d+1); if(frame == d) { for(int e=0; e<STRIDE; e++) arrVd[tid*STRIDE + e] = arrUd[tid*STRIDE + e]; __syncthreads(); break; } amtV[tid] = 0.0f; amtUV[tid] = 0.0f; for(int e=0; e<STRIDE; e++) { float v = arrVd[tid*STRIDE + e]; float u = arrUd[tid*STRIDE + e]; amtV[tid] += v * v; amtUV[tid] += u * v; } __syncthreads(); if(tid == 0) { float nrmV = 0.0f; float dotUV = 0.0f; for(int t=0; t<THREAD; t++) { nrmV += amtV[t]; dotUV += amtUV[t]; } amtV[0] = (float)sqrt(nrmV); amtUV[0] = dotUV; } __syncthreads(); float scalerA = ((float)frame - 1.0) / frame; float scalerB = amtUV[0] / ((float)frame * amtV[0]); for(int e=0; e<STRIDE; e++) { arrVd[tid*STRIDE + e] = arrVd[tid*STRIDE + e] * scalerA + arrUd[tid*STRIDE + e] * scalerB; } __syncthreads(); if( d >= iterate ) continue; ///// dotUV, nrmV の再計算 amtV[tid] = 0.0; amtUV[tid] = 0.0; for(int e=0; e<STRIDE; e++) { float v = arrVd[tid*STRIDE + e]; float u = arrUd[tid*STRIDE + e]; amtV[tid] += v * v; amtUV[tid] += u * v; } __syncthreads(); if(tid == 0) { float nrmV = 0; float dotUV = 0; for(int t=0; t<THREAD; t++) { nrmV += amtV[t]; dotUV += amtUV[t]; } amtV[0] = sqrt(nrmV); amtUV[0] = dotUV; } __syncthreads(); float scalerC = amtUV[0] / (amtV[0]*amtV[0]); for(int e=0; e<STRIDE; e++) { arrUd1[tid*STRIDE + e] = arrUd[tid*STRIDE + e] - arrVd[tid*STRIDE + e] * scalerC; } __syncthreads(); } }