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

#define N 4

double print(int n, double y[n][n])
{
    int i,j,k;

    for (i=0;i<n;i++)
    {
        for (j=0;j<n;j++)
        {
            printf("%8.3lf",y[i][j]);

            if (j==(n-1))
                printf("\n");
            else 
                printf(" ");
        }
    }
}

double printvec(int n, double x[n])
{
    int i;

    for (i=0;i<N;i++)
    {
        printf("%8.3lf",x[i]);

        if (i==(n-1))
            printf("\n");
        else 
            printf(" ");    
    }
}

int main(void)
{
    double A[N][N]={{2,1,4,1},{3,4,-1,-1},{1,-4,1,5},{2,-2,1,3}}, b[N]={-4,3,9,7}, l[N][N]={0}, y[N]={0}, x[N]={0};
    int i,j,k;

    printf("Antes do escalonamento: \n");
    print(N,A);

    for (i=0;i<N;i++)
    {

        l[i][i] = 1.0;

        for (j=i+1;j<N;j++)
        {
            l[j][i]=A[j][i]/A[i][i];

            for (k=i;k<N;k++)
            {
                A[j][k] = A[j][k] - l[j][i]*A[i][k];
            }
        }
    }

    printf("U= \n");
    print(N,A);

    printf("L= \n");
    print(N,l);

    for (i=0;i<N;i++) // SUBSTITUIÇÃO PROGRESSIVA
    {
        double soma = b[i];

        for (j=0;j<i;j++)
        {
            soma-=l[i][j]*y[j];
        }

        y[i]=soma;
    }

    for (i=N-1;i>=0;i--) // SUBSTITUIÇÃO REGRESSIVA
    {
        double soma = y[i];

        for (j=i+1;j<N;j++)
        {
            soma-=A[i][j]*x[j];
        }

        x[i]=soma/A[i][i];
    }


    printf("Solução: \n");
    printvec(N,x);

    return 0;
}