#include <stdio.h>

#define DIM 3

void intercambiar_filas(float[DIM][DIM+1], int, int);
void dividir_fila(float[DIM][DIM+1], int, float);
void restar_filas(float[DIM][DIM+1], int, int);
void print_matrix(float [DIM][DIM+1]);

int main(int argc, char **argv)
{
    float matrix[DIM][DIM+1]={{1,3,2,-19},{-5,-1,-5,10},{-1,3,1,-20}};

    int f, c, i;

    printf("Partimos de esta matriz:\n");
    print_matrix(matrix);
    printf("Let's go...\n");

    for(c=0; c<DIM; c++)
    {
        /* Si el elemento pivote ([0][0], [1][1], etc...) es cero, buscamos una fila donde no lo sea y las intercambiamos (si no lo encontramos, entonces toda la columna es cero y quiere decir que el sistema está sobredimensionado) */
        f=0;
        while(matrix[c][c]==0 && matrix[f][c]==0 && f<DIM)
            f++;

        if(f>0 && f<DIM)
        {
            printf("\n--> Intercambiando filas [%d] y [%d]\n", c+1, f+1);
            intercambiar_filas(matrix, c, f);
            print_matrix(matrix);
        }

        /* Convertimos el pivote en 1 */
        printf("\n--> Diviendo fila [%d] entre <%.2f>...\n", c+1, matrix[f][c]);
        dividir_fila(matrix, c, matrix[c][c]);
        print_matrix(matrix);

        /* Y ahora nos cargamos los elementos por encima y por debajo de los pivotes */
        for(f=0; f<DIM; f++)
        {
            if(matrix[f][c]!=0 && matrix[f][c]!=1)
            {
                printf("\n--> Diviendo fila [%d] entre <%.2f>...\n", f+1, matrix[f][c]);
                dividir_fila(matrix, f, matrix[f][c]);
                print_matrix(matrix);
            }

            if(f!=c && matrix[f][c]!=0)
            {
                printf("\n--> Restando fila [%d] menos [%d]...\n", f+1, c+1);
                restar_filas(matrix, f, c);
                print_matrix(matrix);
            }
        }

        printf("\n");
    }

    printf("Y la matriz escalonada es:\n");
    print_matrix(matrix);
    printf("\n");

    printf("X = %.2f\n", matrix[0][DIM]/matrix[0][0]);
    printf("Y = %.2f\n", matrix[1][DIM]/matrix[1][1]);
    printf("Z = %.2f\n", matrix[2][DIM]/matrix[2][2]);

    return 0;
}

void dividir_fila(float matrix[DIM][DIM+1], int fila, float divisor)
{
    int c;

    for(c=0; c<DIM+1; c++)
        matrix[fila][c]/=divisor;
}

void intercambiar_filas(float matrix[DIM][DIM+1], int fila1, int fila2)
{
    float aux;
    int c;

    for(c=0; c<DIM+1; c++)
    {
        aux=matrix[fila1][c];
        matrix[fila1][c]=matrix[fila2][c];
        matrix[fila2][c]=aux;
    }
}

void restar_filas(float matrix[DIM][DIM+1], int f_restando, int f_restador)
{
    int c;

    for(c=0; c<DIM+1; c++)
        matrix[f_restando][c]-=matrix[f_restador][c];
}

void print_matrix(float matrix[DIM][DIM+1])
{
    int f,c;

    printf("***********************************");

    for(f=0; f<DIM; f++)
        for(printf("\n"), c=0; c<DIM+1; c++)
            printf("%8.2f ", matrix[f][c]);

    printf("\n***********************************\n");
}