Simple LU factorization examples in Python

lu_fact.py

import numpy as np

def LU_recursive(A):
    """LU factorization without pivoting. Will fail if any diagonal is 0.
    Recursive implementation.
    """
    L = np.zeros_like(A,dtype=np.float_)
    U = np.eye(*A.shape,dtype=np.float_)
    
    a = A[0,0]
    
    L[0,0] = 1.0
    U[0,0] = a
    
    if A.shape[0] > 1:
        L[1:,0] = A[1:,0]
        U[0,1:] = A[0,1:]
        
        L[1:,0] *= np.power(a,-1.0)
        L[1:,1:] = A[1:,1:] - np.outer(L[1:,0],A[0,1:])
        
        L[1:,1:], U[1:,1:] = LU_recursive(L[1:,1:])
    return L,U
  
  def LU_loop(A):
    """LU factorization without pivoting. Will fail if any diagonal is 0.
    Loop implementation.
    """
    L = np.copy(A)
    U = np.eye(*A.shape,dtype=np.float_)
    
    i = 0
    M = A.shape[0]
    while( True ):

        a = L[i,i]

        L[i,i] = 1.0
        U[i,i] = a

        if(i == M-1):
            break
        else:
            
            U[i,i+1:] = L[i,i+1:]
            L[i,i+1:] = 0.0
            
            L[i+1:,i] *= np.power(a,-1.0)
            U[i+1:,i] = 0.0
            

            L[i+1:,i+1:] -= np.outer(L[i+1:,i],U[i,i+1:])
        i=i+1

    return L,U