inverse.cpp

基于高斯消元法的可以进行复数矩阵求逆子程序

void inverse(int c)//复数矩阵求逆的子函数
{
	int i,j,k,m,row;
	row=3*c;
	float t;
	float **matrix=NULL;//增广矩阵
	float **cr_real=NULL,**cr_image=NULL;
	extern float **re,**im;

	matrix=new float*[row];
	for(i=0;i<row;i++)
		matrix[i]=new float[row];
	

	for(i=0;i<3*c;i++)
		for(j=0;j<3*c;j++)
	{
			matrix[i][j]=0;
	}

	
	
    for(i=0;i<c;i++)  
	{
		for(j=0;j<c;j++)
			matrix[i][j]=re[i][j];
		for(j=c;j<2*c;j++)
			matrix[i][j]=-1.0*im[i][j-c];
		matrix[i][i+2*c]=1;
	}
	for(i=c;i<2*c;i++)
	{
		for(j=0;j<c;j++)
			matrix[i][j]=im[i-c][j];
		for(j=c;j<2*c;j++)
			matrix[i][j]=re[i-c][j-c];
	}
	//对M做2c次列主元Gauss消去
	for(i=0;i<c;i++)
		matrix[i+2*c][i]=1;
	for(j=0;j<2*c;j++)
	{
        for(i=j;i<2*c;i++)
		{
			if(fabs(matrix[i][j])>fabs(matrix[j][j]))
			{	
				m=i;
			break;
			}
			else 
				m=j;  
		}                            //找到主元所在行
		for(k=0;k<3*c;k++)
		{
			t=matrix[j][k];                       //主元所在行与对角元所在行交换
            matrix[j][k]=matrix[m][k];
		    matrix[m][k]=t;
		}
	    for(i=j+1;i<2*c;i++)
			for(k=3*c-1;k>=j;k--)
				matrix[i][k]=matrix[i][k]-matrix[j][k]*matrix[i][j]/matrix[j][j];//化为上三角阵
	}
		
		
	for(j=2*c-1;j>0;j--)
		for(i=0;i<j;i++)
			for(k=3*c-1;k>=0;k--)
				if(matrix[j][j]!=0)
				matrix[i][k]=matrix[i][k]-matrix[j][k]*matrix[i][j]/matrix[j][j];//化为单位阵
		
	for(i=0;i<2*c;i++)
		for(k=3*c-1;k>=0;k--)
		matrix[i][k]=matrix[i][k]/matrix[i][i];         //对角元化为1
	for(i=0;i<c;i++)
		for(j=0;j<c;j++)
		{
			re[i][j]=matrix[i][j+2*c];
			im[i][j]=matrix[i+c][j+2*c];
		}
		for(i=0;i<row;i++)
		    delete []matrix[i];
		   delete []matrix;
		   }