# input a square matrix!
LUdecomposition <- function(m) {

  # throw an error if input is not square
  dimensions <- dim(m)
  if (dimensions[1] != dimensions[2]) {
    return("Matrix not square!")

  } else {
  # calculate lower and upper triangles
    d <- dimensions[1]
    L <- diag(d)
    U <- m

    # iterate over rows 2 through d
    for (x in 1:d) {

      # make sure m[1,1] = 1
      if (U[1,1] != 1) {
        U[1,] <- U[1,] / U[1,1]
      }

      if (x > 1) {
      # continue decomposing the upper triangle
        for (y in 1:(x-1)) {

          # avoid dividing by zero
          nonzero <- 1
          for (z in (x-1):1) {
            if (U[z,y] != 0) {
              nonzero <- z
              break
            }
          }

          alpha = U[x,y] / U[nonzero,y]
          U[x,] <- U[x,] - (alpha * U[nonzero,])

        # add the multiplier to the lower triangle
          L[x,y] <- alpha
        }
      }
    }
  }
    return(list(L, U))
}