# coding: utf-8

# # LU with Partial Pivoting

# In[131]:

import numpy as np
import numpy.linalg as la

np.set_printoptions(precision=3, suppress=True)


# # Set-up

# Let's grab a (admittedly well-chosen) sample matrix `A`:

# In[132]:

n = 4

np.random.seed(235)
A = np.round(5*np.random.randn(n, n))
A[0,0] = 0
A[2,1] = 17
A[0,2] = 19
A


# ----------------
# ## Permutation matrices
# 
# Now define a function `row_swap_mat(i, j)` that returns a permutation matrix that swaps row i and j:

# In[133]:

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


# What do these matrices look like?

# In[134]:

row_swap_mat(0,1)


# Do they work?

# In[135]:

row_swap_mat(0,1).dot(A)


# --------------

# # Part I

# `U` is the copy of `A` that we'll modify:

# In[136]:

U = A.copy()


# ## First column

# Create P1 to swap up the right row: 

# In[137]:

P1 = row_swap_mat(0, 3)
U = P1.dot(U)
U


# In[138]:

M1 = np.eye(n)
M1[1,0] = -U[1,0]/U[0,0]
M1[2,0] = -U[2,0]/U[0,0]
M1


# In[139]:

U = M1.dot(U)
U


# ## Second column
# 
# Create `P2` to swap up the right row:

# In[140]:

P2 = row_swap_mat(2,1)
U = P2.dot(U)
U


# Make the second-column elimination matrix `M2`:

# In[141]:

M2 = np.eye(n)
M2[2,1] = -U[2,1]/U[1,1]
M2[3,1] = -U[3,1]/U[1,1]
M2


# In[142]:

U = M2.dot(U)
U


# ## Third column
# 
# Create `P3` to swap up the right entry:

# In[143]:

P3 = row_swap_mat(3, 2)
U = P3.dot(U)
U


# Make the third-column elimination matrix `M3`:

# In[144]:

M3 = np.eye(n)
M3[3,2] = -U[3,2]/U[2,2]
M3


# In[145]:

U = M3.dot(U)
U


# ## Wrap-up

# So we've built $M3P_3M_2P_2M_1P_1A=U$.

# In[150]:

M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1).dot(A)


# ---------------------
# That left factor is anything but lower triangular:

# In[151]:

M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1)


# # Part II

# Now try the reordering trick:

# In[160]:

L3 = M3
L2 = P3.dot(M2).dot(la.inv(P3))
L1 = P3.dot(P2).dot(M1).dot(la.inv(P2)).dot(la.inv(P3))


# In[155]:

L3.dot(L2).dot(L1).dot(P3).dot(P2).dot(P1)


# --------------
# We were promised that all of the `L`*n* are still lower-triangular:

# In[168]:

print(L1)
print(L2)
print(L3)


# So their product is, too:

# In[172]:

Ltemp = L3.dot(L2).dot(L1)
Ltemp


# ----
# `P` is still a permutation matrix (but a more complicated one):

# In[174]:

P = P3.dot(P2).dot(P1)
P


# -----------------
# So to sum up, we've made:

# In[175]:

Ltemp.dot(P).dot(A) - U


# --------------
# Multiply from the left by `Ltemp`${}^{-1}$, which is *also* lower triangular:

# In[179]:

L = la.inv(Ltemp)
L


# In[180]:

P.dot(A) - L.dot(U)