请教一个复数矩阵求逆的问题。。。

最近在写一个复数矩阵求逆的kernel函数,利用高斯消元法做的,(A|I)---->(I|INVA),能变换出单位矩阵,但是求出的逆矩阵不对,不知道是不是受精度的影响。原来的矩阵是共轭对称的,变换出的逆矩阵也应该是共轭对称,结果不对,有路过的高手指点下迷津吧,多谢多谢!
代码如下:
// 矩阵求逆之归一操作
global void MatrixInverse_Normalized(Complex* A,Complex* B,int n,int i)
{
int tix = blockDim.x * blockIdx.x + threadIdx.x;
int tiy = blockDim.y * blockIdx.y + threadIdx.y;
Complex temp;//归一化值
float abs_temp;
float p,q,s,m,l,k;
temp = A[i*n+i];
abs_temp = temp.x * temp.x + temp.y * temp.y;
//对A,B矩阵第i行做归一化操作
if(tiy<n && tix<n)
{
if (tix == i)
{
if(tiy == i)
{
A[tiyn+tix].x = 1.0f;
A[tiy
n+tix].y = 0.0f;
}
else
{
p = A[tiyn+tix].xtemp.x/abs_temp;
q = -A[tiyn+tix].ytemp.y/abs_temp;
s = (A[tiyn+tix].x+A[tiyn+tix].y)(temp.x-temp.y)/abs_temp;
A[tiy
n+tix].x = p-q;
A[tiyn+tix].y = s-p-q;
}
m = B[tiy
n+tix].xtemp.x/abs_temp;
l = -B[tiy
n+tix].ytemp.y/abs_temp;
k = (B[tiy
n+tix].x+B[tiyn+tix].y)(temp.x-temp.y)/abs_temp;
B[tiyn+tix].x = m-l;
B[tiy
n+tix].y = k-m-l;
}
}
}
// 矩阵求逆之相消操作
global void MatrixInverse_Elimination(Complex* A,Complex* B,int n,int i)
{
//共享存储矩阵块
shared Complex As[BLOCK_SIZE_1][BLOCK_SIZE_1];
shared Complex Bs[BLOCK_SIZE_1][BLOCK_SIZE_1];
int tix = blockDim.x * blockIdx.x + threadIdx.x;
int tiy = blockDim.y * blockIdx.y + threadIdx.y;

float a,b,c,d;
if(tiy<n && tix<n)
{
//拷贝block到共享存储区
As[threadIdx.y][threadIdx.x] = A[tiyn+tix];
Bs[threadIdx.y][threadIdx.x] = B[tiy
n+tix];
__syncthreads();
a = A[i*n+tix].x * A[tiyn+i].x - A[i*n+tix].y * A[tiyn+i].y;
b = A[i*n+tix].x * A[tiyn+i].y + A[i*n+tix].y * A[tiyn+i].x;
c = A[i*n+tix].x * B[tiyn+i].x - A[i*n+tix].y * B[tiyn+i].y;
d = A[i*n+tix].x * B[tiyn+i].y + A[i*n+tix].y * B[tiyn+i].x;

//用归一化的第i行与非i行相消
As[threadIdx.y][threadIdx.x].x = As[threadIdx.y][threadIdx.x].x - a;
As[threadIdx.y][threadIdx.x].y = As[threadIdx.y][threadIdx.x].y - b;
Bs[threadIdx.y][threadIdx.x].x = Bs[threadIdx.y][threadIdx.x].x - c;
Bs[threadIdx.y][threadIdx.x].y = Bs[threadIdx.y][threadIdx.x].y - d;
__syncthreads();
if(tix != i) //将相消结果存入显存
{
A[tiyn+tix] = As[threadIdx.y][threadIdx.x];
B[tiy
n+tix] = Bs[threadIdx.y][threadIdx.x];
}
}
}
主函数中调用
for(int i=0;i<50;i++)
{
MatrixInverse_Normalized(A,B,i);
MatrixInverse_Elimination(A,B,i);
}

你仔细检查一下你相消的地方,对不准就很容易消错的。

相消都验证了,能化为单位矩阵,但是逆矩阵和matlab上的结果对不是,怀疑是不是精度的问题,matlab是用double型来算的。。。。

有一定的误差是很正常的,要看你误差是否可接受,呵呵

建议先将单线程版本(kernel只用一个线程)调对,再做并行版本。

谢谢,这倒是个好办法。。。现在是直接用CPU计算的,再把值传回来,感觉时间也不慢~呵呵