sessionInfo()

library(Matrix)

set.seed(123) # seed
n <- 5
A <- Matrix(rnorm(n * n), nrow = n)
b <- rnorm(n)
# a random matrix
A

Al <- lower.tri(A) # Returns a matrix of logicals the same size of a given matrix with entries TRUE in the lower triangle
Aupper <- A
Aupper[Al] <- 0
Aupper

backsolve(Aupper, b) 

A <- t(Matrix(c(2.0, -4.0, 2.0, 4.0, -9.0, 7.0, 2.0, 1.0, 3.0), ncol=3))  # R uses column major ordering
A

b <- c(6.0, 20.0, 14.0)

# R's way tp solve linear equation

solve(A, b)

E21 <- t(Matrix(c(1.0, 0.0, 0.0, -2.0, 1.0, 0.0, 0.0, 0.0, 1.0), nrow=3))  # Step 1, inv(L1)
E21

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

E31 <- t(Matrix(c(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0), nrow=3))  # Step 2, find the corresponding matrix!
E31

# zero (3, 1) entry    
E31 %*% E21 %*% A  # Step2, A''

E32 <- t(Matrix(c(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 5.0, 1.0), nrow=3))  # Step 3, find the corresponding matrix!
E32

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

M1 <- E31 %*% E21  # inv(L1 * L2)
M1

M2 <- E32  # inv(L3)
M2

solve(M1)   # L2; solve(A) gives the inverse of A, but use it with caution (see below)

solve(M2)   # L3

A  # toy example, just to recall

alu <- Matrix::lu(A)
class(alu)

ealu <- expand(alu) # expand the decomposition into Factors

ealu$L

ealu$U

ealu$P

with(ealu, L %*% U)

with(ealu, P %*% L %*% U)

det(A) # this does LU!

det(ealu$U) # this is cheap

solve(A) # this does LU!

with(ealu, solve(U) %*% solve(L) %*% t(P) ) # this is cheap???