求助矩阵求逆

加速计算 CUDA
1

kenerl函数
global void matrixInv( float* A,float* B,int i,int wA)
{
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;

const int index_tx=bxblockDim.x+tx;
const int index_ty=byblockDim.y+ty;
/int index=index_ty+(index_tx)wA;/
const float temp=A[i*wA+i];
if(index_ty<wA&&index_tx<wA)
{
B[i*wA+index_tx]=B[i*wA+index_tx]/temp;
A[i*wA+index_tx]=A[i*wA+index_tx]/temp;}
float Avalue=0.0;
float Bvalue=0.0;
// C[index_tywA+index_tx]=A[index_tywA+i]A[i*n+index_tx];
// D[index_tywA+index_tx]=A[index_tywA+i]*B[i*n+index_tx];

shared float As[BLOCK_DIM][BLOCK_DIM];
shared float Bs[BLOCK_DIM][BLOCK_DIM];
shared float Cs[BLOCK_DIM][BLOCK_DIM];
shared float Ds[BLOCK_DIM][BLOCK_DIM];
shared float Ms[BLOCK_DIM][BLOCK_DIM];
if(index_ty<wA&&index_tx<wA)
for(int m=0;m<wA/TILE_WIDTH;++m)
{
Bs[ty][tx] = B[index_tywA+(mTILE_WIDTH+tx)];
As[ty][tx] = A[index_tywA+(mTILE_WIDTH+tx)];
Cs[ty][tx] = A[(m*TILE_WIDTH+ty)wA+i];
Ds[ty][tx] = B[i*wA+(m*TILE_WIDTH+tx)];
Ms[ty][tx] = A[i*wA+(m*TILE_WIDTH+tx)];
__syncthreads();
/ int k=i-(i/TILE_WIDTH)TILE_WIDTH;/

// const int k=i%wA;
// const int j=i/wA;

Avalue=A[index_ty*wA+i]A[i*wA+index_tx];
Bvalue=A[index_tywA+i]*B[i*wA+index_tx];

 Bs[ty][tx]=Bs[ty][tx]-Bvalue;

As[ty][tx]=As[ty][tx]-Avalue;
__syncthreads();
if(index_ty!=i)
{ B[ index_tywA+index_tx]=Bs[ty][tx];
A[ index_tywA+index_tx]=As[ty][tx];
}
}

调用kernel函数
dim3 mygrid((wA+TILE_WIDTH-1) / TILE_WIDTH, (wA+TILE_WIDTH-1) / TILE_WIDTH);
dim3 myblock(TILE_WIDTH, TILE_WIDTH);
start = clock();
for(int i=0;i<wA;i++)
{
matrixInv<<<mygrid,myblock>>>(d_A,d_B,i,wA);

}

2

LZ您好：

请您叙述您的问题，单纯贴代码，谁都不知道您想干什么。

请您予以补充。

3

额，这个我疏忽了！不好意思！应该是调用的时候，这个函数有没有错误呀？我编写的是利用高斯消元法矩阵求逆，但是结果却错误，是kernel函数的问题？还是调用的问题呢？

4

LZ您好：

我不懂相关的线性代数问题及其数值求法，因此如果是算法逻辑上的问题，我无法为您建议。

如果您的代码执行本身有问题，比如kernel挂掉，请使用nsight调试，这样可以立即发现问题在哪里；或者您可以先在kernel后面执行以下cudaDeviceSynchronize()，并检查这个函数的返回值，看看kernel是否成功执行。

只能建议这些了。

祝您好运~

5

感谢你的回复！谢谢，这个问题我做过调试，在求解1616的时候没有问题，但是在3232以上的时候，就有问题了，应该是分块的问题，但是我不清楚为什么分块后就有问题了！算法的逻辑上我觉得应该没问题呀！我现在给你看看串行的吧！其实就是讲这些并行起来呀！

6

for(int i=0;i<m;i++){
float temp=P[i*m+i];
for(int j=0;j<m;j++){
I[i*m+j]=I[i*m+j]/temp;
P[i*m+j]=P[i*m+j]/temp;
}

  for (int j =0;j<m;j++){
	 if(j !=i)
 { 
	{ float a=P[j*m+i];
	 for (int k=0;k<m;k++)
	 {
	 I[j*m+k]  =I[j*m+k]-I[i*m+k]*a;
	 P[j*m+k]  =P[j*m+k]-P[i*m+k]*a;
	 }
 }
	}

}
}
这些串行代码实现了高斯矩阵求逆

7

LZ您好：

请问您分块尺寸改为32*32之后，kernel执行是否成功？

如不成功，此时结果必然不正确，请您提供报错信息，或者直接调试解决问题。

如kernel执行是成功的，但是结果不正确，则我无法为您找出算法逻辑错误，请见谅。

祝您好运~

8

分成3232的以后是可行的，但是以后的呢？像6464,128*128的呢？我不能每一个矩阵都分吧！再说了shared memery 的内存是一定的呀！所以，这应该和分的块没关系吧？谢谢~

9

LZ您好：

shared memory是容量有限的片上高速缓存，如果您需要使用shared memory，请配合硬件限制合理使用。数据量较大时，将所有的数据都缓冲进shared memory是不现实的。

本人无法对您的实际线性代数问题的算法实现给出建议，以及这个也超出了本版的支持范围。
请其他熟悉该问题的网友/斑竹/NV原厂支持人员继续与您讨论。

祝您编码顺利~

10

谢谢斑竹！~~那我应该去哪个版区呢？