最近在写一个复数矩阵求逆的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[tiyn+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[tiyn+tix].x = p-q;
A[tiyn+tix].y = s-p-q;
}
m = B[tiyn+tix].xtemp.x/abs_temp;
l = -B[tiyn+tix].ytemp.y/abs_temp;
k = (B[tiyn+tix].x+B[tiyn+tix].y)(temp.x-temp.y)/abs_temp;
B[tiyn+tix].x = m-l;
B[tiyn+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[tiyn+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[tiyn+tix] = Bs[threadIdx.y][threadIdx.x];
}
}
}
主函数中调用
for(int i=0;i<50;i++)
{
MatrixInverse_Normalized(A,B,i);
MatrixInverse_Elimination(A,B,i);
}