C/C++言語実現マトリックス求逆演算—ガウス約化/消元法
4973 ワード
サイトで検索してみると、C言語コードはオリジナルのバージョンかもしれませんが、コンパイル時に多くのエラーが報告されます.
でも大体読めるから、曲がってもいいかどうか見てみましょう.
原文のコードはマトリクス形式で保存され、実際のコードは1次元配列形式で書かれています.
直接使用できるコードは次のとおりです.
仕事でこんなに複雑な数学をやらなければならなかったので、線形代数をよく勉強しなかったのが悔しかった......
次に、変更後の結果を示します.入力パラメータはdouble**で、使用前に元のマトリクスをバックアップする必要がある場合は、マトリクスをコピーする必要があります.
でも大体読めるから、曲がってもいいかどうか見てみましょう.
原文のコードはマトリクス形式で保存され、実際のコードは1次元配列形式で書かれています.
直接使用できるコードは次のとおりです.
#include
#include
#include
// ， a
int brinv(double a[], int n)
{
int *is,*js,i,j,k,l,u,v;
double d,p;
is = new int[n];
js = new int[n];
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
");
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[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);
}
// C = A.*B
void brmul(double a[], double b[],int m,int n,int k,double c[])
{
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] += a[i*n+l]*b[l*k+j]; // ，
}
}
int main()
{
int i,j;
static double a[16]={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[16],c[16]; // a b c
for (i=0; i<=3; i++)
for (j=0; j<=3; j++)
b[4*i+j]=a[4*i+j];
cout << (i=brinv(a,4)) << endl;
if (i!=0)
{
printf("MAT A IS:
");
for (i=0; i<=3; i++)
{
for (j=0; j<=3; j++)
printf("%13.4f",b[4*i+j]);
printf("
");
}
printf("
MAT A- IS:
");
for (i=0; i<=3; i++)
{
for (j=0; j<=3; j++)
printf("%13.4f",a[4*i+j]);
printf("
");
}
printf("
MAT AA- IS:
");
brmul(b,a,4,4,4,c);
for (i=0; i<=3; i++)
{
for (j=0; j<=3; j++)
printf("%13.4f",c[4*i+j]);
printf("
");
}
}
}
仕事でこんなに複雑な数学をやらなければならなかったので、線形代数をよく勉強しなかったのが悔しかった......
次に、変更後の結果を示します.入力パラメータはdouble**で、使用前に元のマトリクスをバックアップする必要がある場合は、マトリクスをコピーする必要があります.
int matrixInversion(double **a, int n)
{
int *is = new int[n];
int *js = new int[n];
int i,j,k;
double d,p;
for ( k = 0; k < n; k++)
{
d = 0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{
p=fabs(a[i][j]);
if (p>d) { d=p; is[k]=i; js[k]=j;}
}
if ( 0.0 == d )
{
free(is); free(js); printf("err**not inv
");
return(0);
}
if (is[k]!=k)
for (j=0; j<=n-1; j++)
{
p=a[k][j];
a[k][j]=a[is[k]][j];
a[is[k]][j]=p;
}
if (js[k]!=k)
for (i=0; i<=n-1; i++)
{
p=a[i][k];
a[i][k]=a[i][js[k]];
a[i][js[k]]=p;
}
a[k][k] = 1.0/a[k][k];
for (j=0; j<=n-1; j++)
if (j!=k)
{
a[k][j] *= a[k][k];
}
for (i=0; i<=n-1; i++)
if (i!=k)
for (j=0; j<=n-1; j++)
if (j!=k)
{
a[i][j] -= a[i][k]*a[k][j];
}
for (i=0; i<=n-1; i++)
if (i!=k)
{
a[i][k] = -a[i][k]*a[k][k];
}
}
for ( k = n-1; k >= 0; k--)
{
if (js[k]!=k)
for (j=0; j<=n-1; j++)
{
p = a[k][j];
a[k][j] = a[js[k]][j];
a[js[k]][j]=p;
}
if (is[k]!=k)
for (i=0; i<=n-1; i++)
{
p = a[i][k];
a[i][k]=a[i][is[k]];
a[i][is[k]] = p;
}
}
free(is); free(js);
return(1);
}
int main()
{
//prior_testMatrix(); // ， main
int i,j;
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 = new double *[4]; // a
for (i=0; i< 4; i++)
{
b[i] = new double[4];
for (j=0; j< 4; j++)
b[i][j]=a[i][j]; // a
}
cout << (i=matrixInversion(b,4)) << endl; // ， b
if (i!=0)
{
printf("MAT A IS:
");
for (i=0; i<=3; i++)
{
for (j=0; j<=3; j++)
printf("%13.4f",a[i][j]);
printf("
");
}
printf("
MAT A- IS:
");
for (i=0; i<=3; i++)
{
for (j=0; j<=3; j++)
printf("%13.4f",b[i][j]);
printf("
");
}
}
}