/**
 * @file matriceTest.c
 * @brief Contients l'implémentation de la méthode de Gauss.
 * @author Valentin PORTAIL et Jérémie VILLEPREUX
 * @version 1.0
 * @date 02/10/2022
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"matriceTest.h"



/** \def VECTEUR_COLONNE
 * @brief Correspond à l'indice des colonnes d'un vecteur colonne.
 */
#define VECTEUR_COLONNE 0


/** \def TAILLE_MATRICE
 * @brief Correspond à la taille de la matrice le plus souvent considérer.
 */
#define TAILLE_MATRICE 3



/** 
 * @brief Permet de calculer la valeur absolue d'un flottant.
 *
 * @param[in] nbr      La valeur absolue du nombre à calculer.
 *
 * @return Renvoie la valeur absolue du nombre passé en argument.
 */
double valeurAbsolue(double nbr){

    if (nbr<0){
	return -1*nbr;
    }
    return nbr;
}


/** 
 * @brief Permet de permuter 2 lignes d'une même matrice.
 *
 * @param[in] matrice      Matrice à afficher.
 * @param[in] taille       Correspond à la taille de la matrice.
 * @param[in] Li           Correspond à la ligne à permuter avec **Lj**.
 * @param[in] Lj           Correspond à la ligne à permuter avec **Li**.

 * @return Renvoie 0 si la focntion se déroule normalement  jusqu'à la fin.
 */
int permutationLigne(double ** matrice, unsigned int taille, unsigned int Li, unsigned int Lj){

    double * ligneTempo;

    if (Li<taille && Lj<taille){
	
	ligneTempo=matrice[Li];
	matrice[Li]=matrice[Lj];
	matrice[Lj]=ligneTempo;	
    }

    return EXIT_SUCCESS;
}


/** 
 * @brief Fonction qui appliquer la méthode de Gauss.
 *
 * @param[in] A            Matrice du système (qui sera triangulariser).
 * @param[in] B            Matrice associée au second membre du système.
 * @param[in] n            Correspond à la taille de **A** (matrice carrée).
 * @param[in] epsilon      Correspond à la valeur que l'on se fixe, dans la méthode du **pivot maximal**, pour considérer un pivot presque nul.

 * @return Renvoie la matrice solution du sytème, **NULL** si problème.
 */
double ** methodeDeGauss(double **A, unsigned int n, double **B, double
epsilon){


    if (A==NULL || B==NULL){
        printf("Problème lors de l'allocation dynamique de A ou B, veuillez rééssayer\n");
	exit(-1);
    }

    for (unsigned int k=0; k<n-1; ++k){  //parcours par rapport aux colonnes

    //Stratégie choix de pivot maximal
	if (valeurAbsolue(A[k][k]) < epsilon){ // avec un epsilon aussi petit que possible, on a un pivot preque nul

	    unsigned int indicePivotMax=k;
	    double pivotMaxAbs = valeurAbsolue(A[k][k]); //Car à cette étape, k peut aussi bien etre considerer ligne comme colonne (élément de la diagonale)

	    for (unsigned int l=k+1; l<n; ++l){ //Recherche du pivot max (en valeur absolue)
		if (valeurAbsolue(A[k][l]) > pivotMaxAbs && valeurAbsolue(A[k][l])>=epsilon){
		    indicePivotMax=l;
		    pivotMaxAbs=valeurAbsolue(A[k][l]);
		}
	    }

	    if (indicePivotMax==k){
		/* On n'a pas trouver un pivot valide: donc le système ne peut pas être résolu à l'aide de la méthode de Gauss.
		Le système à soit alors une infinté de solution (si 2 lignes sont proportionnels), 
		soit il n'y a de solution (si le système est incompatible)*/
		
		printf("La matrice A ne peut pas être ramener à une matrice parfaitement triangulaire. Le système admet alors soit une infinité de solution, soit aucune solution.\n ");
		return NULL;
	    }

	    permutationLigne(A, n, indicePivotMax, k); //sinon un pivot max existe
	    
	}


	
	for (unsigned int i=k+1; i<n; ++i){  //parcours par rapport aux lignes
	    
	    double pivot=A[i][k]/A[k][k];

	    for (unsigned int j=k; j<n; ++j){
		A[i][j]=A[i][j] - pivot*A[k][j];
	    }
	    
	    B[i][VECTEUR_COLONNE]=B[i][VECTEUR_COLONNE] -pivot*B[k][VECTEUR_COLONNE];
	}
    }

    double **X=creationMatriceDynamique(n, VECTEUR_COLONNE); //Vecteur colonne solution du système

    if (X==NULL){ //Si l'allocation de X a échoué
	return NULL;
    }


    else{
	
	X[n-1][VECTEUR_COLONNE] = B[n-1][VECTEUR_COLONNE]/A[n-1][n-1];
    
	for (int i=n-2; i>=0; --i) {
            double somme=0;
            for (unsigned int j=i+1; j<n; ++j) {
                somme=somme+A[i][j]*X[j][VECTEUR_COLONNE];
            }
	    X[i][VECTEUR_COLONNE] = (B[i][VECTEUR_COLONNE]-somme)/A[i][i];
        }

	
    }

    

    return X;
    
}
