CUDA新手入门，一个简单矩阵问题

[shinson2008]

今天CUDA技术群里meteor兄提了个问题如下

引用

x x x x y y y y
x x x x y y y y
x x x x y y y y
x x x x y y y y
z z z z a a a a
z z z z a a a a
z z z z a a a a
z z z z a a a a

比如这个矩阵
我想对x那些数加一,对y那些数加二
对z那些书加三,对a那些数加四

介于meteor兄是新手，本着互相学习，下面是我写的一段很很简单的程序，完成上述操作，希望对meteor兄有所帮助

引用

#include <stdio.h>

__global__ void testkernel(int *d_A, size_t size)
{
    int dx = blockDim.x * blockIdx.x + threadIdx.x;
    int dy = blockDim.y * blockIdx.y + threadIdx.y;

    if( blockIdx.x == 0 && blockIdx.y == 0 )
       d_A[dx*size+dy] += 1;
    if( blockIdx.x == 0 && blockIdx.y == 1 )
       d_A[dx*size+dy] += 2;
    if( blockIdx.x == 1 && blockIdx.y == 0 )
       d_A[dx*size+dy] += 3;
    if( blockIdx.x == 1 && blockIdx.y == 1 )
       d_A[dx*size+dy] += 4;
}

int main( int argc, char** argv) 
{
int h_A[8][8] = {{1,1,1,1,2,2,2,2},
                 {1,1,1,1,2,2,2,2},
                 {1,1,1,1,2,2,2,2},
                 {1,1,1,1,2,2,2,2},
                 {3,3,3,3,4,4,4,4},
                 {3,3,3,3,4,4,4,4},
                 {3,3,3,3,4,4,4,4},
                 {3,3,3,3,4,4,4,4}};
 
int  *d_A, *h_B;
size_t size = 8 * 8 * sizeof(int);
size_t rsize = 8;
dim3 dimgrid(2,2);
dim3 dimblock(4,4);

h_B = (int*)malloc(size);

cudaMalloc( (void **) &d_A, size );
cudaMemcpy( d_A, h_A, size, cudaMemcpyHostToDevice );

testkernel<<<dimgrid,dimblock>>>(d_A,rsize);

cudaMemcpy( h_B, d_A, size, cudaMemcpyDeviceToHost );

for(int i = 0; i < 8; i++)
{
  for(int j = 0;j < 8; j++)
      printf("%2d ",h_B[i*rsize+j]);
printf("\n");
}

cudaFree(d_A);
free(h_B);
}

介于meteor兄不理解blockDim.x和threadIdx.x，下面借上面这个例子解释，具体的请参见Programme Guide
blockDim就是指block的维度,这里每个block是4*4的,所以blockDim.x=4 blockDim.y = 4
threadIdx就是指block里的线程的索引号,这里每block是4*4维的,每个block里有16个thread,每个thread的threadIdx.x从0到3,threadIdx.y从0到3,和数组一样,这样解释行吗?

以上程序测试通过。。。。

[OpenHero]

[Quote=引用 4 楼 FallingStar08 的回复:]
请问你们群还收人吗
[/Quote]收

[chenglouguanshanjing]

历害

[Muscle000]

直接加四个稀疏矩阵不就行了吗？

[riftvalley]

不好意思，是不是我C语言也没学好，我怎么感觉
if( blockIdx.x == 0 && blockIdx.y == 0 )
d_A[dx*size+dy] += 1;
应该改为
if( blockIdx.x == 0 && blockIdx.y == 0 )
d_A[dx+dy*size] += 1;

[riftvalley]

哦，我想错了，坐标思想太严重

[NaWait]

__global__ void testkernel(int *d_A, size_t size)
{
    int dx = blockDim.x * blockIdx.x + threadIdx.x;
    int dy = blockDim.y * blockIdx.y + threadIdx.y;

    if( blockIdx.x == 0 && blockIdx.y == 0 )
       d_A[dx*size+dy] += 1;
    if( blockIdx.x == 0 && blockIdx.y == 1 )
       d_A[dx*size+dy] += 2;
    if( blockIdx.x == 1 && blockIdx.y == 0 )
       d_A[dx*size+dy] += 3;
    if( blockIdx.x == 1 && blockIdx.y == 1 )
       d_A[dx*size+dy] += 4;
}

if 用這麼多 效率應該很差吧

不如用矩陣加法加一加快多了

[porscher]

[Quote=引用 6 楼 OpenHero 的回复:]
引用 4 楼 FallingStar08 的回复:
请问你们群还收人吗
收
[/Quote]

你们的群收新人么，我刚开始接触CUDA的
群号是多少啊？

[电力信息系统]

运行结果不对啊？
1 1 1 1 2 2 2 2
1 1 1 1 2 2 2 2
1 1 1 1 2 2 2 2
1 1 1 1 2 2 2 2
3 3 3 3 4 4 4 4
3 3 3 3 4 4 4 4
3 3 3 3 4 4 4 4
3 3 3 3 4 4 4 4
3868016 3868016 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
请按任意键继续. . .

[ziyuanyatou]

为什么在我的电脑上运行总是有错误啊，新手上路，麻烦多多关照！

[seagullor]

假设lz给的矩阵没有规律性，即：
int h_A[8][8] = {{1,1,1,1,2,2,2,2},
{1,1,1,1,2,2,2,2},
{1,1,1,1,2,2,2,2},
{1,1,1,1,2,2,2,2},
{3,3,3,3,4,4,4,4},
{3,3,3,3,4,4,4,4},
{3,3,3,3,4,4,4,4},
{3,3,3,3,4,4,4,4}};
这其中1，2，3，4的位置和个数都是混乱的，不如上面那么规划整一
是不是就必须得要if来判断当前元素的多少啊

[bewarm]

__global__ void testkernel(int *d_A, size_t size)
{
	int dx = blockDim.x * blockIdx.x + threadIdx.x;
	int dy = blockDim.y * blockIdx.y + threadIdx.y;

	if( blockIdx.x == 0 && blockIdx.y == 0 )
		d_A[dx*size+dy] += 1;
	if( blockIdx.x == 0 && blockIdx.y == 1 )
		d_A[dx*size+dy] += 2;
	if( blockIdx.x == 1 && blockIdx.y == 0 )
		d_A[dx*size+dy] += 3;
	if( blockIdx.x == 1 && blockIdx.y == 1 )
		d_A[dx*size+dy] += 4;
}


/************************************************************************/
/* HelloCUDA                                                            */
/************************************************************************/
int main(void)
{
	int h_A[8][8] = {{1,1,1,1,2,2,2,2},
	{1,1,1,1,2,2,2,2},
	{1,1,1,1,2,2,2,2},
	{1,1,1,1,2,2,2,2},
	{3,3,3,3,4,4,4,4},
	{3,3,3,3,4,4,4,4},
	{3,3,3,3,4,4,4,4},
	{3,3,3,3,4,4,4,4}};

	int  *d_A, *h_B;
	size_t size =  8 * sizeof(int);
	size_t pitch =  sizeof(int);
	size_t rsize = 8;
	dim3 dimgrid(2,2);
	dim3 dimblock(4,4);

	h_B = (int*)malloc(size);

	//cudaMalloc( (void **) &d_A, size );
	cudaMallocPitch( (void **) &d_A, &pitch,size,rsize );
	//cudaMemcpy( d_A, h_A, size, cudaMemcpyHostToDevice );
	cudaMemcpy2D( d_A, pitch,h_A,pitch, size, rsize,cudaMemcpyHostToDevice );

	testkernel<<<dimgrid,dimblock>>>(d_A,rsize);

	//cudaMemcpy( h_B, d_A, size, cudaMemcpyDeviceToHost );
	cudaMemcpy2D( h_B, pitch,d_A,pitch, size,rsize, cudaMemcpyDeviceToHost );

	for(int i = 0; i < 8; i++)
	{
		for(int j = 0;j < 8; j++)
			printf("%2d ",h_B[i*rsize+j]);
		printf("\n");
	}

	cudaFree(d_A);
	free(h_B);

我想把代码优化一下用了一个cudaMallocPitch和cudaMemcpy2D
然后调试的时候，中间会中断一下
这是VC给的信息：HEAP[CUDAWinApp2.exe]: HEAP: Free Heap block be1800 modified at be1810 after it was freed
出来了这样结果：


是在看不出哪里不对啊

[bewarm]

[Quote=引用 16 楼 bewarm 的回复:]
C/C++ code__global__ void testkernel(int *d_A, size_t size)
{
int dx = blockDim.x * blockIdx.x + threadIdx.x;
int dy = blockDim.y * blockIdx.y + threadIdx.y;

if( blockIdx.x == 0 && blockIdx.y == 0 )
d_A[dx*size+dy] += 1;
if( blockIdx.x == 0 && blockIdx.y == 1 )
d_A[dx*size+dy] += 2;
if( blockIdx.x == 1 && blockIdx.y == 0 )
d_A[dx*size+dy] +=…
[/Quote]
解决了

[Cyrosly]

哈哈.这么多人,又这么简单的问题哎.......,要学会利用硬件架构的知识.否则...~~

--完全可以消除转移指令

原排列:
form<0>
x  x  x  x .. x  y  y  y  y .. y

x  x  x  x .. x  y  y  y  y .. y

x  x  x  x .. x  y  y  y  y .. y

x  x  x  x .. x  y  y  y  y .. y 

.  .  .  . .. .  .  .  .  . .. .

x  x  x  x .. x  y  y  y  y .. y

z  z  z  z .. z  a  a  a  a .. a

z  z  z  z .. z  a  a  a  a .. a

z  z  z  z .. z  a  a  a  a .. a

z  z  z  z .. z  a  a  a  a .. a

.  .  .  . .. .  .  .  .  . .. .

z  z  z  z .. z  a  a  a  a .. a

还可以变化为下面2个排列:

form<1>
x  x  x  x .. x  y  y  y  y .. y  z  z  z  z .. z  a  a  a  a .. a

x  x  x  x .. x  y  y  y  y .. y  z  z  z  z .. z  a  a  a  a .. a

x  x  x  x .. x  y  y  y  y .. y  z  z  z  z .. z  a  a  a  a .. a

x  x  x  x .. x  y  y  y  y .. y  z  z  z  z .. z  a  a  a  a .. a

.  .  .  . .. .  .  .  .  . .. .  .  .  .  . .. .  .  .  .  . .. .

x  x  x  x .. x  y  y  y  y .. y  z  z  z  z .. z  a  a  a  a .. a

form<2>
x  x  x  x  x  x  x  x  x  x  x  x  x  x  x .. x
.  .  .  .  .  .  .  .  .  .  .  .  .  .  . .. . 
x  x  .  .  .  .  .  .  .  .  .  .  .  .  x .. x

y  y  y  y  y  y  y  y  y  y  y  y  y  y  y .. y
.  .  .  .  .  .  .  .  .  .  .  .  .  .  . .. . 
y  y  .  .  .  .  .  .  .  .  .  .  .  .  y .. y

z  z  z  z  z  z  z  z  z  z  z  z  z  z  z .. z
.  .  .  .  .  .  .  .  .  .  .  .  .  .  . .. .
z  z  .  .  .  .  .  .  .  .  .  .  .  .  z .. z

a  a  a  a  a  a  a  a  a  a  a  a  a  a  a .. a
.  .  .  .  .  .  .  .  .  .  .  .  .  .  . .. .
a  a  .  .  .  .  .  .  .  .  .  .  .  .  a .. a


如果这4个区域每个区域中的每个数都相同(那作者的意思是?),则可以如下实施:

//const_cache[0]=x+1
//const_cache[1]=y+2
//const_cache[2]=z+3
//const_cache[3]=a+4
__constant__ unsigned int const_cache[4];

对于form<0>:
__global__ void kernel_derisible_form_0(unsigned int* dst)
{
    unsigned int tidx =__umul24(blockDim.x,blockIdx.x)+threadIdx.x;
    unsigned int tidy =__umul24(blockDim.y,blockIdx.y)+threadIdx.y;
    unsigned int sizex=__umul24(gridDim.x,blockDim.x);
    unsigned int gloc =__umul24(size,tidy)+tidx;
    sizex>>=1;
    unsigned int sizey=__umul24(gridDim.y,blockDim.y);
    sizey>>=1;

    dst[gloc]=const_cache[((tidy/sizey)<<1)+(tidx/sizex)];
}
     

对于form<1>:
__global__ void kernel_derisible_form_1(unsigned int* dst)
{
    unsigned int tidx=__umul24(blockDim.x,blockIdx.x)+threadIdx.x;
    unsigned int tidy=__umul24(blockDim.y,blockIdx.y)+threadIdx.y;
    unsigned int size=__umul24(gridDim.x,blockDim.x);
    unsigned int gloc=__umul24(size,tidy)+tidx;

    size>>=2;

    dst[gloc]=const_cache[tidx/size]; //broadcast to each thread which the should for
}


对于form<2>:
__global__ void kernel_derisible_form_2(unsigned int* dst)
{
    unsigned int tidx=__umul24(blockDim.x,blockIdx.x)+threadIdx.x;
    unsigned int tidy=__umul24(blockDim.y,blockIdx.y)+threadIdx.y;
    unsigned int gloc=__umul24(__umul24(gridDim.x,blockDim.x),tidy)+tidx;

    unsigned int cloc=threadIdx.y/(__umul24(gridDim.y,blockDim.y)>>2);

    dst[gloc]=const_cache[cloc]; //same the up annotation
}

如果每个区域中的每个元素不一定相同,则可如下处理(事实上更通用了,但对于上面说的特殊情况,上面利用const cache的广播机制的方法可以更高效些):
__global__ void kernel_derisible_form_0_general(unsigned int* dst)
{
    unsigned int tidx=__umul24(blockDim.x,blockIdx.x)+threadIdx.x;
    unsigned int tidy=__umul24(blockDim.y,blockIdx.y)+threadIdx.y;
    unsigned int sizex=__umul24(gridDim.x,blockDim.x);
    unsigned int gloc=__umul24(size,tidy)+tidx;
    sizex>>=1;
    unsigned int sizey=__umul24(gridDim.y,blockDim.y);
    sizey>>=1;

    dst[gloc]+=((tidy/sizey)<<1)+(tidx/sizex); //注意这里和form<0>的代码中最后的const_cache的索引是一样的
}

最后需要置仪的是half warp的一致性,我所说的一致性是指half warp内的所有thread应该尽可能的访问同一个const cache的address,否则会有bank confilicts.这就看选择哪中数据分布以及数据的数量了(比如要求是16的倍数,但是如果矩阵足够的小(如:矩阵中的4个区域分布在同一个half warp内,那么就另当别论了(那就尽量减少BC,但是这样的话也没有必要使用CUDA了)

[101monster]

不错...

[yanghangjun]

学习学习

[yaoyaoflyfreely]

群号多少啊？？

[lessmin]

我想知道 如果把testkernel 单独写个.cu文件 该怎么写 是写在 源文件 资源文件 还是头文件里 为什么我不管放在那里都是错误的

[lessmin]

我想知道 如果把testkernel 单独写个.cu文件 该怎么写 是写在 源文件 资源文件 还是头文件里 为什么我不管放在那里都是错误的

[chjunhua]

mark

[goodbeybey]

学习了~~~