matrixinverse.txt

矩阵求逆 采用追赶法实现矩阵求逆的快速算法

#include<stdio.h>
#define rows 4
#define columns 4
#define E 0.000001
float input_Real[rows][columns]={{0.1934f,0.1509f,0.8537f,0.8216f},{0.6822f,0.6979f,0.5936f,0.6449f},{0.3028f,0.3784f,0.4966f,0.8180f},{0.5417f,0.8600f,0.8998f,0.6602f}};
float input_Imag[rows][columns]={{0.3420f,0.7271f,0.3704f,0.6946f},{0.2897f,0.3093f,0.7027f,0.6213f},{0.3412f,0.8385f,0.5466f,0.7948f},{0.5341f,0.5681f,0.4449f,0.9568f}};
float output_Real[rows][columns];
float output_Imag[rows][columns];
//**Print Complex Matrix**//
void printComplexMatrix(float re[][columns] ,float im[][columns])
{                                                              
    int i=0,j=0;
    for(i=0;i<rows;i++) 
	{
		for(j=0;j<columns;j++)
			printf("%.4f+i(%.4f) ",re[i][j],im[i][j]);
		printf("\n");
	}
}

//*to get absolute value
float abs(float Number)
{
	if(Number>=0.0f)
		return Number;
	else return -1*Number;
}
//* to get inv(L),inv(U)*//
void inverseLU(float L[][2*columns],float U[][2*columns],float invL[][2*columns],float invU[][2*columns])
{
	for(int n=0;n<2*rows;n++)
	{
		invL[n][n]=1;invU[n][n]=1;
	}
	for(n=0;n<2*rows;n++)
		for (int i=n+1;i<2*rows;i++)
			for(int j=0;j<i;j++)
				invL[i][j]=invL[i][j]-L[i][n]*invL[n][j];

	for(int i=0;i<2*rows;i++)
		invU[i][i]=invU[i][i]/U[i][i];
	for(n=2*rows-1;n>-1;n--)
		for(i=0;i<n;i++)
			for(int j=n;j<2*rows;j++)
				invU[i][j]=invU[i][j]-invU[n][j]*U[i][n]/U[i][i];
}

//*Get the inverse of a matrix *//
int inverseComplexMatrix(float re[][columns] ,float im[][columns],float invRe[][columns],float invIm[][columns])
{
	int M=2*rows,N=2*columns,i=0,j=0,m=0,k=0,Imax;
	float temp;
	float A[2*rows][2*columns]={0.0f},L[2*rows][2*columns]={0.0f},
		invL[2*rows][2*columns]={0.0f},invU[2*rows][2*columns]={0.0f},
		invPA[2*rows][2*columns]={0.0f},invA[2*rows][2*columns]={0.0f};
	int IP[2*rows][2];for(i=0;i<M;i++){IP[i][0]=i;IP[i][1]=i;}
	for(i=0;i<rows;i++)
		for(j=0;j<columns;j++)
		{
			A[i][j]=re[i][j];
			A[i][j+columns]=-1*im[i][j];
			A[i+rows][j]=im[i][j];
			A[i+rows][j+columns]=re[i][j];
		}
	//*A is decomposed by using LU*//
	for(m=0;m<M;m++)
	{
		L[m][m]=1.0f;
		temp=A[m][m];
		for(i=m+1;i<M;i++)
		{
			if(abs(temp)<abs(A[i][m]))
			{
				k=i;
				temp=A[i][m];
			}
		}
		if(abs(temp)<E)
			return 0;
		//*to get IP(Swap Matrix)*
		Imax=IP[m][1];
		IP[m][1]=IP[k][1];
		IP[k][1]=Imax;
		//* end of geting IP
		//* to slect the main-element
		for(j=m;j<N;j++)
		{
			temp=A[m][j];	A[m][j]=A[k][j];	A[k][j]=temp;
		}
		for(i=m+1;i<M;i++)
		{
			L[i][m]=A[i][m]/A[m][m];
			for(int j=m;j<2*rows;j++)
				A[i][j]=A[i][j]-A[m][j]*L[i][m];
		}
	}
	//* LU decomposing is end//
	//* to get inverse of matrix L and U//
	inverseLU(L,A,invL,invU);
	//* invPA=inv(U)*inv(L)
	for(i=0;i<M;i++)
		for(j=0;j<N;j++)
		{
			m=i;
			if(i<j)m=j;
			while(m<M)
			{
				invPA[i][j]=invPA[i][j]+invU[i][m]*invL[m][j];
				m++;
			}
		}
	//* invA=invPA*P
	for(i=0;i<M;i++)
	{
		if(IP[i][0]!=IP[i][1])
			for(j=0;j<M;j++)
			{temp=invPA[j][i];invPA[j][i]=invPA[j][IP[i][1]];invPA[j][IP[i][1]]=temp;}
	}
	for(i=0;i<rows;i++)
		for(j=0;j<columns;j++)
		{
			invRe[i][j]=invPA[i][j];
			invIm[i][j]=invPA[i+columns][j];
		}
	return 1;	
}
int main()
{
	printf("The original Matrix A:\n");
	printComplexMatrix(input_Real,input_Imag);
	if(inverseComplexMatrix(input_Real,input_Imag,output_Real,output_Imag))
	{	
		printf("The inverse of Matrix A:\n");
		printComplexMatrix(output_Real,output_Imag);
	}
	else printf("The input matrix is badly mobid");
	return 0;
}