115.n阶方阵求逆

#include <stdio.h>
#include <stdlib.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[]; /*矩阵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(" >> 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[]; /*矩阵a*/
{
  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");
  }
}
main()
{
  int i,j,n=0;
  double A[MAX],B[MAX],C[MAX];
    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];
    clrscr();
    puts("**********************************************************");
    puts("*    This program is to inverse a square matrix A(nxn).  *");
    puts("**********************************************************");
    while(n<=0)
    {
      printf(" >> Please input the order n of the matrix (n>0): ");
      scanf("%d",&n);
    }
    printf(" >> Please input the elements of the matrix one by one:\n >> ");
    for(i=0;i<n*n;i++)
    {
      scanf("%lf",&A[i]);
      B[i]=A[i];
    }
    for(i=0;i<4;i++)
  for(j=0;j<4;j++)
    b[i][j]=a[i][j];
    i=brinv(A,n);
    if (i!=0)
    {
  printf("    Matrix A:\n");
  print_matrix(B,n);
  printf("\n");
  printf("    A's Inverse Matrix A-:\n");
  print_matrix(A,n);
  printf("\n");
  printf("    Product of A and A- :\n");
  MatrixMul(B,A,n,n,n,C);
  print_matrix(C,n);
    }
    printf("\n Press any key to quit...");
    getch();
}