c언어 역행렬 구하는 함수

#include <stdlib.h>
#include <stdio.h>

int InvMatrix(int n, const double* A, double* b)  // 역행렬 구하는 함수
{
    double m;
    register int i, j, k;
    double* a = new double[n*n];

    if(a==NULL) 
  return 0;
    for(i=0; i<n*n; i++) 
  a[i]=A[i];
 for(i=0; i<n; i++) 
 {
  for(j=0; j<n; j++)
  {
            b[i*n+j]=(i==j)?1.:0.;
        }
    }
    for(i=0; i<n; i++)
 {
  if(a[i*n+i]==0.) 
  {
   if(i==n-1) 
   {
    delete[] a;
    return 0;
            }
            for(k=1; i+k<n; k++)
   {
                if(a[i*n+i+k] != 0.) 
     break;
            }
            if(i+k>=n)
   {
                delete[] a;
    return 0;
            }
            for(j=0; j<n; j++) 
   {
                m = a[i*n+j];
                a[i*n+j] = a[(i+k)*n+j];
                a[(i+k)*n+j] = m;
                m = b[i*n+j];
                b[i*n+j] = b[(i+k)*n+j];
                b[(i+k)*n+j] = m;
            }
        }
        m = a[i*n+i];
        for(j=0; j<n; j++) 
  {
            a[i*n+j]/=m;
            b[i*n+j]/=m;
        }
        for(j=0; j<n; j++) 
  {
            if(i==j) 
    continue;

            m = a[j*n+i];
            for(k=0; k<n; k++)   
   {
                a[j*n+k] -= a[i*n+k]*m;
                b[j*n+k] -= b[i*n+k]*m;
            }
        }
    }
    delete[] a;
    return 1;
}

 

void main()
{
 double a[4][4], b[4][4];

 // .....

 //InvMatrix(4, (double*)a, (double*)b); // a 배열의 역행렬 구해서 b에 저장
}

 

 

 

 

 

--------------------------------------------------------------------------------

 
 
호출 예)
 
double a[4][4], b[4][4];
 
// .....
 
InvMatrix(4, (double*)a, (double*)b);  // a 배열의 역행렬 구해서 b에 저장