提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档
@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);