---
title: "Gaussian Elimination and LU Decomposition"
subtitle: Advanced Statistical Computing
author: Joong-Ho Won
date: today
date-format: "MMMM YYYY"
institute: Seoul National University
execute:
  echo: true  
format:
  revealjs:
    toc: false
    theme: dark
    code-fold: false    
jupyter: julia    
---
versioninfo()

using Pkg
Pkg.activate(pwd())
Pkg.instantiate()
Pkg.status()

A = [2.0 -4.0 2.0; 4.0 -9.0 7.0; 2.0 1.0 3.0]

b = [6.0, 20.0, 14.0]

# Julia way to solve linear equation
# equivalent to `solve(A, b)` in R
A \ b

E21 = [1.0 0.0 0.0; -2.0 1.0 0.0; 0.0 0.0 1.0]

# zero (2, 1) entry
E21 * A   # Step 1, A'

E31 = [1.0 0.0 0.0; 0.0 1.0 0.0; -1.0 0.0 1.0]

# zero (3, 1) entry
E31 * E21 * A  # Step 2, A''

E32 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 5.0 1.0]

# zero (3, 2) entry
E32 * E31 * E21 * A   # Step 3, A'''

M1 = E31 * E21   # inv(L2 * L1)

M2 = E32    # inv(L3)

inv(M1)    # L2 * L1. inv(A) give the inverse of A, but use with caution (see below)

inv(M2)    # L3

A

using LinearAlgebra

# second argument indicates partial pivoting (default) or not
alu = lu(A)
typeof(alu)

fieldnames(typeof(alu))

alu.L

alu.U

alu.p

alu.P

alu.L * alu.U

A[alu.p, :]

# this is doing two triangular solves, 2n^2 flops
alu \ b

det(A) # this does LU!

det(alu) # this is cheap

inv(A) # this does LU!

inv(alu) # this is cheap