高维矩阵求逆：Gauss-Jordan算法

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

@TOC高维矩阵求逆：Gauss-Jordan算法©

---

前言

亲测好用~~~

代码

int MatInv_GaussJordan(float *a, int n)
{
int *is,*js,i,j,k,l,u,v;
float d,p;

 is = malloc(n * sizeof(int));
 js = malloc(n * sizeof(int));
 for (k=0; k<=n-1; k++)
 {
     d=0.0f;
     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 (fabs(d) < 0.0001f) 
     { 
         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.0f/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);

}

测试用例

Mat a, b;

// MatCreate(&a, 4, 4);
// MatCreate(&b, 4, 4);
// a.element[0][0]=0.2368f; a.element[0][1]=0.2471f; a.element[0][2]=0.2568f; a.element[0][3]=1.2671f;
// a.element[1][0]=1.1161f; a.element[1][1]=0.1254f; a.element[1][2]=0.1397f; a.element[1][3]=0.1490f;
// a.element[2][0]=0.1582f; a.element[2][1]=1.1675f; a.element[2][2]=0.1768f; a.element[2][3]=0.1871f;
// a.element[3][0]=0.1968f; a.element[3][1]=0.2071f; a.element[3][2]=1.2168f; a.element[3][3]=0.2271f;

MatCreate(&a, 12, 12);
MatCreate(&b, 12, 12);
MatZeros(&a);
a.element[0][0]=181.7958f; a.element[0][1]=0.0837f;   a.element[0][2]=1.7910f;   a.element[0][9]=0.0749f;
a.element[1][0]=0.0837f;   a.element[1][1]=180.9509f; a.element[1][2]=1.8846f;   a.element[1][9]=0.5655f;
a.element[2][0]=1.7910f;   a.element[2][1]=1.8846f;   a.element[2][2]=120.8461f; a.element[2][9]=19.8508f;

a.element[3][3]=181.7958f; a.element[3][4]=0.0837f;   a.element[3][5]=1.7910f;   a.element[3][10]=0.0749f;
a.element[4][3]=0.0837f;   a.element[4][4]=180.9509f; a.element[4][5]=1.8846f;   a.element[4][10]=0.5655f;
a.element[5][3]=1.7910f;   a.element[5][4]=1.8846f;   a.element[5][5]=120.8461f; a.element[5][10]=19.8508f;

a.element[6][6]=181.7958f; a.element[6][7]=0.0837f;   a.element[6][8]=1.7910f;   a.element[6][11]=0.0749f;
a.element[7][6]=0.0837f;   a.element[7][7]=180.9509f; a.element[7][8]=1.8846f;   a.element[7][11]=0.5655f;
a.element[8][6]=1.7910f;   a.element[8][7]=1.8846f;   a.element[8][8]=120.8461f; a.element[8][11]=19.8508f;

a.element[9][0] =0.0749f;   a.element[9][1] =0.5655f;   a.element[9][2] =19.8508f;  a.element[9][9]  =5.0f;
a.element[10][3]=0.0749f;   a.element[10][4]=0.5655f;   a.element[10][5]=19.8508f;  a.element[10][10]=5.0f;
a.element[11][6]=0.0749f;   a.element[11][7]=0.5655f;   a.element[11][8]=19.8508f;  a.element[11][11]=5.0f;

int i,j;
success = MatInv_GaussJordan(&a, &b);
MatDelete(&a);

static float c[12][12];
for (i=0; i<12; i++)
{
	for (j=0; j<12; j++)
	{
		c[i][j]=b.element[i][j];
	}
}
MatDelete(&b);

---

总结