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