/* fichero resolve.c */
#include <stdio.h>
#define   N     10
#define   EPS   1.e-12
void main(void)
{
     int    n, v_ret;
     double a[N][N], x[N], b[N];
     void   lect(int *n, double a[][N], double b[]);
     int    triang(int n, double (*a)[N], double *b);
     void   sust(int n, double a[][N], double b[], double x[]);
     void   escribir(int n, double x[]); 
     
     /* Lectura de datos */
     lect(&n, a, b);
     /* Triangularizacion de la matriz */
     v_ret = triang(n, a, b);
     if(v_ret != 1) {
          printf("Ha aparecido un pivot nulo o muy pequeño. Agur! \n");
          exit(1);
     }
     /* Vuelta atras */
     sust(n, a, b, x);
     /* Escritura de resultados */
     escribir(n, x); 
}

void lect(int *n, double mat[][N], double rhs[]) 
{
     int i, j;
     printf("Numero de ecuaciones: ");
     scanf("%d", n);
     printf("\nMatriz del sistema:\n ");
     for (i=0; i<*n; i++)
          for(j=0; j<*n; j++) {
               printf("\na(%d, %d): ", i+j, j+1);
               scanf("%lf", &mat[i][j]);
          }
     printf("\nTermino independiente:\n ");
     for (i=0; i<*n; i++) {
          printf("\nb(%d): ", i+1);
          scanf("%lf", &rhs[i]);
     }
     return;
}

int triang(int nec, double a[][N], double b[]) 
{
    int i, j, k, error;
    double fac;
    double va( double );
    for(k=0; k<nec-1; k++)
        for(i=k+1; i<nec; i++) {
            if(va(a[k][k]) < EPS) 
                return error = -1;
            fac = -a[i][k]/a[k][k];
            for(j=k; j<nec; j++)
                a[i][j] += a[k][j]*fac;
            b[i] += b[k]*fac;
       }
    return error=1;
}
void sust(int n, double a[][N], double b[], double sol[]) 
{
    int i, k;
    double sum;
    for(k=n-1; k>=0; k--) {
        sum=0.0;
        for(i=k+1; i<n; i++)
            sum += a[k][i]*sol[i];
        sol[k] = (b[k]-sum)/a[k][k];
    }
}

void escribir(int n, double solucion[]) 
{
    int i;
    printf("El vector solucion es: \n");
    for(i=0; i<n; i++)
         printf("solucion(%d)= %lf\n", i, solucion[i]);
}
double va( double u ) 
{
    if( u<0.0 )
        return -u;
    else
        return u;
}