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

# ChEn-3170: Computational Methods in Chemical Engineering Fall 2018 UMass Lowell; Prof. V. F. de Almeida **24Sep2018**
# 
# # 05. Computational Linear Algebra Fundamentals
# $  
#   \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}}
#   \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}}
#   \newcommand{\Cmtrx}{\boldsymbol{\mathsf{C}}}
#   \newcommand{\Dmtrx}{\boldsymbol{\mathsf{D}}}
#   \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}}
#   \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}}
#   \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}}
#   \newcommand{\Qmtrx}{\boldsymbol{\mathsf{Q}}}
#   \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}}
#   \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}}
#   \newcommand{\xvec}{\boldsymbol{\mathsf{x}}}
#   \newcommand{\avec}{\boldsymbol{\mathsf{a}}}
#   \newcommand{\bvec}{\boldsymbol{\mathsf{b}}}
#   \newcommand{\cvec}{\boldsymbol{\mathsf{c}}}
#   \newcommand{\rvec}{\boldsymbol{\mathsf{r}}}
#   \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert}
#   \DeclareMathOperator{\rank}{rank}
# $

# ---
# ## Table of Contents
# * [Theory](#theory)
# * [Matrix-vector and matrix-matrix product operations](#product)
# * [NumPy and SciPy Linear Algebra](#pylinalg)
#  + [Matrix solve](#pysolve)
#  + [$\Lmtrx$ forward solve](#pyl)
#  + [$\Umtrx$ backward solve](#pyu)
#  + [$\Amtrx = \Pmtrx\,\Lmtrx\,\Umtrx$ factorization](#[pyplu)
# * [Course Linear Algebra](#courselinalg)
#  + [$\Lmtrx$ forward solve](#l)
#  + [$\Umtrx$ forward solve](#u)
#  + [$\Amtrx = \Lmtrx\,\Umtrx$ factorization](#lu)
#  + [$\Pmtrx\,\Amtrx = \Lmtrx\,\Umtrx$ factorization](#plu)
#  + [$\Pmtrx\,\Amtrx\,\Qmtrx = \Lmtrx\,\Umtrx$ factorization](#pqlu)
# ---

# ## Theory<a id="theory"></a>
# The course notes (OneNote [ChEn-3170-linalg](https://studentuml-my.sharepoint.com/:o:/g/personal/valmor_dealmeida_uml_edu/EhHdq8qIW_JNiMBu60tKHi8BV6Nv8jonllo__NOtX7JgkQ?e=jA5xj8)) cover basic elements of lineary system of algebraic equations. Their representation in row and column formats and the insight gained from them when searching for a solution of the system.
# 
# Some basic theoretical aspects of solving for $\xvec$ in the matrix equation $\Amtrx\,\xvec = \bvec$ are covered. $\Amtrx$ is a matrix, $\bvec$ and $\xvec$ are vectors.

# ## Matrix-vector and matrix-matrix product operations<a id="product"></a>
# The following operations between vectors and matrices are obtained directly from the buil-in functions in the `numpy` package.

# In[ ]:


'''Import the NumPy package as usual'''

import numpy as np


# Scalar product of two vectors: $\avec \cdot \bvec$.

# In[ ]:


'''Vector scalar product or dot product of vectors'''

a_vec = np.array( np.random.random(3) )
b_vec = np.array( np.random.random(3) )

a_vec_dot_b_vec = np.dot( a_vec, b_vec ) # clear linear algebra operation
print('a.b =', a_vec_dot_b_vec)

a_vec_x_b_vec = a_vec @ b_vec   # consistent linear algebra multiplication
print('a@b =', a_vec_x_b_vec )


# Matrix vector product: $\Amtrx\,\bvec$.

# In[ ]:


'''Matrix-vector product'''

a_mtrx = np.array( [ [ 2.,  1., 1.],   # per course notes (p 04)
                     [ 4., -6., 0.],
                     [-2.,  7., 2.]])

b_vec = np.array( [5., -2., 9.])       # per course notes

a_mtrx_x_b_vec = a_mtrx @ b_vec # linear algebra matrix-vector product

print('A x b =', a_mtrx_x_b_vec)


# Matrix-vector product: $\Imtrx\,\bvec = \bvec$. Note: $\begin{pmatrix}
# 1 & 0 & 0 \\
# 0 & 1 & 0 \\
# 0 & 0 & 1
# \end{pmatrix} \, \begin{pmatrix} b_1\\b_2\\b_3 \end{pmatrix} = \begin{pmatrix} b_1\\b_2\\b_3 \end{pmatrix} $.

# In[ ]:


'''Identity-matrix vector product'''

i_mtrx = np.eye(3)

i_mtrx_x_b_vec = i_mtrx @ b_vec # linear algebra matrix-vector product

print('I x b =', i_mtrx_x_b_vec)
print('b     =', b_vec)


# Matrix-matrix product: $\Imtrx\,\Amtrx = \Amtrx$. Note: $\begin{pmatrix}
# 1 & 0 & 0 \\
# 0 & 1 & 0 \\
# 0 & 0 & 1
# \end{pmatrix} \, 
# \begin{pmatrix} 
# A_{1,1} & A_{1,2} & A_{1,3} \\
# A_{2,1} & A_{2,2} & A_{2,3} \\
# A_{3,1} & A_{3,2} & A_{3,3}
# \end{pmatrix} = 
# \begin{pmatrix} 
# A_{1,1} & A_{1,2} & A_{1,3} \\
# A_{2,1} & A_{2,2} & A_{2,3} \\
# A_{3,1} & A_{3,2} & A_{3,3}
# \end{pmatrix}
# $.

# In[ ]:


'''Matrix-matrix product IA = A'''

i_mtrx_x_a_mtrx = i_mtrx @ a_mtrx   # linear algebra matrix-matrix product

print('I x A =\n', i_mtrx_x_a_mtrx)
print('A     =\n', a_mtrx)


# Matrix-matrix product: $\Amtrx\,\Bmtrx = \Cmtrx$. Note: 
# $\begin{pmatrix}
# A_{1,1} & A_{1,2} & A_{1,3} \\
# A_{2,1} & A_{2,2} & A_{2,3} \\
# A_{3,1} & A_{3,2} & A_{3,3}
# \end{pmatrix} \, 
# \begin{pmatrix} 
# B_{1,1} & B_{1,2} & B_{1,3} \\
# B_{2,1} & B_{2,2} & B_{2,3} \\
# B_{3,1} & B_{3,2} & B_{3,3}
# \end{pmatrix} = 
# \begin{pmatrix} 
# C_{1,1} & C_{1,2} & C_{1,3} \\
# C_{2,1} & C_{2,2} & C_{2,3} \\
# C_{3,1} & C_{3,2} & C_{3,3}
# \end{pmatrix}
# $ where each $C_{i,j}$ is a vector product of the $i$th row of $\Amtrx$ and the $j$th column of $\Bmtrx$, *i.e.* 
# $C_{i,j} = \sum\limits_{k=1}^3 A_{i,k}\, B_{k,j}$.

# In[ ]:


'''Matrix-matrix product AB = C'''

b_mtrx = np.array( [[5. , 5. , 5.],
                    [-2., -2., -2.],
                    [9. , 9. , 9.]]
                 )
c_mtrx = a_mtrx @ b_mtrx   # linear algebra matrix-matrix product

print('A =\n', a_mtrx)
print('B =\n', b_mtrx)
print('C =\n', c_mtrx)


# <div class="alert alert-block alert-danger">
# The matrix-matrix product: $\Bmtrx\,\Amtrx = \Dmtrx \ne \Cmtrx$, does not commute in general. Note 
# $D_{i,j} = \sum\limits_{k=1}^3 B_{i,k}\, A_{k,j}$.
# </div>

# In[ ]:


'''Matrix-matrix product BA = D'''

b_mtrx = np.array( [[5. , 5. , 5.],
                    [-2., -2., -2.],
                    [9. , 9. , 9.]]
                 )
d_mtrx = b_mtrx @ a_mtrx   # linear algebra matrix-matrix product

print('A =\n', a_mtrx)
print('B =\n', b_mtrx)
print('D =\n', d_mtrx)


# ## NumPy and SciPy Linear Algebra<a id="pylinalg"></a>
# [NumPy](http://www.numpy.org/) has extensive support for [linear algebra](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html?highlight=linear%20algebra) arrays. We collect here the relevant operations for this course.
# However additional resources are instead added to [SciPy](https://docs.scipy.org/doc/scipy-1.1.0/reference/) for general scientific computing including [linear algebra](https://docs.scipy.org/doc/scipy-1.1.0/reference/tutorial/linalg.html).
# 
# Linear algebra operations are obtained from the `linalg` sub-package of the `numpy` package, and the `linalg` sub-package of `scipy`.

# In[ ]:


'''Import the NumPy linear algebra sub-package as usual'''

#import numpy.linalg as linalg
#from numpy import linalg       # often used alternative
'''or leave it commented since the usage of np.linalg is self-documenting'''


# The 2-norm or norm (magnitude) of a vector $\bvec$ is indicated as $\norm{\bvec}$ and computed as follows:

# In[ ]:


'''Vector norm (or magnitude)'''

norm_b_vec = np.linalg.norm( b_vec ) # default norm is the 2-norm 
print('||b|| =', norm_b_vec)      # same as magnitude


# ### Matrix solve <a id="pysolve"></a>
# 
# *Solve*<a id="solve"></a> for $\xvec$ in the matrix equation $\Amtrx\,\xvec = \bvec$, where $\Amtrx = 
# \begin{pmatrix}
# 2 & 1 & 1 \\
# 4 & -6 & 0 \\
# -2 & 7 & 2
# \end{pmatrix}
# $
# and $\bvec = \begin{pmatrix} 5\\ -2\\ 9 \end{pmatrix}$.

# In[ ]:


'''Matrix solver (this is short for solution of a linear algebraic system of equations)'''

x_vec = np.linalg.solve( a_mtrx, b_vec ) # solve linear system for A, b; per course notes 04
print('solution x =', x_vec)


# The residual vector defined as $\rvec = \bvec - \Amtrx\,\xvec$ is of importance. So is its norm $\norm{\rvec}$.

# In[ ]:


'''Verify the accuracy of the solution'''

res_vec = b_vec - a_mtrx @ x_vec
print('b - A x =',res_vec)
print('||b - A x|| =',np.linalg.norm( res_vec ))


# The rank of a matrix of coefficients, $\rank(\Amtrx)$, of a linear algebraic system of equations determines weather the solution is unique or singular.

# In[ ]:


'''Matrix rank'''

k = np.linalg.matrix_rank( a_mtrx )  # rank; per course notes 14
print('rank(A) =',k)
print('shape(A) =',a_mtrx.shape)

if k == a_mtrx.shape[0] and k == a_mtrx.shape[1]:     # flow control
    print('A is non-singular; solution is unique ')
else: 
    print('A is singular')


# Why is this matrix $\Bmtrx$ singular?

# In[ ]:


b_mtrx = np.array( [ [ 2.,  1.,  3.],   # singular
                     [ 4., -6., -2.],
                     [-2.,  7.,  5.]])

k = np.linalg.matrix_rank( b_mtrx )  # rank 
print('rank(B) =',k)
print('shape(B) =',b_mtrx.shape)

if k == b_mtrx.shape[0] and k == b_mtrx.shape[1]:     # flow control
    print('B is non-singular; solution is unique ')
else: 
    print('B is singular')


# In[ ]:


'''Matrix determinant'''

det_a_mtrx = np.linalg.det( a_mtrx ) # determinant; course notes 16
print('det(A) =', det_a_mtrx)


# The inverse matrix is denoted as $\Amtrx^{-1}$ and is computed as the matrix that multiplies $\bvec$ and produces the solution $\xvec$, that is, $\xvec = \Amtrx^{-1}\,\bvec$.

# In[ ]:


'''Matrix inverse'''

a_mtrx_inv = np.linalg.inv( a_mtrx ) # matrix inverse; per course notes 17
print('A^-1 =\n', a_mtrx_inv)


# Recall $\Amtrx^{-1}\,\Amtrx = \Imtrx$ where $\Imtrx$ is the identity matrix.

# In[ ]:


'''Identity matrix'''

i_mtrx = a_mtrx_inv @ a_mtrx  # identity matrix; per course notes 17
print('A^-1 A =\n',i_mtrx)


# Using the inverse, the same solution will be found: $\xvec = \Amtrx^{-1}\,\bvec$.

# In[ ]:


'''Solution using the inverse'''

x_vec_again = a_mtrx_inv @ b_vec  # matrix-vector multiply; per course notes 17
print('solution x =', x_vec_again)


# This is the element-by-element reciprocal of the matrix $(\Amtrx)^{-1}$, which is very different than the inverse.

# In[ ]:


'''Inverse power of a matrix'''

#a_mtrx_to_negative_1 = a_mtrx**(-1) # this will cause an error (division by zero)


# Let's look at the determinant of a larger matrix, say $\Mmtrx$.

# In[ ]:


'''Generate a larger matrix from an image'''

from matplotlib import pyplot as plt     # import the pyplot function of the matplotlib package
plt.rcParams['figure.figsize'] = [20, 4] # extend the figure size on screen output
mtrx = plt.imread('https://raw.githubusercontent.com/dpploy/chen-3170/master/notebooks/images/cermet.png',format='png')
m_mtrx = mtrx[:,:]
plt.figure(1)              # create a figure placeholder
plt.imshow( m_mtrx,cmap='gray')
plt.title('Cermet',fontsize=14)
plt.xlabel('x pixels',fontsize=12)
plt.ylabel('y pixels',fontsize=12)
print('M shape =', m_mtrx.shape)


# In[ ]:


'''Larger matrix determinant'''

det_m_mtrx = np.linalg.det( m_mtrx )    # determinant
print('max(M) =',m_mtrx.max())
print('min(M) =',m_mtrx.min())

print('det(M) = %10.3e (not an insightful number)'%det_m_mtrx) # formatting numeric output

print('rank(M) = ',np.linalg.matrix_rank( m_mtrx, tol=1e-5 ) )


# Let's *solve* for this matrix with $\cvec$ as the right side vector, that is, $\Mmtrx\,\xvec = \cvec$.

# In[ ]:


'''Solve M x = c and plot x'''

c_vec = np.random.random(mtrx.shape[0]) # any c will do it

sol = np.linalg.solve( m_mtrx, c_vec ) # solve linear system for A, b

plt.figure(2)
plt.plot(range(c_vec.size),sol,'k')
plt.title('M x = c',fontsize=20)
plt.xlabel('n',fontsize=18)
plt.ylabel('$c_j$',fontsize=18)
print('')


# In[ ]:


res_vec = c_vec - m_mtrx @ sol
#print('c - M x =',res_vec)
print('||c - M x|| =%12.4e'%np.linalg.norm( res_vec ))


# ### $\Lmtrx$ <a id="pyl"></a> forward solve
# A lower triangular matrix allows for a forward solve.

# In[ ]:


'''L forward solve'''

l_mtrx = np.array( [[1., 0., 0.],   # per course notes 
                    [2., 3., 0.],
                    [4., 5., 6.]] )

b_vec = np.array( [1.,2.,3.] )

x_vec = np.linalg.solve( l_mtrx, b_vec )

np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)


# ### $\Umtrx$ <a id="pyu"></a> backward solve
# An upper triangular matrix allows for a backward solve.

# In[ ]:


'''U backward solve'''

u_mtrx = np.array( [[1., 2., 3.],   # per course notes
                    [0, 4., 5.],
                    [0., 0., 6.]] )

b_vec = np.array( [1.,2.,3.] )

x_vec = np.linalg.solve( u_mtrx, b_vec )

np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)


# ### $\Pmtrx\,\Lmtrx\,\Umtrx$ <a id="pyplu"></a> factorization
# The factors: $\Pmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Pmtrx\,\Lmtrx\,\Umtrx = \Amtrx$ can be obtained from the SciPy linear algebra package. $\Pmtrx$ is a permutation matrix if the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.

# In[ ]:


'''Import only the linear algebra package'''

import scipy.linalg
import numpy as np


# In[ ]:


'''P L U factors of A'''

a_mtrx = np.array( [[1, 2,  3],
                    [4, 5,  6],
                    [7, 8, 10]] )

(p_mtrx, l_mtrx, u_mtrx) = scipy.linalg.lu( a_mtrx )

print('P =\n',p_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('Checking...')
print('PLU - A =\n', p_mtrx @ l_mtrx @ u_mtrx - a_mtrx)


# In[ ]:


'''P^-1 = P^T (i.e. the transpose of a permutation matrix is its inverse)'''

pinv_mtrx = np.linalg.inv(p_mtrx)
print('P^-1 =\n', pinv_mtrx)
print('Checking...')
print('P^-1 - P^T =\n', pinv_mtrx - p_mtrx.transpose())


# In[ ]:


'''PLU x = b; that is: Forward: L y = P^-1 b, Backward: U x = y '''

b_vec = np.array([1.,2.,3.])
y_vec = scipy.linalg.solve(l_mtrx, p_mtrx.transpose() @ b_vec) # L y = P^T b

x_vec = scipy.linalg.solve(u_mtrx, y_vec)                      # U x = y

print('x =', x_vec)

x_vec_gold = scipy.linalg.solve( a_mtrx, b_vec )    # solution using A x = b
print('||x - x_gold|| =',scipy.linalg.norm(x_vec-x_vec_gold))


# In[ ]:


'''Deterninant of U or L: product of the diagonal'''

det_u = np.linalg.det(u_mtrx)
print('det(U) = %8.3e'%det_u)

diag_vec = np.diagonal(u_mtrx)
prod = np.prod(diag_vec)
print('diag(U) product = %8.3e'%prod )


# In[ ]:


'''Determinant of P (always +1 or -1)'''

det_p = np.linalg.det(p_mtrx)
print('det(P) = %8.3e'%det_p)


# In[ ]:


'''Determinant of A = det(PLU)'''

det_l = np.prod( np.diagonal(l_mtrx) )
det_plu = det_p * det_l * det_u # last term is det of L

print('det(PLU) = %8.3e'%det_plu)
print('det(A)   = %8.3e'%np.linalg.det(a_mtrx))


# ## ChEn-3170 Linear Algebra<a id="courselinalg"></a>
# In this course various algorithms need to be programmed. These should be compared to `SciPy` and/or `NumPy`.

# ### $\Lmtrx$ <a id="l"></a> forward solve
# A lower triangular matrix allows for a forward solve.
# The algorithm for $\Lmtrx\,\xvec=\bvec$ is as follows: 
# 
# \begin{equation*}
# x_i = \Bigl(b_i - \sum\limits_{j=1}^{i-1} L_{i,j}\,x_j \Bigr)\,L^{-1}_{i,i} \quad\ \forall \quad\ i=1,\ldots,m
# \end{equation*}
# 
# for $i$ and $j$ with offset 1. Recall that `NumPy` and `Python` have offset 0 for their sequence data types.

# In[ ]:


'''L forward solve'''

l_mtrx = np.array( [[1., 0., 0.],   # per course notes 
                    [2., 3., 0.],
                    [4., 5., 6.]] )

b_vec = np.array( [1.,2.,3.] )

from chen_3170.help import forward_solve 

x_vec = forward_solve( l_mtrx, b_vec )

np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)


# ### $\Umtrx$ <a id="u"></a> backward solve
# A upper triangular matrix allows for a backward solve.
# The algorithm for $\Umtrx\,\xvec=\bvec$ is as follows: 
# 
# \begin{equation*} x_i = \Bigl(b_i - \sum\limits_{j=i+1}^{m} U_{i,j}\,x_j \Bigr)\,U^{-1}_{i,i} \quad\ \forall \quad\ i=m,\ldots,1
# \end{equation*}
# 
# for $i$ and $j$ with offset 1. Recall that `NumPy` and `Python` have offset 0 for their sequence data types.

# In[ ]:


'''U backward solve'''

u_mtrx = np.array( [[1., 2., 3.],   # per course notes
                    [0, 4., 5.],
                    [0., 0., 6.]] )

b_vec = np.array( [1.,2.,3.] )

try:    
    from chen_3170.toolkit import backward_solve   
except ModuleNotFoundError:
    assert False, 'You need to provide your own backward_solve function here. Bailing out.'

x_vec = backward_solve( u_mtrx, b_vec )

np.set_printoptions(precision=3) # one way to control printing of numpy arrays
print('x = ',x_vec)


# ### $\Amtrx = \Lmtrx\,\Umtrx$ <a id="lu"></a> factorization
# $\Lmtrx\,\Umtrx$ factorization algorithm (without using pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ computes the $\Lmtrx\,\Umtrx$ factors. The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that 
# 
# \begin{equation*}
#  A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \quad\ \forall\ i=k+1,\ldots,m \quad\ \text{and} \quad\ j=k,\ldots,m
# \end{equation*}
# 
# where the multipliers $m_{i,k}$ are given by $m_{i,k} = \frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \ \forall \ i>j$,  $L_{i,i}=1$, and $L_{i,j}=0 \ \forall \ i<j$.

# In[ ]:


'''L U factors of A'''

a_mtrx = np.array( [[1., 2.,  3.],
                    [4., 5.,  6.],
                    [7., 8., 10.]] )

try:    
    from chen_3170.toolkit import lu_factorization   
except ModuleNotFoundError:
    assert False, 'You need to provide your own lu_factorization function here. Bailing out.'

(l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx )

print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - A =\n', l_mtrx @ u_mtrx - a_mtrx)


# ### $\Pmtrx\,\Amtrx = \Lmtrx\,\Umtrx$ <a id="plu"></a> factorization
# The factors: $\Pmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Lmtrx\,\Umtrx = \Pmtrx\,\Amtrx$ can be obtained from **partial pivoting** strategy to the $\Lmtrx\,\Umtrx$ factorization algorithm shown above. $\Pmtrx$ is a row permutation matrix obtained by the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.
# 
# Program a $\Lmtrx\,\Umtrx$ factorization algorithm (using partial pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ and compute the $\Pmtrx\,\Lmtrx\,\Umtrx$ factors. 
# 
# The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \ \forall\ i=k+1,\ldots,m \ \text{and}\ j=k,\ldots,m$ where the multipliers $m_{i,k}$ are given by $m_{i,k} = \frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \ \forall \ i>j$,  $L_{i,i}=1$, and $L_{i,j}=0 \ \forall \ i<j$. However, every $k$-step selects a pivot $A^{(k)}_{k,k}$ of maximum absolute value via row
# exchanges recorded in the permutation matrix $\Pmtrx$.

# In[ ]:


'''P L U factors of A'''

a_mtrx = np.array( [[1., 2.,  3.],
                    [4., 5.,  6.],
                    [7., 8., 10.]] )

try:    
    from chen_3170.toolkit import lu_factorization   
except ModuleNotFoundError:
    assert False, 'You need to provide your own lu_factorization function here. Bailing out.'

(p_mtrx, l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx, pivoting_option='partial' )

print('P =\n',p_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - PA =\n', l_mtrx @ u_mtrx - p_mtrx @ a_mtrx)


# ### $\Pmtrx\,\Amtrx\,\Qmtrx = \Lmtrx\,\Umtrx$ <a id="pqlu"></a> factorization
# The factors: $\Pmtrx$, $\Qmtrx$, $\Lmtrx$, and $\Umtrx$ where $\Lmtrx\,\Umtrx = \Pmtrx\,\Amtrx\,\Qmtrx$ can be obtained from a user-developed algorithm with **complete pivoting**. $\Pmtrx$ is a row permutation matrix and $\Qmtrx$ is a column permutation matrix, obtained by the underlying Gaussian elimination used to construct the $\Lmtrx$ and $\Umtrx$ factors.
# 
# Program a $\Lmtrx\,\Umtrx$ factorization algorithm (using complete pivoting) for a square matrix $\overset{(m \times m)}{\Amtrx}$ and compute the $\Pmtrx\,\Qmtrx\,\Lmtrx\,\Umtrx$ factors. The factorization is obtained by elimination steps $k = 1,\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\, m_{i,k} \ \forall\ i=k+1,\ldots,m \ \text{and}\ j=k,\ldots,m$ where the multipliers $m_{i,k}$ are given by $m_{i,k} = \frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \ \forall \ i>j$,  $L_{i,i}=1$, and $L_{i,j}=0 \ \forall \ i<j$. However, every $k$-step selects a pivot $A^{(k)}_{k,k}$ of maximum absolute value via row
# exchanges recorded in the permutation matrix $\Pmtrx$ and column exchanges recorded in the permutation matrix $\Qmtrx$.

# In[ ]:


'''P Q L U factors of A'''

a_mtrx = np.array( [[1., 2.,  3.],
                    [4., 5.,  6.],
                    [7., 8., 10.]] )

try:    
    from chen_3170.toolkit import lu_factorization   
except ModuleNotFoundError:
    assert False, 'You need to provide your own lu_factorization function here. Bailing out.'

(p_mtrx, q_mtrx, l_mtrx, u_mtrx, rank) = lu_factorization( a_mtrx, pivoting_option='complete' )

print('P =\n',p_mtrx)
print('Q =\n',q_mtrx)
print('L =\n',l_mtrx)
print('U =\n',u_mtrx)
print('rank =',rank)
print('')
print('Checking...')
print('LU - PAQ =\n', l_mtrx @ u_mtrx - p_mtrx @ a_mtrx @ q_mtrx)


# In[ ]: