ex12.c

采用全选主元的高斯—若当方法求矩阵的逆矩阵

//  file <RINV.c>全选主元高斯-若尔当消去法实矩阵求逆

  #include "stdlib.h"
  #include "math.h"
  #include "stdio.h"

  int rinv(n,a)
	  //返回一个整型标志。若返回0，表示矩阵奇异；否则表示正常返回；
  int n;		//矩阵阶数
  double 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;}	//记忆行列交换信息
          }
        if (d+1.0==1.0)		//矩阵奇异错误返回
          { free(is); free(js); printf("err**not inv\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];
        for (j=0; j<=n-1; j++)		//归一化
          if (j!=k)
            { u=k*n+j; a[u]=a[u]*a[l];}
        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];
                }
        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);		//正常返回
  }

 void main()
  { int i,j;
	
    static double a[4][4]=
           { {1,3,2,13},
             {7,2,1,-2},
             {9,15,3,-2},
             {-2,-2,11,5} };
    static double b[4][4];	//保存原矩阵

    for (i=0; i<=3; i++)
		for (j=0; j<=3; j++)
			b[i][j]=a[i][j];
    
	i=rinv(4,a);
    
	if (i!=0)
		{ printf("MAT A IS:\n");
		  for (i=0; i<=3; i++)
			{for (j=0; j<=3; j++)
				 printf("%10.1e ",b[i][j]);
				 printf("\n");
			}
        printf("\n");
        printf("MAT A- IS:\n");
        for (i=0; i<=3; i++)
			{for (j=0; j<=3; j++)
             printf("%15.5e ",a[i][j]);
             printf("\n");
			}
        printf("\n");
		}
  }