R/LU.R

In rTensor2: MultiLinear Algebra

LU <- function (A)
{
  n <- nrow(A)
  m <- ncol(A)
  decomp <- function(A) {
    n <- nrow(A)
    m <- ncol(A)
    if (n == 1) {
      L <- A
      U <- diag(n)
      result <- list(L = L, U = U)
      return(result)
    }
    else {
      L <- matrix(0, nrow = n, ncol = n)
      U <- L
      a11 <- A[1, 1]
      a12 <- A[1, 2:n]
      a21 <- A[2:n, 1]
      a22 <- A[2:n, 2:n]
      l11 <- 1
      u11 <- a11
      l12 <- matrix(0, nrow = (n - 1), ncol = (n - 1))
      u12 <- a12
      l21 <- a21/u11
      u21 <- matrix(0, nrow = (n - 1), ncol = (n - 1))
      s22 <- a22 - l21 %o% u12
      L[1, 1] <- 1
      L[2:n, 1] <- l21
      U[1, 1] <- u11
      U[1, 2:n] <- u12
      result <- list(L = L, U = U, S = s22)
      return(result)
    }
  }
  if (rankMatrix(A)[[1]] != n)
    stop("matrix is singular")
  if (m != n)
    stop("argument x is not a square matrix")
  L <- matrix(0, nrow = n, ncol = n)
  U <- L
  for (i in 1:n) {
    d <- decomp(A)
    Li <- d$L
    Ui <- d$U
    Si <- d$S
    L[i:n, i] <- Li[, 1]
    U[i, i:n] <- Ui[1, ]
    A <- Si
  }
  return(list(L = L, U = U))
}