大数据量计算,请教大家个问题

我现在有个2000x2000的矩阵,想要对所有元素求和,该怎么做效率最高呢?谢谢!
我试过先处理每行,再求所有列的和,但感觉效率不理想。

而且我的卡是GTX650,只能支持1024个线程,无法一次处理一行

楼主您好,您的卡同时支持更多线程的。您说的1024,是一个block里的线程最多数目。但您可以使用更多blocks的。

例如一次启动<<<>>>中的第一个参数,则可以指定您需要blocks数目。

其二,如果求和。
楼主您这个问题是典型的卡在访存上的,直接求和即可,举个例子说:
//array是您的矩阵,result是您要求的和。简单示范而已。
global void sum(int *array, int *result)
{
int id = blockIdx.x * blockDim.x + threadIdx.x;
int step = gridDim.x * blockDim.x;
int sum = 0;
while (id < 4000000)
{
sum += array[id];
}
atomicAdd(result, sum);
}

然后<<<128,256>>>之类的启动即可。
你的算法必然卡在贵卡的访存上,只要能达到访存最高效率即可。无需可以优化什么“计算”。

以及,需要指出的是,如果数据原本是在CPU memory上的,建议就地求和。不要移动到Global memory的。除非数据已经在global memory上,这样才有意义。

因为CPU求和一遍的时间,也只有读取全部数据一遍的时间。
而您如果复制了,可能需要2倍读取的时间+PCI-E传输时间+kernel运行时间。不一定合算的。

非常感谢版主的答复:
1、我的这个矩阵本身就是在global memory里计算出来的,并没有跟CPU memory之间交换,所以不存在重复的传输。
2、我知道可以用多个block,可能昨晚发帖没讲清楚(手机发的,呵呵),我现在的做法是每个block分为1000个thread,对1000个元素归约求和,然后再把每个block的结果归约一遍。

版主说的问题,我会仔细考虑并测试一下,有结果了就上来通报,谢谢!

您客气了。

另外,出于精度要求,我的数组元素都是double型的,怎么原子加呢?

double建议如下操作:

(1)每个线程的执行多次加法,然后block内部对线程们的结果进行规约,然后每个block的结果写入global memory.
(2)再执行个kernel, 对(1)中的这些结果规约。

我用的就是这种方法,只是想知道有没有更高效的方法,代码如下(出于项目原因,删除了一些内部处理,但大体流程可以看出来):
pdInRefPoint输入矩阵
pdOutK输出矩阵,存放每个block的求和结果

__global__ void test_kernel(double2 *pdInRefPoint, size_t stInRefPointPitch, 
   double *pdOutK, size_t stOutRefPointPitch, 
   _DWORD dwWidth, _DWORD dwHeight)
{
   //shared memory首地址
   extern __shared__ __align__(8) uchar4 pstSharedMemoryDev[];

   double2 *adInRefPoint = (double2 *)pstSharedMemoryDev;
   double3 *adKWithZ1 = (double3 *)&adInRefPoint[dwWidth >> 1];
   double3 *adKWithZ2 = &adKWithZ1[dwWidth >> 2];

   _DWORD dwXIndex = threadIdx.x;
   _DWORD dwYIndex = blockIdx.x;
   _DWORD i;

   _DWORD dwIndex = dwYIndex * ((stInRefPointPitch / sizeof(_DOUBLE)) / sizeof(uchar2)) + dwXIndex;

   //把参考点搬移到shared memeory
   adInRefPoint[dwXIndex] = pdInRefPoint[dwIndex];
   __syncthreads();

   //以下计算拟合方程所需要用的矩阵中与参考点高度有关的数据
   if ( 0 == (dwXIndex & 1) ) //能被2整除
   {
   _DOUBLE *pdRefPoint = &adInRefPoint[dwXIndex].x;

   adKWithZ1[dwXIndex >> 1].x = 0;
   adKWithZ1[dwXIndex >> 1].y = 0;
   adKWithZ1[dwXIndex >> 1].z = 0;
   adKWithZ2[dwXIndex >> 1].x = 0;
   adKWithZ2[dwXIndex >> 1].y = 0;
   adKWithZ2[dwXIndex >> 1].z = 0;

   for ( i = 0; i < 4; i++, pdRefPoint++ ) //每个线程计算四个加法
   {
   if ( INVALID_COMPENSATE_HEIGHT != *pdRefPoint )
   {
   dwIndex = dwXIndex * sizeof(uchar2) + i;

   adKWithZ1[dwXIndex >> 1].x += (_DOUBLE)dwIndex * (_DOUBLE)dwIndex * (*pdRefPoint);
   adKWithZ1[dwXIndex >> 1].y += (_DOUBLE)dwYIndex * (_DOUBLE)dwYIndex * (*pdRefPoint);
   adKWithZ1[dwXIndex >> 1].z += (_DOUBLE)dwIndex * (_DOUBLE)dwYIndex * (*pdRefPoint);
   adKWithZ2[dwXIndex >> 1].x += (_DOUBLE)dwIndex * (*pdRefPoint);
   adKWithZ2[dwXIndex >> 1].y += (_DOUBLE)dwYIndex * (*pdRefPoint);
   adKWithZ2[dwXIndex >> 1].z += (*pdRefPoint);
   }
   }
   }
   __syncthreads();

   //转换为2的幂,便于归约
   for ( i= (m_adwPowerTable[blockDim.x] >> 2); i > 0; i >>= 1)
   {
   if ( dwXIndex < i && dwXIndex + i < (dwWidth >> 2) )
   {
   adKWithZ1[dwXIndex].x += adKWithZ1[dwXIndex + i].x;
   adKWithZ1[dwXIndex].y += adKWithZ1[dwXIndex + i].y;
   adKWithZ1[dwXIndex].z += adKWithZ1[dwXIndex + i].z;
   adKWithZ2[dwXIndex].x += adKWithZ2[dwXIndex + i].x;
   adKWithZ2[dwXIndex].y += adKWithZ2[dwXIndex + i].y;
   adKWithZ2[dwXIndex].z += adKWithZ2[dwXIndex + i].z;
   }

   __syncthreads();
   }

   if ( 0 == dwXIndex )
   {
   pdOutK[CUDA_COMPENSATE_MATRIX_LENGTH * dwYIndex + CUDA_COMPENSATE_K_WITHOUT_Z] = adKWithZ1[dwXIndex].x;
   pdOutK[CUDA_COMPENSATE_MATRIX_LENGTH * dwYIndex + CUDA_COMPENSATE_K_WITHOUT_Z + 1] = adKWithZ1[dwXIndex].y;
   pdOutK[CUDA_COMPENSATE_MATRIX_LENGTH * dwYIndex + CUDA_COMPENSATE_K_WITHOUT_Z + 2] = adKWithZ1[dwXIndex].z;
   pdOutK[CUDA_COMPENSATE_MATRIX_LENGTH * dwYIndex + CUDA_COMPENSATE_K_WITHOUT_Z + 3] = adKWithZ2[dwXIndex].x;
   pdOutK[CUDA_COMPENSATE_MATRIX_LENGTH * dwYIndex + CUDA_COMPENSATE_K_WITHOUT_Z + 4] = adKWithZ2[dwXIndex].y;
   pdOutK[CUDA_COMPENSATE_MATRIX_LENGTH * dwYIndex + CUDA_COMPENSATE_K_WITHOUT_Z + 5] = adKWithZ2[dwXIndex].z;
   }
}

然后再用另一个kernel对上面的结果再次归约求和。
第二个kernel时间很短,0.0x个ms而已,而第一个kernel,需要好几个ms(大约6、7个ms),怎么相差如此之大,我就是觉得效率应该还有提升空间,所以想请教下大家,谢谢!

首先说,楼主您使用的不是“我这种方法”。

我们的主要区别在于,举个例子说,有100M个数据和256个blocks,每个里面256个线程。
我的算法在上述数据的情况下是:
kernel 1:
(1)每个block里的每个线程读取并累加1600个数据。(在寄存器中!)
(2)然后整个block对256个线程进行规约8次,得到1个值。写入global memory.
kernel 2:
(3)对256个blocks产生的256个数据进行规约,得到1个最终数据。

而不是,kernel 1尝试对大量数据直接规约,反复折磨shared memory。(请注意我刚才的步骤1是在寄存器中的)。

这可能是您“感觉慢“的原因之一。

其二,我们知道,累加是个典型的访存密集型应用,无论怎么写,随便写总是容易到峰值的(卡在访存上)。您可能只是心里作用。贵kernel的几个ms可能对于您的数据规模,已经是峰值了。
则此时您无需修改。:slight_smile:

(此外,我看到了大量的DWORD类似字样,您以前是win32开发员吗?)

以及,我的上文只是建议,不排除您有其他地方的问题(也不排除您完全无问题)。

您可以从如下方面估算您大致需要的时间,
假设数据类型是double, 贵卡是2Tflopsd float能力/1TFlops double能力(FMA), 访存global memory是在100GB/S,数据规模有100M个。

那么您需要的理论计算:100M / 0.5 T次/s = 0.2 ms
那么您需要的理论访存:100M * sizeof(double ) / 100GB/s = 8ms
取瓶颈8ms。

然后您可以根据这个,折算您的卡的具体参数,和您的数据规模,得到理论最小时间。
然后根据理论最小时间估算您是否需要修改算法或者进行优化。

[
多谢,我还是有点不太明白,怎么在寄存器里对若干个double进行累加,并对block可见呢?请指点。
关于性能,我算算看。
恰好相反,我一直是做嵌入式和linux方向,近两年才开始做windows开发,一周前开始接触CUDA,只是这么多年保留下来的编码规范而已,呵呵。

先每个线程单独累加是为了避免上去就进行昂贵的规约操作。
如何累加?您可以采用多种加法方式。最简单的就是反复+=。

(例如1个block有256个线程,每个线程进行1000次加法。这样25600个数据将得到256个和)

然后您可以使用shared memory, 进行您惯常的(例如您上文)中的规约。

最后进行总规约。

这样可以极大的减少规约操作的数目。

建议先上去单独加法是因为:
(1)我们的N卡不能就地在shared memory里将2个数据加起来的,需要读取2个数-加法-回写的。大量这么来代价高昂。
(2)而单独加法可以有效的避免这个过程。

以及,DWORD之所以叫这个名字是因为很多原因。

任何在非win32 api交互的场合使用该类型,都是不推荐的。
我建议您使用int32_t之类的标准化的名字,而不是DWORD这个带有太多历史负担的名字。

使用DWORD可能是不好的“规范”。

以上均是建议。而且和CUDA无关。仅供参考。

明白了,我一开始理解错您的回复,把问题想复杂了。。。

我把这里稍微说明下,虽然估计LZ已经明白是怎么回事了,不过为了看帖的其他人能快速明白,还是写下吧。

这里其实不需要维护对block可见,各个线程自己加自己的局部和即可,其结果对自己可见就足够了。
当各个线程各自累加部分完成的时候,此时才需要线程间合作进行常见的那种规约,此时才需要维护对block内部其他线程可见,而block内可见的有global memory和shared memory,一般可以考虑后者,速度高延迟短。

大致这样,横扫斑竹已经完整地回答了LZ的问题,我在这稍微补充下,以方便看帖的人。

祝LZ编码顺利~