Inverse Matrix

/*
1,1,0,2
3,-1,1,0
0,2,2,1
2,0,1,3
の逆行列を求める
*/

#include 
#include 

#define N 4
#define N2 2*N
double minv(double a[][N],double ainv[][N],int n);
double minv_row_trans(double b[][N2],int n,int n2,int i);
void out_res(double a[][N],int m,int n);

main()
{
 double a[][N]={{1,1,0,2},{3,-1,1,0},{0,2,2,1},{2,0,1,3}},ainv[N][N],d;
 d=minv(a,ainv,N);
 if(d==0)printf("Singular Matrix! \n");
 else{
  printf("Inverse Matrix \n");
  out_res(ainv,N,N);
  }
}

double minv(double a[][N],double ainv[][N],int n)
{
 double b[N][N2],d;
 int i,j,n2;
 n2=2*n;
 for(i=0;i<=n-1;i++){
  for(j=0;j<=n-1;j++)
    b[i][j]=a[i][j];
      for(j=n;j<=n2-1;j++)b[i][j]=0;
       b[i][i+n]=1;
  }
  for(i=0;i<=n-1;i++){
    d=minv_row_trans(b,n,n2,i);
     if(d==0)break;
  }
  for(i=0;i<=n-1;i++)
   for(j=0;j<=n-1;j++)
      ainv[i][j]=b[i][j+n];
  return d;
}

double minv_row_trans(double b[][N2],int n,int n2,int i)
{
 int i1,j,k;
 double b1,b2;
 i1=i-1;
 do{
  i1++;
  b1=b[i1][i];
  }while(b1==0 && i1＜n-1);

 if(b1!=0){
   if(i1!=i)
     for(j=0;j＜=n2-1;j++){
       b2=b[i][j];
       b[i][j]=b[i1][j];
       b[i1][j]=b2;
    }
   b1=b[i][i];
   for(j=0;j＜=n2-1;j++)
     b[i][j]=b[i][j]/b1;
   for(j=0;j＜=n-1;j++){
     if(i!=j){
        b1=b[j][i];
        for(k=0;k＜=n2-1;k++)
           b[j][k]=b[j][k]-b1*b[i][k];
       }
    }
 }
 return b1;
}

void out_res(double x[][N],int m,int n)
{
 int i,j;
 for(i=0;i＜=m-1;i++){
  for(j=0l;j＜=n-1;j++) printf("%12.5f",x[i][j]);
  printf("\n");
  }
}

[出力結果]

Inverse Matrix 
     0.45833     0.37500    -0.04167    -0.29167
     0.70833     0.12500     0.20833    -0.54167
    -0.66667     0.00000     0.33333     0.33333
    -0.08333    -0.25000    -0.08333     0.41667

Back to C Language