结果为什么不对呀？

system · 2012 年11 月 22 日 09:14

#define XROW 200
#define XCOL 100 
__global__ void KernelV(int *dev_XMatrix,double *dev_V,double *dev_VK)
{
int tid=threadIdx.x;
int bid=blockIdx.x;
int xIndex=bid*blockDim.x+tid;
   for(int i=0;i<XROW;i++)
  {
   for(int k=0;k<XROW;k++)
   {
   dev_V[i*XROW+k]=dev_V[i*XROW+k]+dev_XMatrix[xIndex+i*XCOL]*dev_XMatrix[xIndex+k*XCOL]*dev_VK[xIndex];
   }
  }
  }

int main(void)
{
clock_t t1,t2;
int *XMatrix=(int*)malloc(XROW*XCOL*sizeof(int));
double *V=(double*)malloc(XROW*XROW*sizeof(double));
double *VK=(double*)malloc(XCOL*sizeof(double));
double *V2=(double*)malloc(XROW*XROW*sizeof(double));
for(int i=0;i<XCOL;i++)
{
  VK=2.0;
}
for(int i=0;i<XROW*XCOL;i++)
{
  XMatrix=1.0;
}
for(int i=0;i<XROW*XROW;i++)
{
  V=0.0;
}
int *dev_XMatrix;
double *dev_V,*dev_VK;
cudaMalloc((void**)&dev_XMatrix,XROW*XCOL*sizeof(int));
cudaMalloc((void**)&dev_V,XROW*XROW*sizeof(double));
cudaMalloc((void**)&dev_VK,XCOL*sizeof(double));
cudaMemcpy(dev_XMatrix,XMatrix,XROW*XCOL*sizeof(int),cudaMemcpyHostToDevice);
cudaMemcpy(dev_V,V,XROW*XROW*sizeof(double),cudaMemcpyHostToDevice);
cudaMemcpy(dev_VK,VK,XCOL*sizeof(double),cudaMemcpyHostToDevice);
   t1=clock();

KernelV<<<1,100>>>(dev_XMatrix,dev_V,dev_VK);
t2=clock();
cudaMemcpy(V2,dev_V,XROW*XROW*sizeof(double),cudaMemcpyDeviceToHost);
for(int i=0;i<10;i++)
{
  for(int j=0;j<6;j++)
  {
   printf("%f ",V2[i*6+j]);
  }
  printf("\n");
}
printf("CPU所用时间：%d毫秒\n",(t2-t1));
free(V);
free(XMatrix);
free(VK);
free(V2);
}

我的目的是将以下代码改为并行:

for(int j=0;j<XCOL;j++)
{
  for(int i=0;i<XROW;i++)
  {
   for(int k=0;k<XROW;k++)
   {
   V[i*XROW+k]=V[i*XROW+k]+XMatrix[j+i*XCOL]*XMatrix[j+k*XCOL]*VK[j];
   }
  }
}

我是想j个元素同时运算，但是运行结果是错误的~

system · 2012 年11 月 22 日 09:57

大致看了一下您的代码，您是打算将最后给出的那个三重循环的最外面那一重拆开，每个thread完成一个j对应的计算。

我觉得这里面可能有两个问题
1：您kernel中的赋值语句是给“dev_V[i*XROW+k]”赋值，您所有的100个线程都在给同样的一个数组赋值（并不保证赋值的顺序），这一行为可能和您原先串行的实现不等价，会产生预料之外的结果。
2：您使用了双精度浮点数，这需要您编译时设定计算能力为1.3或更高。请也检查这一点。

一个可能的建议：
如果您不是将j的循环拆开，而是拆开i，k的二重循环，或许是合适的方案。

感谢您莅临cudazone，祝您编码愉快！

system · 2012 年11 月 22 日 11:02

我的显卡是GTX580，计算能力是2.0，支持双精度是没有问题的。好的，我试试看。

system · 2012 年11 月 22 日 11:39

580是计算能力2.0的卡，但依然需要您正确配置。

祝您编码愉快。

system · 2012 年11 月 22 日 11:54

修改为如下代码，可以得到正确的结果，XROW=600，XCOL=200时进行测试，结果显示，串行用了480ms，并行用了1ms(包括了数据传输)。但是又出现了新的问题，当我将XROW=2000，XCOL=1000时，出现了显卡驱动重启，GPU透过WDDM进行超时检测和恢复，我按照网上的说法修改了注册表，还是无效啊！！！
global void KernelV(int *dev_XMatrix,double *dev_V,double dev_VK)
{
int tid=threadIdx.x;
int bid=blockIdx.x;
int xIndex=bidblockDim.x+tid;
for(int j=0;j<XCOL;j++)
{
for(int i=0;i<XROW;i++)
{

		dev_V[i*XROW+xIndex]=dev_V[i*XROW+xIndex]+dev_XMatrix[j+i*XCOL]*dev_XMatrix[j+xIndex*XCOL]*dev_VK[j];
		
	}
}

}

system · 2012 年11 月 22 日 12:50

现在可以了，我将注册表修改了，将WIN7超时检测功能关闭了，不知道会不会对其他功能有影响~

system · 2012 年11 月 22 日 13:10

看上去不错的样子，恭喜恭喜，一般无不良影响，另外修改超时检测后要重启下系统才能生效。

system · 2012 年11 月 22 日 13:22

但是我有个问题一直不明白，当我的XROW=2000，XCOL=1000的时候，就会黑屏一下，显卡驱动重启，而用CUBLAS的矩阵相乘，60006000以上才会出现这种状况，用CULA的矩阵求逆1200012000才会出现这种状况，这是为什么啊？

system · 2012 年11 月 22 日 13:23

另外这么简单的代码有必要用share memory吗？希望给解释一下。

system · 2012 年11 月 22 日 13:38

可能因为你的代码效率低，人家60006000才超时，而你的代码20001000就到时间了。

system · 2012 年11 月 22 日 13:45

我这里十分卡顿了，两楼一起简单说下：

1：这可能是因为同等规模下，您的代码比CUBLAS的实现工作量多，耗时长。或者CUBLAS可额能进行了某种分块，致使一个kernel运行时间较短。不过您的工作已经胜利完成，这些不是很重要了。

2：如果您需要block内通信，或者可以读进来反复用以减少访存压力，那么可以考虑shared memory。否则可以不管。

再次恭喜您解开“BUG退散，风轻云淡”成就，祝您编码愉快！

system · 2012 年11 月 22 日 14:09

好的。提醒了我当block内需要线程通信的时候要用到share memory.O(∩_∩)O谢谢

system · 2012 年11 月 22 日 14:31

善哉，EXP +5！

祝您编码愉快！

system · 2012 年11 月 23 日 15:43

我发现我时间函数弄错了，我写的kernel并没有这么快，10050花了98ms，而串行才用了68ms;600100并行才比串行快，好像是137ms,串行用了200ms。这个for循环还能不能再优化一下呢？我就是那个发了一段SAS程序问能否并行的同学，版主跟我说整体不能并行，只能细类度的并行，我的细类度并行大都如刚才这个for循环，求高手优化指点~

system · 2012 年11 月 24 日 03:05

1:确实需要较大规模才能体现出对CPU的优势，小规模问题CPU更快。

2：结合您顶楼的代码和后面给出的并行实现，可以看出，您只将最内层的k循环进行了展开，如果将i循环也进行展开，可能效果会更好一些，此时可以安排更多的线程同时为您服务。

祝您编码愉快~

system · 2012 年11 月 25 日 02:35

修改完了，这次的确加速明显了，但是第一个数据dev_V[0]的值是正确值的2倍，其他值都是正确的，这是什么原因？
global void KernelV(int *dev_XMatrix,double *dev_V,double dev_VK)
{
int tid=threadIdx.x;
int bid=blockIdx.x;
int xIndex=bidblockDim.x+tid;

if(xIndex<XROW)
{

	for(int j=0;j<XCOL;j++)
	{
		dev_V[xIndex]=dev_V[xIndex]+dev_XMatrix[j]*dev_XMatrix[j+xIndex*XCOL]*dev_VK[j];
		__syncthreads();
		dev_V[xIndex*XROW]=dev_V[xIndex*XROW]+dev_XMatrix[j+xIndex*XCOL]*dev_XMatrix[j]*dev_VK[j];
	}
	__syncthreads();
}

}

system · 2012 年11 月 25 日 04:13

根据您的代码，当threadIdx.x==0 且 blockIdx.x==0时，您的xIndex==0。

此时j循环内的两行赋值总是在干同样的事情，所以会导致第一个值变为2倍。

如果根据您的算法，这样仅仅是第一个值变为2倍，而没有其他错误的话（这一点需要您自己根据算法保证，我无法判断。）那么您可以：

1：在kernel里面加入判断，if （xIndex==0）循环A ，else 循环B。对xIndex==0的情况特别处理。

2：或者，既然我们知道最终的结果仅仅是dev_V[0]变为正常值的2倍，其他无妨，而且这是个确定的行为，那么也可以不修改device端的任何代码，只最后在host端收集最终结果的时候，人为将第一个值除以2修正即可。

祝您编码愉快！

system · 2012 年11 月 25 日 08:26

谢谢，运行成功了，第2种方法我也想到过，不过我觉得第一种方法比较好~

system · 2012 年11 月 25 日 11:40

恭喜恭喜，祝您编码愉快~