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

#define eps 1e-5

double SolveSystem(int dim, double* matr, double* fr, double* result)
{
  double maxEl;
  double koef, det=1;

  int i , j, k, m;

  for(k=0; k<(dim-1); k++)
  {
    // search of line with max element
    maxEl = fabs( matr[k*dim + k] );
    m = k;
    for(i=k+1; i<dim; i++)
    {
      if(maxEl < fabs(matr[i*dim + k]) )
      {
        maxEl = matr[i*dim + k];
        m = i;
      }
    }

    if(m != k)
    {
      det *= -1;
      for(i=k; i<dim; i++)
      {
        koef                = matr[k*dim + i];
	matr[k*dim + i] = matr[m*dim + i];
	matr[m*dim + i] = koef;
      }
      koef = fr[k];
      fr[k] = fr[m];
      fr[m] = koef;
    }

    if( fabs(matr[k*dim + k]) < eps) return 0; //Virodhz matr


    for(j=(k+1); j<dim; j++)
    {
      koef = - matr[j*dim + k] / matr[k*dim + k];
      for(i=k; i<dim; i++)
      {
	matr[j*dim + i] = matr[j*dim + i] + koef*matr[k*dim + i];
      }
      fr[j] = fr[j] + koef*fr[k]; // free member recalculation
    }
  }

  for(k=(dim-1); k>=0; k--)
  {
    result[k] = fr[k];
    for(i=(k+1); i<dim; i++)
    {
      result[k] -= (matr[k*dim + i]*result[i]);
    }
    result[k] = result[k] / matr[k*dim + k];
  }

  for(k=0; k<dim; k++)
  {
    det *= matr[k*dim+k];
  }

  return det;

}

#define nDim 5

int main(void)
{ //                    clrscr();
  double matr[nDim*nDim] =
  {
  6.37, 1.24, 0.83, 1.225,0.09,
  0.96, 4.05, 1.3 ,-1.63, 0.72,
  1.11,-2.46, 5.55, 2.1,  1.183,
  1.255,0.16, 2.1,  5.99, -14,
  0.39,-1.38,-0.917,7,    -6
  };
  double m_back[nDim*nDim], odinichn[nDim*nDim];
  double fr[nDim] = {2.1, 0.12, 3.01, 7.36, -0.03};
  double additional[nDim], result[nDim];
  double res;
  int i,j,k;

  res = SolveSystem(nDim, matr, fr, result);

  for(i=0; i<nDim; i++)
  {
    for(j=0; j<nDim; j++)
      additional[j] = (i==j);

    res = SolveSystem(nDim, matr, additional, result);


    if( fabs(res)<eps )
    {
      puts("Matr virodzh");
      getchar();
      return 1;
    } else {
      for(j=0; j<nDim; j++)
	m_back[j*nDim + i] = result[j];
    }
  }

  res = SolveSystem(nDim, matr, fr, result);

  puts("A: ");
  for( i = 0 ; i < nDim ; i++)
  {
    for( j = 0 ; j < nDim ; j++)
      printf("%.6f ", matr[i*nDim + j]);
    puts("");
  }


  if( fabs(res) < eps)
  {
    printf("No solution!\n");
    getchar();
    return 1;
  }
  else
  {
    puts("\nA^(-1):");
    for( i = 0 ; i < nDim; i++)
    {
      for( j = 0 ; j < nDim; j++)
	printf("%.6f ", m_back[i*nDim +j]);
      puts("");
    }

    printf("\ndetA = %.6f\n", res);

    printf("Solution:\n");
    for(i=0; i<nDim; i++)
    {
      printf("%.6f ", result[i]);
    }
    puts("");
  }

  for ( i = 0 ; i<nDim ; i++)
    for ( j = 0 ; j < nDim ; j++)
    {
      odinichn[i*nDim + j] = 0;
      for( k = 0 ; k<nDim ; k++)
	odinichn[i*nDim + j] += matr[i*nDim + k] * m_back[k*nDim + j];
    }

  puts("\nA * A^(-1)");
  for( i = 0 ; i < nDim; i++)
  {
    for( j = 0 ; j < nDim; j++)
      printf("%.6e ", odinichn[i*nDim +j]);
    puts("");
  }

  //nevyazki
  for( i = 0 ; i<nDim; i++){
    additional[i] = fr[i];
    for ( j = 0 ; j<nDim ; j++)
    {
    additional[i] -= matr[i*nDim + j] * result[j];
    }
  }
  printf("\nNevyazki:\n");
  for(i=0; i<nDim; i++)
  {
    printf("%.6e ", additional[i]);
  }
  puts("");
  getchar();
  return 0;
}