//////////////////////////////////////////////////////
//
// 逆行列を算出する関数
//
//////////////////////////////////////////////////////

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

//---matrix inversion & determinant---//
//逆行列の計算//
void minv(double xx[512],double eps,int row)
{
    int i, j, k, n, kaisu;
    double **a, **c;

//    printf( "元の数=" );
//    scanf( "%d", &n);

	n=row;

    /* メモリ確保 */
    a = (double**)malloc(n*sizeof(double*));
    c = (double**)malloc(n*sizeof(double*));
    for( i = 0; i < n; i++ ){
        a[i] = (double*)malloc(n*sizeof(double));
        c[i] = (double*)malloc(n*sizeof(double));
    }

    /* 連立方程式の入力 */
	kaisu=0;
    for( i = 0; i < n; i++ ){
        for( j = 0; j < n; j++ ){
			a[i][j] = xx[kaisu];
			kaisu++;
//            scanf( "%lf", &a[i][j] );
        }
    }

    /* LU分解 */
    for( i = 0; i < n; i++ ){
        for( j = i+1; j < n; j++ ){
            a[j][i] /= a[i][i];
            for( k = i+1; k < n; k++ ){
                a[j][k] -= a[i][k] * a[j][i];
            }
        }
    }

    /* 逆行列を求める */
    for( k = 0; k < n; k++ ){
        /* 初期化 */
        for( i = 0; i < n; i++ ){
            if( i == k ){ c[i][k] = 1; }
            else{ c[i][k] = 0; }
        }
        /* 解を求める */
        for( i = 0; i < n; i++ ){
            for( j = i+1; j < n; j++ ){
                c[j][k] -= c[i][k] * a[j][i];
            }
        }
        for( i = n-1; i >= 0; i-- ){
            for( j = i+1; j<n; j++ ){
                c[i][k] -= a[i][j] * c[j][k];
            }
            c[i][k] /= a[i][i];
        }
    }

    /* 出力 */
	kaisu=0;
    for( i = 0; i < n; i++ ){
        for( j = 0; j < n; j++ ){
//            printf( "%10lf ", c[i][j] );
			xx[kaisu]=c[i][j];
			kaisu++;
        }
//        putchar( '\n' );
    }
}