LU-Decomposition

# install.packages("pracma") # Only required if you haven't installed it before
# library(pracma)

A <- matrix(c(2, 1, -2, 3, 2, 2, 5, 4, 3), nrow = 3, byrow = TRUE)
b <- matrix(c(1, 7, 12), ncol = 1)
A0 <- A
b0 <- b
n <- nrow(A)
I <- diag(n)
k <- 0  # Initialize phase counter

k <- k + 1  # Reset/update phase counter
# Locate position of element largest modulus
# in the pivot element
r <- which.max(abs(A[k:n, k]))
# Create pivot matrix to interchange
# rows r & k
P1 <- I
P1[c(k, r), ] <- P1[c(r, k), ]
# Interchange rows r & k
A <- P1 %*% A
# Create the Gaussian Transformation matrix
i <- (k + 1):n
M1 <- I
M1[i, k] <- -A[i, k] / A[k, k]
# Zero out elements below the pivot element
A <- M1 %*% A

print(M1)
print(A)

k <- k + 1
# Locate position of element largest modulus in the pivot element
r <- which.max(abs(A[k:n, k]))
r <- r + k - 1 # Correct position in pivot column
# Create pivot matrix to interchange rows r & k
P2 <- I
P2[c(k, r), ] <- P2[c(r, k), ]
# Interchange rows r & k of A
A <- P2 %*% A
# Interchange rows/columns r & k of M1
M1 <- P2 %*% M1 %*% P2

print(M1)
# Create the Gaussian Transformation matrix
i <- (k + 1):n
M2 <- I
M2[i, k] <- -A[i, k] / A[k, k]
# Zero out elements below the pivot element
A <- M2 %*% A
print(M2)
print(A)



A0 <- A
b0 <- b
n <- nrow(A)
I <- diag(n)
k <- 0  # Initialize phase counter

# Calculate P, L, and U matrices
k <- k + 1
r <- which.max(abs(A[k:n, k]))
r <- r + k - 1
P1 <- I
P1[c(k, r), ] <- P1[c(r, k), ]
A <- P1 %*% A
i <- (k + 1):n
M1 <- I
M1[i, k] <- -A[i, k] / A[k, k]
A <- M1 %*% A
P2 <- I
P2[c(k, r), ] <- P2[c(r, k), ]
P <- P2 %*% P1
L <- solve(M2 %*% M1)
U <- A

# LU decomposition using pracma package
LU <- lu(A0)
Lm <- LU$L
Um <- LU$U
Pm <- LU$P

# Solve L*y = P*b
y <- solve(L, P %*% b)

# Solve U*x = y
x <- solve(U, y)

# Solve A0\b0
xm <- solve(A0, b0)

# Output the results
print("Permutation matrix P:")
print(P)
print("LU decomposition:")
print("L:")
print(Lm)
print("U:")
print(Um)

print("Solving L*y = P*b:")
print(y)
print("Solving U*x = y:")
print(x)

print("Solving A0\\b0:")
print(xm)


# First part of the code

k <- k + 1
r <- which.max(abs(A[k:n, k]))
r <- r + k - 1
P3 <- I
P3[c(k, r), ] <- P3[c(r, k), ]
A <- P3 %*% A
M1 <- P3 %*% M1 %*% P3
i <- (k + 1):n
M3 <- I
M3[i, k] <- -A[i, k] / A[k, k]
A <- M3 %*% A

# Calculate permutation matrix P, lower triangular matrix L, and upper triangular matrix U
P <- P3 %*% P2 %*% P1
L <- solve(M3 %*% M2 %*% M1)
U <- A

# Display results
print("Permutation matrix P:")
print(P)
print("Lower triangular matrix L:")
print(L)
print("Upper triangular matrix U:")
print(U)


# Second part of the code

A <- matrix(c(2, 4, 3, 2, 3, 6, 5, 2, 2, 5, 2, -3, 4, 5, 14, 14), nrow = 4, byrow = TRUE)
n <- nrow(A)
I <- diag(n)
k <- 0  # Initialise phase counter
L <- matrix(0, nrow = n, ncol = n)  # Initialise Lower triangular matrix

k <- k + 1
r <- which.max(abs(A[k:n, k]))
r <- r + k - 1
P1 <- I
P1[c(k, r), ] <- P1[c(r, k), ]
A <- P1 %*% A
i <- (k + 1):n
M1 <- I
M1[i, k] <- -A[i, k] / A[k, k]
L[i, k] <- A[i, k] / A[k, k]
A <- M1 %*% A


# Your previous code

k <- k + 1
r <- which.max(abs(A[k:n, k]))
r <- r + k - 1
P2 <- I
P2[c(k, r), ] <- P2[c(r, k), ]
A <- P2 %*% A
i <- (k + 1):n
M2 <- I
M2[i, k] <- -A[i, k] / A[k, k]
L <- P2 %*% L
L[i, k] <- A[i, k] / A[k, k]
A <- M2 %*% A

k <- k + 1
r <- which.max(abs(A[k:n, k]))
r <- r + k - 1
P3 <- I
P3[c(k, r), ] <- P3[c(r, k), ]
A <- P3 %*% A
i <- (k + 1):n
M3 <- I
M3[i, k] <- -A[i, k] / A[k, k]
L <- P3 %*% L
L[i, k] <- A[i, k] / A[k, k]
A <- M3 %*% A

L <- I + L
U <- A
P <- P3 %*% P2 %*% P1

# Display results
print("Permutation matrix P:")
print(P)
print("Lower triangular matrix L:")
print(L)
print("Upper triangular matrix U:")
print(U)


# Second part of the code

A <- matrix(c(2, 4, 3, 2, 3, 6, 5, 2, 2, 5, 2, -3, 4, 5, 14, 14), nrow = 4, byrow = TRUE)
U <- A
n <- nrow(A)
L <- matrix(0, nrow = n, ncol = n)
I <- diag(n)
P <- I

for (k in 1:(n - 1)) {
  r <- which.max(abs(U[k:n, k]))
  r <- r + k - 1
  # Permutation matrix to swap rows r & k
  Pk <- I
  Pk[c(k, r), ] <- Pk[c(r, k), ]
  # Swap rows r & k
  U <- Pk %*% U
  # Create kth Gaussian Transformation matrix
  i <- (k + 1):n
  m <- U[i, k] / U[k, k]
  Mk <- I
  Mk[i, k] <- -m
  # Update elements of L
  L <- Pk %*% L
  L[i, k] <- m
  # Update permutation matrix
  P <- Pk %*% P
  # Zero out elements below pivot element
  U <- Mk %*% U
}
L <- I + L

# Display results
print("Permutation matrix P:")
print(P)
print("Lower triangular matrix L:")
print(L)
print("Upper triangular matrix U:")
print(U)



# LU decomposition of A0
LU <- lu(A0)
Lm <- as.matrix(LU$L)
Um <- as.matrix(LU$U)
Pm <- as.matrix(LU$P)

# Display results
print("Permutation matrix P:")
print(P)
print("Permutation matrix from LU decomposition, Pm:")
print(Pm)
print("Lower triangular matrix L:")
print(L)
print("Lower triangular matrix from LU decomposition, Lm:")
print(Lm)
print("Upper triangular matrix U:")
print(U)
print("Upper triangular matrix from LU decomposition, Um:")
print(Um)

# Solve L*y = P*b
y <- solve(L, P %*% b)

# Solve U*x = y
x <- solve(U, y)

# Solve A0\b0
xm <- solve(A0, b0)

# Output the solutions
print("Solution:")
print(xm)


# Your second part of the code

A <- matrix(c(pi, -exp(1), sqrt(2), -sqrt(3), pi^2, exp(1), -exp(1)^2, 3/7, sqrt(5), -sqrt(6), 1, -sqrt(2), pi^3, exp(1)^2, -sqrt(7), 1/9), nrow = 4, byrow = TRUE)
U <- A
n <- nrow(A)
L <- matrix(0, nrow = n, ncol = n)
I <- diag(n)
P <- I

for (k in 1:(n - 1)) {
  i <- (k + 1):n
  m <- U[i, k] / U[k, k]
  Mk <- I
  Mk[i, k] <- -m
  L[i, k] <- m
  U <- Mk %*% U
}
L <- I + L

# Display results
print("Permution matrix, P:")
print(P)
print("Lower triangular matrix, L:")
print(L)
print("Upper triangular matrix, U:")
print(U)

# Solve the system using back substitution
solution <- solve(U, solve(L, P %*% b))

# Output the solutions
print("Solution:")
print(solution)