#!/usr/bin/env python
# coding: utf-8

# # LU Factorization

# In[2]:


import numpy as np
import numpy.linalg as la

# ## Part 1: One Column of LU

# In[3]:


n = 3

np.random.seed(15)
A = np.round(5*np.random.randn(n, n))

A

# Initialize `L` and `U` with zeros:

# In[4]:


L = np.zeros((n,n))
U = np.zeros((n,n))

# Set `U` to be the first row of `A`:

# In[5]:


U[0,:] = A[0,:]
U

# Compute the first column of `L`:

# In[6]:


L[:,0] = A[:,0]/U[0,0]
L

# Compare what we have to `A`:

# In[7]:


print(A)
print(L@U)

# Perform the Schur complement update and store the result in `A1`:

# In[8]:


A1 = A - L @ U
A1

# Take the second row of `U` to be the second row of `A1`:

# In[9]:


U[1,1:] = A1[1,1:]
U

# We can now compute the next column of `L`:

# In[10]:


L[1:,1] = A1[1:,1]/U[1,1]
L

# And finally, compute the bottom right elements of `L` and `U`

# In[11]:


U[2,2] = A1[2,2] - L[2,1]*U[1,2]
L[2,2] = 1.0

# In[12]:


print(L)
print(U)

# In[13]:


print(A)
print(L@U)

# ## Part 2: The Full Algorithm
# 
# Implement the general LU factorization algorithm

# In[14]:


n = 4
A = np.random.random((n,n)) 
L = np.zeros((n,n)) 
U = np.zeros((n,n)) 
M = A.copy()

# In[15]:


for i in range(n):
    U[i,i:] = M[i,i:]
    L[i:,i] = M[i:,i]/U[i,i]
    M[i+1:,i+1:] -= np.outer(L[i+1:,i:i+1],U[i:i+1,i+1:])   

# In[16]:


print(L)
print(U)
print(A-L@U)

# ## Part 3: LU with Partial Pivoting
# 
# When a divisor (U_ii) becomes zero, while nonzeros remain in the trailing matrix, the LU factorization does not exist. Further, if U_ii is small, the magnitude of elements in L and U can blow up, increasing round-off error. Partial pivoting circumvents both problems.

# In[31]:


n = 4
A = np.random.random((n,n)) 
L = np.zeros((n,n)) 
U = np.zeros((n,n)) 
P = np.eye(n,n)
M = A.copy()

def perm(i,j):
    P = np.eye(n,n)
    P[i,i] = 0
    P[j,j] = 0
    P[i,j] = 1
    P[j,i] = 1
    return P

for i in range(n):
    imax = np.argmax(np.abs(M[i:,i]))
    Pi = perm(i,i+imax)
    P = Pi @ P
    L = Pi @ L
    M = Pi @ M
    U[i,i:] = M[i,i:]
    L[i:,i] = M[i:,i]/U[i,i]
    M[i+1:,i+1:] -= np.outer(L[i+1:,i:i+1],U[i:i+1,i+1:])

# Check residual error for the permuted A, i.e., PA-LU. What property do off-diagonal entries in L satisfy?

# In[32]:


print(P@A - L @ U)

# In[33]:


L

# In[34]:


P

# In[ ]: