/* Questo programma risolve un sistema la cui matrice associata 
è non singolare attraverso l'algoritmo dell'eliminazione di Gauss */

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

const int SIZE=51;
typedef float matrix [SIZE][SIZE];
typedef float vector [SIZE];
void errore(int);
void presentazione ();
void leggi_dim (int&);
float modulo(float);
void leggi_mat (int,matrix,vector);
void Gauss (int,matrix,vector);
int find_max (int, matrix, int, int);
void change_lines (int,int,int,matrix,vector);
void solve (int,matrix,vector,vector);
void print_sol (int,vector);


main ()
{
  presentazione ();
  int n; matrix a; vector b,sol;
  leggi_dim (n);
  leggi_mat (n,a,b);
  Gauss (n,a,b);
  solve (n,a,b,sol);
  print_sol (n,sol);
  system ("PAUSE");
  return 0;
}

void presentazione ()
{
  printf ("* * * Questo programma esegue l'algoritmo dell'eliminazione di Gauss * * *\n");
  printf ("\n\n* * * *  su una matrice non singolare di dimensione nxn  * * * *\n\n\n");
  system ("PAUSE");
  return;
}

void leggi_dim (int& dim)
{
  printf ("\nQual e` la dimensione della matrice dei coefficineti associata al sistema?");
  printf("\nInserici la dimensione: ");
  scanf ("%d",&dim);
  while (dim<=0 || dim>=SIZE)
    {
      printf ("\nAttenzione!");
      printf("\nLa dimensione inserita non e' permessa.\n Inserisci dinuovo la dimensione: ");
      scanf ("%d",&dim);
    }
  return;
}

void leggi_mat (int dim,matrix a,vector b)
{
  int i,j;
  printf ("\nInserisci i coefficienti della matrice:");
  for (i=1;i<=dim;i++)
    {
      for (j=1;j<=dim;j++)
	  {
	    printf ("\na[%d,%d]=",i,j); scanf ("%f",&a[i][j]);
	  }
    }
  printf ("\nInserisci i coefficienti del vettore dei termini noti:");
  for (i=1;i<=dim;i++)
    {
      printf ("\nb[%d]=",i); scanf ("%f",&b[i]);
    }
  return;
}

void Gauss (int dim,matrix a,vector b)
{
  int i,j,k,r; float m;
  for (k=1;k<dim;k++)
    {
      for (i=k;i<=dim;i++)
        {
	      r=find_max (dim,a,i,k); if (r!=i) change_lines (dim,i,r,a,b);
        }
      for (i=k+1;i<=dim;i++) 
        {
          if(a[k][k]==0) errore(1);
          m=a[i][k]/a[k][k]; a[i][k]=0;
	      for (j=k+1;j<=dim;j++) a[i][j]=a[i][j]-m*a[k][j];
	      b[i]=b[i]-m*b[k];
        }
     }
  printf ("\n\nMatrice ottenuta dall'eliminazione di Gauss:");
  for (i=1;i<=dim;i++)
    {
      printf ("\n| ");
      for (j=1;j<=dim; j++) printf ("%6.3f\t ",a[i][j]);
      printf ("| %6.3f |",b[i]);
    }
  return;
}

int find_max (int dim,matrix a,int i,int k)
{
  int imax=i;
  for (int j=i+1;j<=dim;j++)
    {
      if (modulo(a[j][k])>modulo(a[imax][k])) imax=j;
    }
  return imax;
}

void change_lines (int dim,int i,int j,matrix a,vector b)
{
  float aux;
  for (int l=1;l<=dim;l++)
    {
      aux=a[i][l]; a[i][l]=a[j][l]; a[j][l]=aux;
    }
  aux=b[i]; b[i]=b[j]; b[j]=aux;
  return;
}

void solve (int dim,matrix a,vector b,vector x)
{
  if(a[dim][dim]==0) errore(2);
  x[dim]=b[dim]/a[dim][dim];
  for (int i=dim-1;i>=1;i--)
    {
      for (int j=i+1;j<=dim;j++) b[i]=b[i]-a[i][j]*x[j];
      x[i]=b[i]/a[i][i];
    }
  return;
}

void print_sol (int dim,vector x)
{
  printf ("\n\nLe soluzioni del sistema sono:");
  for (int i=1;i<=dim;i++)
    printf ("\n\nx[%d]=%f",i,x[i]);
  printf ("\n\n");
  return;
}

void errore (int x)
{
   if(x==1){printf("\nAttenzione, la matrice inserita e' singolare\n");system ("PAUSE");exit(1);}
   if(x==2){printf("\nAttenzione il sistema non e' compatibile\n");system ("PAUSE");exit(1);}
}
float modulo(float a)
{
   if(a>=0) return a;
   else return -a;
}