inverse.c

推荐刚刚开始学C语言的同志们看下

#include <stdio.h>
#include <malloc.h>
#include <math.h>

void MatrixMul(a,b,m,n,k,c)  /*实矩阵相乘*/
int m,n,k; /*m:矩阵A的行数, n:矩阵B的行数, k:矩阵B的列数*/
double a[],b[],c[]; /*a为A矩阵, b为B矩阵, c为结果，即c = AB */
{ 
	int i,j,l,u;
	/*逐行逐列计算乘积*/
	for (i=0; i<=m-1; i++)
		for (j=0; j<=k-1; j++)
		{
			u=i*k+j; c[u]=0.0;
			for (l=0; l<=n-1; l++)
				c[u]=c[u]+a[i*n+l]*b[l*k+j];
		}
		return;
}
int brinv(a,n) /*求矩阵的逆矩阵*/
int n; /*矩阵的阶数*/
double a[]; /*矩阵A*/
{ 
	int *is,*js,i,j,k,l,u,v;
    double d,p;
    is=malloc(n*sizeof(int));
    js=malloc(n*sizeof(int));
    for (k=0; k<=n-1; k++)
	{ 
		d=0.0;
		for (i=k; i<=n-1; i++)
			/*全选主元，即选取绝对值最大的元素*/
			for (j=k; j<=n-1; j++)
			{ 
				l=i*n+j; p=fabs(a[l]);
				if (p>d) { d=p; is[k]=i; js[k]=j;}
			}
		/*全部为0，此时为奇异矩阵*/
		if (d+1.0==1.0)
		{ 
			free(is); free(js); printf("这是奇异矩阵，无法求逆矩阵。\n");
			return(0);
		}
		/*行交换*/
		if (is[k]!=k)
			for (j=0; j<=n-1; j++)
			{ 
				u=k*n+j; v=is[k]*n+j;
				p=a[u]; a[u]=a[v]; a[v]=p;
			}
		/*列交换*/
		if (js[k]!=k)
			for (i=0; i<=n-1; i++)
			{
				u=i*n+k; v=i*n+js[k];
				p=a[u]; a[u]=a[v]; a[v]=p;
			}
		l=k*n+k;
		a[l]=1.0/a[l]; /*求主元的倒数*/
		/* a[kj]a[kk] -> a[kj] */
		for (j=0; j<=n-1; j++)
			if (j!=k)
			{ 
				u=k*n+j; a[u]=a[u]*a[l];
			}
		/* a[ij] - a[ik]a[kj] -> a[ij] */
		for (i=0; i<=n-1; i++)
			if (i!=k)
				for (j=0; j<=n-1; j++)
					if (j!=k)
					{ 
						u=i*n+j;
						a[u]=a[u]-a[i*n+k]*a[k*n+j];
					}
		/* -a[ik]a[kk] -> a[ik] */
		for (i=0; i<=n-1; i++)
			if (i!=k)
			{ 
				u=i*n+k; a[u]=-a[u]*a[l];
			}
	}
    for (k=n-1; k>=0; k--)
	{ 
		/*恢复列*/
		if (js[k]!=k)
			for (j=0; j<=n-1; j++)
			{ 
				u=k*n+j; v=js[k]*n+j;
				p=a[u]; a[u]=a[v]; a[v]=p;
			}
	    /*恢复行*/
		if (is[k]!=k)
			for (i=0; i<=n-1; i++)
			{
				u=i*n+k; v=i*n+is[k];
				p=a[u]; a[u]=a[v]; a[v]=p;
			}
	}
    free(is); free(js);
    return(1);
}

main()
{ 
	int i,j;
    static double a[4][4]={ {0.2368,0.2471,0.2568,1.2671},
	{1.1161,0.1254,0.1397,0.1490},
	{0.1582,1.1675,0.1768,0.1871},
	{0.1968,0.2071,1.2168,0.2271}};
    static double b[4][4],c[4][4];
    for (i=0; i<=3; i++)
		for (j=0; j<=3; j++)
			b[i][j]=a[i][j];
	i=brinv(a,4);
	if (i!=0)
	{ 
		printf("矩阵A:\n");
		for (i=0; i<=3; i++)
		{ 
			for (j=0; j<=3; j++)
				printf("%13.7f\t",b[i][j]);
			printf("\n");
		}
		printf("\n");
		printf("A的逆矩阵A-:\n");
		for (i=0; i<=3; i++)
		{
			for (j=0; j<=3; j++)
				printf("%13.7f\t",a[i][j]);
			printf("\n");
		}
		printf("\n");
		printf("矩阵A与其逆矩阵A-的乘积:\n");
		MatrixMul(b,a,4,4,4,c);
		for (i=0; i<=3; i++)
		{ 
			for (j=0; j<=3; j++)
				printf("%13.7f\t",c[i][j]);
			printf("\n");
		}
	}
}