#!/usr/bin/python3
# La décomposition LU   
   
import numpy as np

def matriceP(n, i, j): # P_ij est la matrice identité mais avec les lignes i et j échangées
    P = np.identity(n)
    P[i, i] = P[j, j] = 0
    P[i, j] = P[j, i] = 1
    return P
    
def matriceF(n, j, r): # F est la matrice identité avec les éléments r_i au-dessous de F_jj 
    F = np.identity(n)   
    F[j+1:, j] = r 
    return F
    
def decomp_LU(A):  # effectue la décomposition LU de la matrice A; renvoie P, L et U
    n = A.shape[0]     # le nombre de lignes (et colonnes)
    U = np.copy(A)     # au début, U est A et P et L sont l'identité
    P, L = np.identity(n), np.identity(n)
    for i in range(n-1): # boucle sur toutes les colonnes (sauf la dernière)
        for k in range(i, n): # trouver un pivot...
            if U[k, i] != 0:        # ...cet élément peut l'être ? si oui:
                U[[i, k], :] = U[[k, i], :]   # échanger lignes i et k...
                P = matriceP(n, k, i) @ P     # ...enrégistrer la transformation dans P...
                L =  L @ matriceP(n, k, i)    # ...et dans L...
                pivot = U[i, i]               # ...mémoriser la valeur du pivot...
                break                         # ...et arrêter la boucle sur k
        else:                   # toute la colonne était zéro ?
            continue            #   alors rien à faire pour celle-ci
        r = np.empty(n-i-1) # vecteur de colonne pour construire la matrice F^-1
        for k in range(i+1, n): # éliminer tous les éléments au-dessous de la diagonale:
            facteur = -U[k, i]/pivot        # ajouter (facteur) * (ligne du pivot)...
            U[k, :] += facteur * U[i, :]    # ...à la ligne k...
            r[k-i-1] = -facteur             # ... et mémoriser le facteur dans r
        L = L @ matriceF(n, i, r)     # enrégistrer la transformation dans L
    return P, P @ L, U  # à la fin il faut toujours multiplier L par P à la gauche
 
# tester le code
if __name__ == "__main__":    
    A = np.array([[1,2,3], [-5,2,1], [0,1,7]], dtype=float)
    print("A = ", A)
    P, L, U = decomp_LU(A)
    print("P = ", P)
    print("L = ", L)
    print("U = ", U)
    print("LU = ", L @ U)
    print("PA = ", P @ A)
    print("\n")
    
    A = np.random.randint(10, size=(5,5)).astype(float)
    P, L, U = decomp_LU(A)
    print("P = ", P)
    print("L = ", L)
    print("U = ", U)
    print("LU-PA = ", L @ U - P @ A)