矩阵相乘出现的问题

system · 2013 年6 月 30 日 03:03

用一个8048的矩阵和48120的矩阵相乘，为了验证算法先把矩阵值都设置为1，那么结果肯定都应该是48，但是在具体实践中遇到问题，代码如下：

__global__ void MatrixMul_Kernel(float *dst_c,const float* src_a,const float* src_b,int width_a,int width_c)
{ //用的是2D的,不然不好分块
	int BLOCK_SIZE=16;
	int bidx=blockIdx.x;
	int bidy=blockIdx.y;
	int tidx=threadIdx.x;
	int tidy=threadIdx.y;
	int i,j;
	int n=(width_a+BLOCK_SIZE-1)/BLOCK_SIZE; //循环次数
	__shared__ float As[16][16];
	__shared__ float Bs[16][16];

	for (i=0;i<n;i++){
		As[tidy][tidx]=1/*src_a[(bidy*blockDim.y+tidy)*width_a+i*BLOCK_SIZE+tidx]*/; //y一样,x(0,width)
		Bs[tidy][tidx]=1/*src_b[(i*BLOCK_SIZE+tidy)*width_c+bidx*blockDim.x+tidx]*/; //x一样,y(0,width) 
		__syncthreads();
		for (j=0;j<BLOCK_SIZE;j++){
			dst_c[(bidy*blockDim.y+tidy)*width_c+bidx*blockDim.x+tidx]+=As[tidy][j]*Bs[j][tidx];
		}
		__syncthreads();
	}
}

在这里面如果把AS、BS的值写成：

As[tidy][tidx]=1//y一样,x(0,width)
Bs[tidy][tidx]=1; //x一样,y(0,width)

结果就是48，但是如果这样写：

As[tidy][tidx]=src_a[(bidy*blockDim.y+tidy)*width_a+i*BLOCK_SIZE+tidx]; //y一样,x(0,width)
Bs[tidy][tidx]=src_b[(i*BLOCK_SIZE+tidy)*width_c+bidx*blockDim.x+tidx]; //x一样,y(0,width)

结果却是1.#QNAN。。但是对于这种情况用Nsight进行调试值确实是1，截图如下：
[attach]3222[/attach]
这是怎么回事呢？谢谢

system · 2013 年6 月 30 日 05:08

楼主您好，因为您没有提供您的src_a, src_b的大小，以及其他一些相关数据，我仅能猜测您的问题可能出现在如下方面：

（1）src_a和src_b中的数据没有初始化，此时您得到NAN是正常的。
（2）假设条件(1)不成立。那么您的一些数据访问可能超出有效范围。
举个例子说，您的width_a 无法被 BLOCK_SIZE 整除，而您使用了向上取整：
int n=(width_a+BLOCK_SIZE-1)/BLOCK_SIZE;
（例如width_a = 1000, BLOCK_SIZE = 256, 则n = 4)
那么在下面的循环中。src_a[ … * width_a + …]就有可能超过总有效元素的范围，即您越界了，
而后面的数据不一定无法访问，可能是有些的地址，但里面的是垃圾数据，从而导致您出现NAN。

建议您检查此2方面，
感谢来访，周末愉快。

system · 2013 年6 月 30 日 06:48

谢谢版主的解答，果然是是src_b不对，在这里：


	checkCudaErrors(cudaMemcpy(Matrix_db,Matrix_hb,bw*bh*sizeof(float),cudaMemcpyHostToDevice));

写成

checkCudaErrors(cudaMemcpy(Matrix_da,Matrix_ha,bw*bh*sizeof(float),cudaMemcpyHostToDevice))

囧。。实际在代码中根本就没有对Matrix_db赋值。。谢谢版主。。

system · 2013 年6 月 30 日 06:49

感谢来访，周末愉快。

（以及上述解答不表示您完全不存在其他问题，而是可能存在哪些问题）

system · 2013 年6 月 30 日 07:07

囧。。果然还有问题。。哎。。版主厉害啊。。后面用了一个800480和4801200的矩阵相乘，同样把矩阵值设置成1，这次结果倒是480，但是这里在CPU上统计的时间是1141.11ms，而核函数统计的时间是0.0245755ms，两个时间不能相差这么多啊。代码如下：
这个是主函数：

void CDICView::OnMatrixmulGpu()
{
	// TODO: Add your command handler code here
	//////////////////////////////////////////////////////////////////////////
	//给CPU分配空间并赋值
	int i,j;
	int BLOCK_SIZE=16;
	int aw=30*BLOCK_SIZE;
	int ah=50*BLOCK_SIZE;
	int bw=80*BLOCK_SIZE;
	int bh=aw;
	int cw=bw;
	int ch=ah;
	float *Matrix_ha=(float*)malloc(aw*ah*sizeof(float));  //尽量用c语言,这样通用性比较强
	float *Matrix_hb=(float*)malloc(bw*bh*sizeof(float));
	float *Matrix_hc=(float*)malloc(cw*ch*sizeof(float));
	for (i=0;i<ah;i++){ //a
		for (j=0;j<aw;j++){
			Matrix_ha[i*aw+j]=1/*(float)rand()/RAND_MAX*/;
		}
	}
	for (i=0;i<bh;i++){ //b
		for (j=0;j<bw;j++){
			Matrix_hb[i*bw+j]=1/*(float)rand()/RAND_MAX*/;
		}
	}
	memset(Matrix_hc,0,cw*ch*sizeof(float)); //c

	ofstream cpu_time("cputim.txt");

	//计算CPU上的数据和计算时间
	_LARGE_INTEGER  cpu_start_time;
	_LARGE_INTEGER cpu_end_time;
	LARGE_INTEGER f;
	QueryPerformanceFrequency(&f);
	double dqFreq=(double)f.QuadPart;

	QueryPerformanceCounter(&cpu_start_time);
	double s1;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			s1=0;
			for (int k=0;k<aw;k++){
				s1=s1+Matrix_ha[i*aw+k]*Matrix_hb[k*bw+j];
			}
			Matrix_hc[i*cw+j]=s1;
		}
	}
	QueryPerformanceCounter(&cpu_end_time);
	cpu_time<<"计算时间为:"<<(cpu_end_time.QuadPart-cpu_start_time.QuadPart)/dqFreq*1000<<"ms"<<endl;

	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			cpu_time<<Matrix_hc[i*cw+j]<<"  ";
		}
		cpu_time<<endl;
	}
	AfxMessageBox("CPU计算完成!");


	//////////////////////////////////////////////////////////////////////////
	//GPU
	MatrixMul(Matrix_hc,Matrix_ha,Matrix_hb);

	//交出空间
	free(Matrix_ha);
	free(Matrix_hb);
	free(Matrix_hc);

	AfxMessageBox("GPU计算完成!");

}

下面是调用函数的：

void MatrixMul(float *Matrix_hc,const float *Matrix_ha,const float *Matrix_hb)
{
	int i,j;
	int BLOCK_SIZE=16;
	int aw=30*BLOCK_SIZE;
	int ah=50*BLOCK_SIZE;
	int bw=80*BLOCK_SIZE;
	int bh=aw;
	int cw=bw;
	int ch=ah;
	float *Matrix_da,*Matrix_db,*Matrix_dc;
	ofstream MatrixSub("MatrixSub.txt");


	checkCudaErrors(cudaMalloc((void**)&Matrix_da,aw*ah*sizeof(float)));
	checkCudaErrors(cudaMalloc((void**)&Matrix_db,bw*bh*sizeof(float)));
	checkCudaErrors(cudaMalloc((void**)&Matrix_dc,cw*ch*sizeof(float)));
	checkCudaErrors(cudaMemcpy(Matrix_da,Matrix_ha,aw*ah*sizeof(float),cudaMemcpyHostToDevice));
	checkCudaErrors(cudaMemcpy(Matrix_db,Matrix_hb,bw*bh*sizeof(float),cudaMemcpyHostToDevice));
	checkCudaErrors(cudaMemset(Matrix_dc,0,cw*ch*sizeof(float)));

	dim3 BKSize((cw+BLOCK_SIZE-1)/BLOCK_SIZE,(ch+BLOCK_SIZE-1)/BLOCK_SIZE);
	dim3 THSize(BLOCK_SIZE,BLOCK_SIZE);
	_LARGE_INTEGER gpu_start_time,gpu_end_time;
	LARGE_INTEGER f;
	QueryPerformanceFrequency(&f);
	double dqFreq=(double)f.QuadPart;
	QueryPerformanceCounter(&gpu_start_time);
	MatrixMul_Kernel<<<BKSize,THSize>>>(Matrix_dc,Matrix_da,Matrix_db,aw,cw);
	QueryPerformanceCounter(&gpu_end_time);
	MatrixSub<<"计算时间为:"<<(gpu_end_time.QuadPart-gpu_start_time.QuadPart)/dqFreq*1000<<"ms"<<endl;

	//传值回来
	float* Matrix_hc1=(float*)malloc(cw*ch*sizeof(float));
	checkCudaErrors(cudaMemcpy(Matrix_hc1,Matrix_dc,cw*ch*sizeof(float),cudaMemcpyDeviceToHost));

	//验证值
	MatrixSub<<"CPU计算值为:"<<endl;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			MatrixSub<<(Matrix_hc[i*cw+j])<<"   ";
		}
		MatrixSub<<endl;
	}
	MatrixSub<<"GPU计算值为:"<<endl;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			MatrixSub<<(Matrix_hc1[i*cw+j])<<"   ";
		}
		MatrixSub<<endl;
	}
	MatrixSub<<"差值为:"<<endl;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			MatrixSub<<(Matrix_hc1[i*cw+j]-Matrix_hc[i*cw+j])<<"   ";
		}
		MatrixSub<<endl;
	}


	//释放空间
	checkCudaErrors(cudaFree(Matrix_da));
	checkCudaErrors(cudaFree(Matrix_db));
	checkCudaErrors(cudaFree(Matrix_dc));
	free(Matrix_hc1);
}

下面这个是核函数：

__global__ void MatrixMul_Kernel(float *dst_c,const float* src_a,const float* src_b,int width_a,int width_c)
{ //用的是2D的,不然不好分块
	int BLOCK_SIZE=16;
	int bidx=blockIdx.x;
	int bidy=blockIdx.y;
	int tidx=threadIdx.x;
	int tidy=threadIdx.y;
	int i,j;
	int n=(width_a+BLOCK_SIZE-1)/BLOCK_SIZE; //循环次数
	__shared__ float As[16][16];
	__shared__ float Bs[16][16];

	for (i=0;i<n;i++){
		As[tidy][tidx]=src_a[(bidy*blockDim.y+tidy)*width_a+i*BLOCK_SIZE+tidx]; //y一样,x(0,width)
		Bs[tidy][tidx]=src_b[(i*BLOCK_SIZE+tidy)*width_c+bidx*blockDim.x+tidx]; //x一样,y(0,width) 
		__syncthreads();
		for (j=0;j<BLOCK_SIZE;j++){
			dst_c[(bidy*blockDim.y+tidy)*width_c+bidx*blockDim.x+tidx]+=As[tidy][j]*Bs[j][tidx];
		}
		__syncthreads();
	}
}

为什么两个统计的时间相差这么多呢？怎么感觉都不是很正常。。谢谢

system · 2013 年6 月 30 日 07:08

囧。。果然还有问题。。哎。。版主厉害啊。。后面用了一个800480和4801200的矩阵相乘，同样把矩阵值设置成1，这次结果倒是480，但是这里在CPU上统计的时间是1141.11ms，而核函数统计的时间是0.0245755ms，两个时间不能相差这么多啊。代码如下：
这个是主函数：

void CDICView::OnMatrixmulGpu()
{
	// TODO: Add your command handler code here
	//////////////////////////////////////////////////////////////////////////
	//给CPU分配空间并赋值
	int i,j;
	int BLOCK_SIZE=16;
	int aw=30*BLOCK_SIZE;
	int ah=50*BLOCK_SIZE;
	int bw=80*BLOCK_SIZE;
	int bh=aw;
	int cw=bw;
	int ch=ah;
	float *Matrix_ha=(float*)malloc(aw*ah*sizeof(float));  //尽量用c语言,这样通用性比较强
	float *Matrix_hb=(float*)malloc(bw*bh*sizeof(float));
	float *Matrix_hc=(float*)malloc(cw*ch*sizeof(float));
	for (i=0;i<ah;i++){ //a
		for (j=0;j<aw;j++){
			Matrix_ha[i*aw+j]=1/*(float)rand()/RAND_MAX*/;
		}
	}
	for (i=0;i<bh;i++){ //b
		for (j=0;j<bw;j++){
			Matrix_hb[i*bw+j]=1/*(float)rand()/RAND_MAX*/;
		}
	}
	memset(Matrix_hc,0,cw*ch*sizeof(float)); //c

	ofstream cpu_time("cputim.txt");

	//计算CPU上的数据和计算时间
	_LARGE_INTEGER  cpu_start_time;
	_LARGE_INTEGER cpu_end_time;
	LARGE_INTEGER f;
	QueryPerformanceFrequency(&f);
	double dqFreq=(double)f.QuadPart;

	QueryPerformanceCounter(&cpu_start_time);
	double s1;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			s1=0;
			for (int k=0;k<aw;k++){
				s1=s1+Matrix_ha[i*aw+k]*Matrix_hb[k*bw+j];
			}
			Matrix_hc[i*cw+j]=s1;
		}
	}
	QueryPerformanceCounter(&cpu_end_time);
	cpu_time<<"计算时间为:"<<(cpu_end_time.QuadPart-cpu_start_time.QuadPart)/dqFreq*1000<<"ms"<<endl;

	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			cpu_time<<Matrix_hc[i*cw+j]<<"  ";
		}
		cpu_time<<endl;
	}
	AfxMessageBox("CPU计算完成!");


	//////////////////////////////////////////////////////////////////////////
	//GPU
	MatrixMul(Matrix_hc,Matrix_ha,Matrix_hb);

	//交出空间
	free(Matrix_ha);
	free(Matrix_hb);
	free(Matrix_hc);

	AfxMessageBox("GPU计算完成!");

}

下面是调用函数的：

void MatrixMul(float *Matrix_hc,const float *Matrix_ha,const float *Matrix_hb)
{
	int i,j;
	int BLOCK_SIZE=16;
	int aw=30*BLOCK_SIZE;
	int ah=50*BLOCK_SIZE;
	int bw=80*BLOCK_SIZE;
	int bh=aw;
	int cw=bw;
	int ch=ah;
	float *Matrix_da,*Matrix_db,*Matrix_dc;
	ofstream MatrixSub("MatrixSub.txt");


	checkCudaErrors(cudaMalloc((void**)&Matrix_da,aw*ah*sizeof(float)));
	checkCudaErrors(cudaMalloc((void**)&Matrix_db,bw*bh*sizeof(float)));
	checkCudaErrors(cudaMalloc((void**)&Matrix_dc,cw*ch*sizeof(float)));
	checkCudaErrors(cudaMemcpy(Matrix_da,Matrix_ha,aw*ah*sizeof(float),cudaMemcpyHostToDevice));
	checkCudaErrors(cudaMemcpy(Matrix_db,Matrix_hb,bw*bh*sizeof(float),cudaMemcpyHostToDevice));
	checkCudaErrors(cudaMemset(Matrix_dc,0,cw*ch*sizeof(float)));

	dim3 BKSize((cw+BLOCK_SIZE-1)/BLOCK_SIZE,(ch+BLOCK_SIZE-1)/BLOCK_SIZE);
	dim3 THSize(BLOCK_SIZE,BLOCK_SIZE);
	_LARGE_INTEGER gpu_start_time,gpu_end_time;
	LARGE_INTEGER f;
	QueryPerformanceFrequency(&f);
	double dqFreq=(double)f.QuadPart;
	QueryPerformanceCounter(&gpu_start_time);
	MatrixMul_Kernel<<<BKSize,THSize>>>(Matrix_dc,Matrix_da,Matrix_db,aw,cw);
	QueryPerformanceCounter(&gpu_end_time);
	MatrixSub<<"计算时间为:"<<(gpu_end_time.QuadPart-gpu_start_time.QuadPart)/dqFreq*1000<<"ms"<<endl;

	//传值回来
	float* Matrix_hc1=(float*)malloc(cw*ch*sizeof(float));
	checkCudaErrors(cudaMemcpy(Matrix_hc1,Matrix_dc,cw*ch*sizeof(float),cudaMemcpyDeviceToHost));

	//验证值
	MatrixSub<<"CPU计算值为:"<<endl;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			MatrixSub<<(Matrix_hc[i*cw+j])<<"   ";
		}
		MatrixSub<<endl;
	}
	MatrixSub<<"GPU计算值为:"<<endl;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			MatrixSub<<(Matrix_hc1[i*cw+j])<<"   ";
		}
		MatrixSub<<endl;
	}
	MatrixSub<<"差值为:"<<endl;
	for (i=0;i<ch;i++){
		for (j=0;j<cw;j++){
			MatrixSub<<(Matrix_hc1[i*cw+j]-Matrix_hc[i*cw+j])<<"   ";
		}
		MatrixSub<<endl;
	}


	//释放空间
	checkCudaErrors(cudaFree(Matrix_da));
	checkCudaErrors(cudaFree(Matrix_db));
	checkCudaErrors(cudaFree(Matrix_dc));
	free(Matrix_hc1);
}

下面这个是核函数：

__global__ void MatrixMul_Kernel(float *dst_c,const float* src_a,const float* src_b,int width_a,int width_c)
{ //用的是2D的,不然不好分块
	int BLOCK_SIZE=16;
	int bidx=blockIdx.x;
	int bidy=blockIdx.y;
	int tidx=threadIdx.x;
	int tidy=threadIdx.y;
	int i,j;
	int n=(width_a+BLOCK_SIZE-1)/BLOCK_SIZE; //循环次数
	__shared__ float As[16][16];
	__shared__ float Bs[16][16];

	for (i=0;i<n;i++){
		As[tidy][tidx]=src_a[(bidy*blockDim.y+tidy)*width_a+i*BLOCK_SIZE+tidx]; //y一样,x(0,width)
		Bs[tidy][tidx]=src_b[(i*BLOCK_SIZE+tidy)*width_c+bidx*blockDim.x+tidx]; //x一样,y(0,width) 
		__syncthreads();
		for (j=0;j<BLOCK_SIZE;j++){
			dst_c[(bidy*blockDim.y+tidy)*width_c+bidx*blockDim.x+tidx]+=As[tidy][j]*Bs[j][tidx];
		}
		__syncthreads();
	}
}

为什么两个统计的时间相差这么多呢？怎么感觉都不是很正常。。谢谢

system · 2013 年6 月 30 日 07:14

楼主您好，您的计时方式有问题，（这也是我为何不建议各位楼主自己测时的原因）。

QueryPerformanceCounter(&gpu_start_time);
MatrixMul_Kernel<<<>>>(…);
QueryPerformanceCounter(&gpu_end_time);

您的2个QueryPerformanceCounter之间，只是kernel开始启动执行了，而不是他完全执行完毕了。换句话说，<<<>>>不等执行完毕就会返回的，您需要某种形式的同步，加入一句：
cudaError_t result = cudaDeviceSynchronize();
是个不错的行为，

建议尝试下。

system · 2013 年6 月 30 日 07:26

嗯。。感谢版主在周末及时准确的回答。。。后面像版主讲的加入了同步操作，结果比较正常了，GPU的测试时间变成了95.9632ms，也体会到了原来讲的异步操作。版主周末快乐。。

system · 2013 年6 月 30 日 07:28

您客气了。服务您是我们的荣幸。

感谢您的来访，也祝您周末愉快。