#include <stdio.h> 
#include <stdlib.h> 
#include <malloc.h> 
#include <math.h> 
#define MAX 255 

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[]; /*矩陣C*/ 
{ 
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(" >> This is a singular matrix, can't be inversed!\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); 
} 

print_matrix(a,n)/*打印方陣a的元素*/ 
int n; /*矩陣的階數*/ 
double a[]; /*矩陣C或D*/ 
{ 
int i,j; 
for (i=0; i<n; i++) 
{ 
for (j=0; j<n; j++) 
printf("%13.7f\t",a[i*n+j]); 
printf("\n"); 
} 
} 

int main() 
{ 
int i,j,k,n=0,m; 
double A[MAX],B[MAX],C[MAX],D[MAX],E[MAX]; /* 題中沒有說明，這里做實數看待 */ 
system("cls"); 
puts("**********************************************************"); 
puts("* This program is to inverse a square matrix(nxn). *"); 
puts("**********************************************************"); 
while(n<=0) 
{ /* 設A是n*m矩陣則B必為m*n矩陣，這樣AB才會是方陣存在逆矩陣 */ 
printf(" >> Please input the number of rows in A, n= "); 
scanf("%d",&n); 
printf(" >> Please input the number of cols in A, m= "); 
scanf("%d",&m); 
} 
printf(" >> Please input the elements of the matrix A one by one:\n >> "); 
for(i=0;i<n*m;i++) 
scanf("%lf",&A[i]); 

printf(" >> Please input the elements of the matrix B one by one:\n >> "); 
for(i=0;i<m*n;i++) 
scanf("%lf",&B[i]); 

MatrixMul(A,B,n,m,n,C);/* A*B的結果存在C里面 */ 
for(i=0;i<n*n;i++) /* D[n*n]是C[n*n](也就是AB)的副本，當C[n*n]變換為本身的逆矩陣后再和D[n*n]做乘積運算 */ 
D[i]=C[i]; /* 如果乘積結果為單位矩陣，則說明計算結果正確 */ 
i=brinv(C,n); 
if (i!=0) 
{ 
printf(" Matrix AB:\n"); 
print_matrix(D,n); 
printf("\n"); 
printf(" AB's Inverse Matrix AB-:\n"); 
print_matrix(C,n); 
printf("\n"); 
printf(" Product of AB and AB- :\n"); 
MatrixMul(D,C,n,n,n,E); 
print_matrix(E,n); 
} 
system("pause"); 
return 0; 
}