R/hw2_2.R

In exp500/MAT8054: Provide Solutions of Homework and midterm exam in MAT8054 Class

#' 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)
}