import numpy as np

#from copy import copy

def forward$\_$subs$($G$,$b$)$:

rows,cols $=$ G$.$shape

x $=$ np$.$zeros$(($rows$,),$dtype $=$ float$)$

for i in range$($rows$)$:

Gx $=0\mathrm{.}0$

for j in range$($i$)$:

Gx $+=$ G$[$i$,$ j$]*$ x$[$j$]$

x$[$i$]=($b$[$i$]-$ Gx$)/$ G$[$i$,$ i$]$

return x

def back$\_$subs$($G$,$b$)$:

rows,cols $=$ G$.$shape

x $=$ np$.$zeros$(($rows$,),$dtype $=$ float$)$

for i in range$($rows $-1,-1,-1)$:

Gx $=0\mathrm{.}0$

for j in range$($i $+1,$ cols$)$:

Gx $+=$ G$[$i$,$ j$]*$ x$[$j$]$

x$[$i$]=($b$[$i$]-$ Gx$)/$ G$[$i$,$ i$]$

return$($x$)$

def cholesky$\_$factor$($A$)$:

rows, cols $=$ A$.$shape

R $=$ np$.$zeros$(($rows$,$ cols$),$ dtype$=$float$)$

for row in range$($rows$)$:

for col in range$($row $+1)$:

if row $==$ col:

# Diagonal elements

R$[$row$,$ col$]=$ np$.$sqrt$($A$[$row$,$ row$]-$ np$.$sum$($R$[$row$,$ :row$]**2))$

else:

# Off$-$diagonal elements

R$[$row$,$ col$]=($A$[$row$,$ col$]-$ np$.$sum$($R$[$row$,$ :col$]*$ R$[$col$,$ :col$]))/$ R$[$col$,$ col$]$

return R

def solve$\_$system$\_$cholesky$($R$,$b$)$:

# Solve R$^$T y $=$ b using forward substitution

y $=$ forward$\_$subs$($R$.$T$,$ b$)$

# Solve R x $=$ y using backward substitution

x $=$ back$\_$subs$($R$,$ y$)$

return x

def LU$\_$decomp$($A$)$:

rows, cols $=$ A$.$shape

L $=$ np$.$zeros$(($rows$,$ cols$),$ dtype$=$float$)$

U $=$ A$.$copy$()$

for i in range$($rows$)$:

L$[$i$,$ i$]=1$

for j in range$($i$+1,$ rows$)$:

factor $=$ U$[$j$,$ i$]/$ U$[$i$,$ i$]$

L$[$j$,$ i$]=$ factor

U$[$j$,$ i:$]-=$ factor $*$ U$[$i$,$ i:$]$

return$($L$,$U$)$

def solve$\_$system$\_$LU$($L$,$U$,$b$)$:

x $=$ np$.$zeros$($b$.$shape$)$

y $=$ forward$\_$subs$($L$,$ b$)$ # Forward substitution to solve Ly $=$ b

x $=$ back$\_$subs$($U$,$ y$)$ # Backward substitution to solve Ux $=$ y

return$($x$)$

$1)$ Use the forward substitution and back substitution functions above

$2)$ define a function

$$ name : PLU$\_$decomp

$$ inputs : A $(2$D numpy array$)$

$$ output : P $(2$d numpy array$),$ L $(2$d numpy array$),$ U $(2$d numpy array$)$

$$ note $1$ : the function should work for square matrices of any size

$$ note $2$ : use LU decomposition with pivoting

$$ note $3$ : PA $=$ LU

def PLU$\_$decomp$($A$)$:

# add code here

return$($P$,$L$,$U$)$

$3)$ define a function

$$ name : solve$\_$system$\_$PLU

$$ inputs : P $(2$d numpy array$),$ L $(2$D numpy array$),$ U $(2$D numpy array$),$ b $(1$d numpy array$)$

$$ output : x $(1$d numpy array$)$

$$ note $1$ : the function should work for square matrices of any size

$$ note $2$ : the function should solve the equation Ax$=$b where

$$ PA$=$LU is the LU decomposition $($with pivoting$)$ of A

$$ note $3$ : call the forward subs and back subs functions to solve

$4)$ def solve$\_$system$\_$PLU$($P$,$L$,$U$,$b$)$:

x $=$ np$.$zeros$($b$.$shape$)$

# add code here

return$($x$)$