#include <stdio.h>
#include <stdlib.h>
#define TAILLE 20

struct matrice {
    double** values;
    int height; // i
    int length; // j
};

typedef struct matrice matrice;

/* fonctions annexes */
matrice* init(int hauteur, int largeur)
{
    int i,j;
    matrice* res=malloc(sizeof(matrice));
    res->height=hauteur;
    res->length=largeur;
    res->values=malloc(hauteur*sizeof(double*));
    for (i=0;i<hauteur;i++)
    {
        res->values[i]=malloc(largeur*sizeof(double));
        for (j=0;j<largeur;j++)
        {
            res->values[i][j]=0;
        }
    }
    return res;
}
/* fonction d'affichage de matrice de taille height*length */
void affichage_matrice(matrice* mat)
{
    int i,j;
    for (i=0;i<mat->height;i++)
    {
        for (j=0;j<mat->length;j++)
        {
            printf("%.3f ",mat->values[i][j]);
        }
        printf("\n");
    }
}


matrice* copier_matrice(matrice* mat)
{
    int i,j;
    matrice* res = init(mat->height,mat->length);
    for(i=0;i<res->height;i++)
    {
        for(j=0;j<res->length;j++)
        {
            res->values[i][j]=mat->values[i][j];
        }
    }
    return res;
}

/** Exercice 1.a
    La matrice A doit être triangulaire inférieure */
matrice* solinf(matrice* L,matrice* b)
{
    matrice* res=init(L->height,1);
    int i,j;
    int error=0;
    if (b->height != L->height || b->length != 1) {error = 1;}
    else
    {
        if (L->height != L->length) {error = 2;}
        else
        {
            for (i=0;i<L->height;i++)
            {
                for (j=i+1;j<L->length;j++)
                {
                    if (L->values[i][j] != 0) error=3;
                }
            }
        }
    }
    switch(error)
    {
        case 1:
          fprintf(stderr,"[Erreur] la matrice n'est pas carrée\n");
          break;
        case 2:
          fprintf(stderr,"[Erreur] La matrice b n'a pas les bonnes dimensions\n");
          break;
        case 3:
          fprintf(stderr,"[Erreur] la matrice n'est pas triangulaire inférieure\n");
          break;
    }

    // La matrice ne doit pas contenir de 0 sur la diagonale
    for (i=0;i<L->height;i++)
    {
        res->values[i][0]=0;
        for (j=0;j<i;j++)
        {
            res->values[i][0]+=L->values[i][j]*res->values[j][0];
        }
        res->values[i][0]=(b->values[i][0]-res->values[i][0])/L->values[i][i];
    }
    return res;
}

/** Exercice 1.b
    La matrice A doit être triangulaire supérieure */
matrice* solsup(matrice* U,matrice* b)
{
    matrice* res=init(U->height,1);
    int i,j;
    int error=0;
    if (b->height != U->height || b->length != 1) {error = 1;}
    else
    {
        if (U->height != U->length) {error = 2;}
        else
        {
            for (j=0;j<U->length;j++)
            {
                for (i=j+1;i<U->height;i++)
                {
                    if (U->values[i][j] != 0) error=3;
                }
            }
        }
    }
    switch(error)
    {
        case 1:
          fprintf(stderr,"[Erreur] la matrice n'est pas carrée\n");
          break;
        case 2:
          fprintf(stderr,"[Erreur] La matrice b n'a pas les bonnes dimensions\n");
          break;
        case 3:
          fprintf(stderr,"[Erreur] la matrice n'est pas triangulaire supérieure\n");
          break;
    }

    // La matrice ne doit pas contenir de 0 sur la diagonale
    for (i=U->height-1;i>=0;i--)
    {
        res->values[i][0]=0;
        for (j=i+1;j<U->height;j++)
        {
            res->values[i][0]+=U->values[i][j]*res->values[j][0];
        }
        res->values[i][0]=(b->values[i][0]-res->values[i][0])/U->values[i][i];
    }
    return res;
}

/* Exercice 2.a */
matrice** trigGauss(matrice* A, matrice* b)
{
    matrice** tab = (matrice**)malloc(2*sizeof(matrice*));
    matrice* mat = copier_matrice(A);
    matrice* vect = copier_matrice(b);
    int i,j,k,c;
    for(k=0;k<mat->height-1;k++)
    {
        if (mat->values[k][k]==0) printf("[Erreur] Une des valeurs sur la diagonale est nulle, utilisez la méthode du pivot\n");
        for(i=k+1;i<mat->height;i++)
        {
            c=mat->values[i][k]/mat->values[k][k];
            vect->values[i][0]=vect->values[i][0]-c*vect->values[k][0];
            mat->values[i][k]=0;
            for(j=k+1;j<mat->length;j++)
            {
                mat->values[i][j]=mat->values[i][j]-c*mat->values[k][j];
            }
        }
    }
    tab[0]=mat;
    tab[1]=vect;
    return tab;
}

/* Exercice 2.b */
matrice* ResolutionGauss(matrice* A, matrice* b)
{
    matrice** gauss=trigGauss(A,b);
    matrice* res=solsup(gauss[0],gauss[1]);
    return res;
}

/* Exercice 3.a */
matrice** FactLU(matrice* A)
{
    matrice** tab = (matrice**)malloc(2*sizeof(matrice*));
    matrice* L=init(A->height,A->length);
    matrice* U=init(A->height,A->length); // Initialisation à la matrice nulle
    int i,j,k;
    double sum;
    for(i=0;i<A->height;i++)
    {
        L->values[i][i]=1; // Initialisation à la matrice identitée
    }
    for(i=0;i<A->height-1;i++)
    {
        for(j=i;j<A->height;j++)
        {
            sum = 0;
            for(k=0;k<i+1;k++)
            {
                sum=sum+L->values[i][k]*U->values[k][j];
            }
            U->values[i][j]=A->values[i][j]-sum;
        }
        for(j=i+1;j<A->height;j++)
        {
            sum = 0;
            for(k=0;k<i+1;k++)
            {
                sum=sum+L->values[j][k]*U->values[k][i];
            }
            L->values[j][i]=(A->values[j][i]-sum)/U->values[i][i];
        }
    }
    sum=0;
    for(k=0;k<A->height-1;k++)
    {
        sum=sum+L->values[L->height-1][k]*U->values[k][U->length-1];
    }
    U->values[U->height-1][U->length-1]=A->values[A->height-1][A->length-1] - sum;
    tab[0]=L;
    tab[1]=U;
    return tab;
}

matrice* ResolutionLU(matrice* A, matrice* b)
{
    matrice** tab=FactLU(A);
    matrice* res=init(A->height,1);
    res=solinf(tab[0],b);
    res=solsup(tab[1],res);
    return res;
}

/* Similaire à ResolutionLU, permet de réduire le nombre de calculs pour inverser une matrice */
matrice* ResolutionLU2(matrice** tab,matrice* b)
{
    matrice* res=solinf(tab[0],b);
    res=solsup(tab[1],res);
    return res;
}

matrice* inverse(matrice* A)
{
    matrice* mat;
    matrice* res=init(A->height,A->length);
    matrice* b=init(A->height,1);
    matrice** tab=FactLU(A);
    int i,j;
    for(i=0;i<A->height;i++)
    {
        b->values[i][0]=1; /* Initialisation du vecteur unitaire */
        mat=ResolutionLU2(tab,b);
        for(j=0;j<A->height;j++)
        {
            res->values[j][i]=mat->values[j][0];
        }
        b->values[i][0]=0;
    }
    return res;
}


/* Initialisation 1.a */
matrice* exo1a_L(matrice* mat)
{
    mat->values[0][0]=1;
    mat->values[0][1]=0;
    mat->values[0][2]=0;
    mat->values[1][0]=2;
    mat->values[1][1]=3;
    mat->values[1][2]=0;
    mat->values[2][0]=1;
    mat->values[2][1]=4;
    mat->values[2][2]=-1;
    return mat;
}

matrice* exo1a_b(matrice* mat)
{
    mat->values[0][0]=1;
    mat->values[1][0]=8;
    mat->values[2][0]=10;
    return mat;
}

/* Initialisation 1.b */
matrice* exo1b_U(matrice* mat)
{
    mat->values[0][0]=1;
    mat->values[0][1]=2;
    mat->values[0][2]=3;
    mat->values[1][0]=0;
    mat->values[1][1]=4;
    mat->values[1][2]=8;
    mat->values[2][0]=0;
    mat->values[2][1]=0;
    mat->values[2][2]=5;
    return mat;
}

matrice* exo1b_b(matrice* mat)
{
    mat->values[0][0]=6;
    mat->values[1][0]=16;
    mat->values[2][0]=15;
    return mat;
}

/* Initialisation 2.a */
matrice* exo2a_A(matrice* mat)
{
    mat->values[0][0]=3;
    mat->values[0][1]=1;
    mat->values[0][2]=2;
    mat->values[1][0]=3;
    mat->values[1][1]=2;
    mat->values[1][2]=6;
    mat->values[2][0]=6;
    mat->values[2][1]=1;
    mat->values[2][2]=-1;
    return mat;
}

matrice* exo2a_b(matrice* mat)
{
    mat->values[0][0]=2;
    mat->values[1][0]=1;
    mat->values[2][0]=4;
    return mat;
}

/* Initialisation 2.b */
matrice* exo2b_A(matrice* mat)
{
    mat->values[0][0]=1;
    mat->values[0][1]=2;
    mat->values[0][2]=3;
    mat->values[1][0]=5;
    mat->values[1][1]=2;
    mat->values[1][2]=1;
    mat->values[2][0]=3;
    mat->values[2][1]=-1;
    mat->values[2][2]=1;
    return mat;
}

matrice* exo2b_b(matrice* mat)
{
    mat->values[0][0]=5;
    mat->values[1][0]=5;
    mat->values[2][0]=6;
    return mat;
}
/* Initialisation 3.a */

matrice* exo3a_A(matrice* mat)
{
    mat->values[0][0]=3;
    mat->values[0][1]=1;
    mat->values[0][2]=2;
    mat->values[1][0]=3;
    mat->values[1][1]=2;
    mat->values[1][2]=6;
    mat->values[2][0]=6;
    mat->values[2][1]=1;
    mat->values[2][2]=-1;
    return mat;
}

/* Initialisation 3.b */

matrice* exo3b_A(matrice* mat)
{
    mat->values[0][0]=1;
    mat->values[0][1]=2;
    mat->values[0][2]=3;
    mat->values[1][0]=5;
    mat->values[1][1]=2;
    mat->values[1][2]=1;
    mat->values[2][0]=3;
    mat->values[2][1]=-1;
    mat->values[2][2]=1;
    return mat;
}

matrice* exo3b_b(matrice* mat)
{
    mat->values[0][0]=5;
    mat->values[1][0]=5;
    mat->values[2][0]=6;
    return mat;
}

/* Initialisation 4 */

matrice* exo4_A(matrice* mat)
{
    mat->values[0][0]=1;
    mat->values[0][1]=2;
    mat->values[0][2]=3;
    mat->values[1][0]=5;
    mat->values[1][1]=2;
    mat->values[1][2]=1;
    mat->values[2][0]=3;
    mat->values[2][1]=-1;
    mat->values[2][2]=1;
    return mat;
}

int main()
{
    printf("   Exercice 1.a : matrice triangulaire inférieure\n   Valeur de x attendue : 1 2 -1\n");
    matrice* L = init(3,3);
    L = exo1a_L(L);
    matrice* b = init(3,1);
    b = exo1a_b(b);
    matrice* x = solinf(L,b);
    affichage_matrice(x);

    printf("   Exercice 1.b : matrice triangulaire supérieure\n   Valeur de x attendue : 1 -2 3\n");
    matrice* U = init(3,3);
    U = exo1b_U(U);
    b = exo1b_b(b);
    x = solsup(U,b);
    affichage_matrice(x);

    printf("   Exercice 2.a : Algorithme d'élimination de Gauss\n   Valeurs de A et x attendues :\n  1  0  0     2\nA 0  3  0  x -3\n  0  0 -1     5\n\n");
    matrice* A = init(3,3);
    A = exo2a_A(A);
    b = exo2a_b(b);
    matrice** tab = trigGauss(L,b);
    affichage_matrice(tab[0]);
    affichage_matrice(tab[1]);

    printf("   Exercice 2.b : Algorithme de résolution de Gauss\n   Valeur de x attendue : 16/30 19/30 32/30\n");
    A = exo2b_A(A);
    b = exo2b_b(b);
    x = ResolutionGauss(A,b);
    affichage_matrice(x);

    printf("   Exercice 3.a : Algorithme de factorisation par méthode Doolittle\n   Valeurs de L et U attendues :\n  1  0  0     3  1  2\nL 1  1  0   U 0  1  4\n  2 -1  1     0  0 -1\n\n");
    A = exo3a_A(A);
    tab=FactLU(A);
    printf("Matrice L : \n");
    affichage_matrice(tab[0]);
    printf("Matrice U : \n");
    affichage_matrice(tab[1]);

    printf("   Exercice 3.b : Alorithme de résolution par méthode Doolittle\n   Valeur de x attendue : 1 -1 2\n");
    A = exo3b_A(A);
    b = exo3b_b(b);
    x = ResolutionLU(A,b);
    affichage_matrice(x);

    printf("   Exercice 4 : Inversion de matrice\n   Valeur de A attendue :\n-0.088  0.147  0.117\n 0.059  0.235 -0.411\n 0.323 -0.206  0.235\n\n");
    A = exo4_A(A);
    A = inverse(A);
    affichage_matrice(A);
    return 0;
}