nni.c

计算n阶矩阵的逆矩阵的函数

#include<stdio.h>
#include<math.h>
void ni(double *a,double *s,int n)
{
	int i,j;
	double b,temp,d;
	int k;

	//建立增广矩阵
	for(i=0;i<n;i++)	
	{
		for(k=0;k<2*n;k++)
			if(k<n)
				s[i*2*n+k]=a[i*n+k];
			else
			{
				if(i==k-n)
					s[i*2*n+k]=1;
				else
					s[i*2*n+k]=0;
			}
	}
	//完成建立增广矩阵
	
	//i从0变化到n
	for(i=0;i<n;i++)
	{
		b=s[i*2*n+i];			
		j=i;
		//2.1查找第i列绝对值最大的元素
		for(k=i;k<n;k++)		
		{	
			if(abs(s[k*2*n+i])>abs(b))
			{
				b=s[k*2*n+i];
				j=k;
			}
		
		
		}
		//2.1查找完成

		//2.3交换第i行和第j行
		if(b!=0)			
		{
			if(j!=i)
			{
				for(k=0;k<2*n;k++)
				{
					temp=s[j*2*n+k];
					s[j*2*n+k]=s[i*2*n+k];
					s[i*2*n+k]=temp;
				}
			}
		}
		//2.3交换完成

		//2.4 b!=1,i行元素除以b
		if(b!=1)
		{
			
			for(k=0;k<2*n;k++)
			{
				s[i*2*n+k]=s[i*2*n+k]/b;	
			}
		}
		//2.4

		//2.5 第i列k行为非0元素s[k][i],则第k行第j列元素为s[k][j]=s[k][j]-s[k][i]*s[k][i]
		for(k=0;k<2*n;k++)
		{
			if(k!=i&&s[k*2*n+i]!=0)
			{
				d=s[k*2*n+i];
				for(j=0;j<2*n;j++)
				{
					s[k*2*n+j]=s[k*2*n+j]-s[i*2*n+j]*d;
				}
			}
		}
		//2.5
	}
}//完成


//矩阵相乘
void mult(double *a,double *b,int n,int m,int x)
{
	int i,j,k;
	double *c;
	c=(double *)malloc(n*m*sizeof(double));
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		{
			c[i*m+j]=0;
			for(k=0;k<x;k++)
			{
				c[i*m+j]+=a[i*x+k]*b[k*m+j];	
			}
		}		
	}
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		{
			printf("%.6lf ",c[i*m+j]);
		}	
		printf("\n");
	}
}






main()
{
	double x[4][4]={{3,5,7,4},{21,323,443,55},{33,42,45,67},{421,3321,211,811}};
	double *a,*s,*at;
	int n=4;
	int i,j;
	a=(double *)malloc(n*n*sizeof(double));
	s=(double *)malloc(2*n*n*sizeof(double));
	at=(double *)malloc(n*n*sizeof(double));
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			a[i*n+j]=x[i][j];
		}		
	}
	ni(a,s,4);
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			at[i*n+j]=s[i*2*n+n+j];
		}
	}
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			a[i*n+j]=x[i][j];
		}		
	}


	for(i=0;i<n;i++){
		for(j=0;j<n;j++)
		{
			printf("%.6lf ",at[i*n+j]);
		}
		printf("\n");

	}
	mult(a,at,n,n,4);
}