#!/usr/bin/python3
#
# Calcul de la matrice inverse avec la décomposition LU

import numpy as np
import LU  # LU.py doit contenir une fonction decomp_LU qui renvoie P, L et U

def resoudre(P, L, U, b): # trouve x dans LUx = Pb
    n = U.shape[0] # les matrices sont n fois n
    bprime = P.dot(b)      # le vecteur b' = Pb
    y = np.zeros(n)     # le vecteur y
    for i in range(n):  # initialiser y
        sig = 0
        for k in range(i): # k entre 0 et i-1
            sig += L[i, k] * y[k]
        y[i] = (bprime[i] - sig)/L[i, i]            
    x = np.zeros(n)     # le vecteur x
    for i in range(n-1, -1, -1): # i entre 0 and n-1,
        sig = 0                   #          à l'envers
        for k in range(i+1, n):   # k entre i+1 et n-1
            sig += U[i, k] * x[k]
        x[i] = (y[i] - sig)/U[i, i]
    return x
        
        
def inverse(A): # trouve l'inverse de A
   P, L, U = LU.decomp_LU(A) 
   n = A.shape[0]
   Ainv = np.zeros((n, n))  # la matrice inverse
   iden = np.identity(n)    # la matrice identité dont la i-ème ligne est le vecteur unitaire e_i
   for i in range(n):
       Ainv[:, i] = resoudre(P, L, U, iden[i])
   return Ainv


# Tester le code
if __name__ == '__main__':
    A = np.array([[1,2,3], [-5,2,1], [0,1,7]], dtype=float)
    print(inverse(A))
    print(A @ inverse(A)) 
    
    import numpy.random as rnd
    B = rnd.rand(20, 20) # matrice aléatoire 20 x 20
    nearzero = B @ inverse(B) - np.identity(20)
    print(np.amax(np.abs(nearzero)))