text: import numpy as np
  import copy
  import scipy
  import scipy.linalg
  
  def ResolverCosasProgresivas(a,b):
      a = np.array(a)
      n = len(b)
  # Eliminación
      for j in range(0,n-1): # colunmas
          for i in range(j+1,n): # renglones
              if a[i,j] != 0.0:
                  lam = a[i,j]/a[j,j]
                  L[i][j] = lam
                  a[i,j:n] = a[i,j:n] - lam*a[j,j:n]
                  b[i] = b[i] - lam*b[j]
      for i in range (0, n-1):
          b[i] = (b[i] - np.dot(a[i,i+1:n], b[i+1:n]))/a[i,i]
      return b
  
  
  
  def ResolverCosasRegresivas(a,b):
      a = np.array(a)
      n = len(b)
  # Eliminación
      for j in range(0,n-1): # colunmas
          for i in range(j+1,n): # renglones
              if a[i,j] != 0.0:
                  lam = a[i,j]/a[j,j]
                  L[i][j] = lam
                  a[i,j:n] = a[i,j:n] - lam*a[j,j:n]
                  b[i] = b[i] - lam*b[j]
      # Sustitución regresiva
      for i in range(n-1,-1,-1):
          b[i] = (b[i] - np.dot(a[i,i+1:n],b[i+1:n]))/a[i,i]
      return b
  
  
  
  def LU(a):
      L = [[0 for i in range(len(a))]for j in range(len(a))] #inicializar matriz de ceros para la trinagular inferior
      for i in range(0,len(L)):
          L[i][i] = 1
  
      n = len(a)
  # Eliminación
      for j in range(0,n-1): # colunmas
          for i in range(j+1,n): # renglones
              if a[i,j] != 0.0:
                  lam = a[i,j]/a[j,j]
                  L[i][j] = lam
                  a[i,j:n] = a[i,j:n] - lam*a[j,j:n]
      return L, a
  
  def Imprime(M):
      for l in M:
          print(l)
      print("\n")
  
  if __name__ == "__main__":
      
      a = np.array([[1,-3,7], [0,3,-2], [-2,6,1.]])
      a_original = copy.deepcopy(a)
      a_inv = []
  
      I = np.identity(len(a))
  
      L = []
  
      L, a = LU(a)
      
      Imprime(L)
      Imprime(a)
      print (f"{np.dot(L,a)}")
      
      for i in range(0,len(a)):
          b = I[i]
          y = ResolverCosasProgresivas(L,b)
          x = ResolverCosasRegresivas(a,y)
          a_inv.append(x.T)
  
      a_inv = np.array(a_inv).T
      print(f"A inversa es \n")
      Imprime (a_inv)
  
      print (f"Comprobacion A * A^-1 \n{np.dot(a_original, a_inv)}")
  
  
  
  #%%
  