#' Solution of Question 2 in HW2
#' @export
#' @param A matrix variable
#' @param b vector variable
#' @param pivot if we use pivoting
hw2_2 <- function(A, b, pivot = TRUE){
if (dim(A)[1] != dim(A)[2]) stop("Matrix is not square.")
n <- ncol(A)
L <- P <- diag(n)
U <- A
for (k in 1:(n-1)){
if (pivot){
piv <- max(abs(U[k:n, k]))
for (j in k:n){
if (abs(U[k, k]) == piv) break
}
U[c(k, j), k:n] <- U[c(j, k), k:n]
L[c(k, j), 1:k-1] <- L[c(j, k), 1:k-1]
P[c(k, j), ] <- P[c(j, k), ]
}
for (j in (k+1):n){
L[j, k] <- U[j, k] / U[k, k]
U[j, k:n] <- U[j, k:n] - L[j, k] * U[k, k:n]
}
}
y <- forwardsolve(L, P%*%b) # Ly = b
x <- backsolve(U, y) # Ux = y
obj=list(L = L, U = U, P = P, x=x)
return(obj)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.