Raw File
lu.R
##
##  l u. R  LU Decomposition
##


lu <- function(A, scheme = c("kji", "jki", "ijk")) {
    stopifnot(is.numeric(A), is.matrix(A))
    n <- nrow(A)
    if (ncol(A) != n || n <= 1)
        stop("Argument 'A' must be a square, positive definite matrix.")

    scheme <- match.arg(scheme)

    if (scheme == "kji") {
        for (k in 1:(n-1)) {
            if (A[k, k] == 0)
                stop("All diagonal elements of matrix 'A' must be non-zero.")
            for (i in (k+1):n) {
                A[i, k] <- A[i, k] / A[k, k]
                A[i, (k+1):n] <- A[i, (k+1):n] - A[i, k] * A[k, (k+1):n]
            }
        }

    } else if (scheme == "jki") {
        if (A[1, 1] == 0)
            stop("All diagonal elements of matrix 'A' must be non-zero.")
        i <- 2:n
        A[i, 1] <- A[i, 1] / A[1, 1]
        for (j in 2:n) {
            if (A[j, j] == 0)
                stop("All diagonal elements of matrix 'A' must be non-zero.")
            for (k in 1:(j-1)) {
                i <- (k+1):n
                A[i, j] <- A[i, j] - A[i, k] * A[k, j]
            }
            if (j < n) {
                i <- (j+1):n
                A[i, j] <- A[i, j] / A[j, j]
            }
        }

    } else if (scheme == "ijk") {       # 'compact' Doolittle scheme
        for (i in 2:n) {                # j in 1:n
            for (j in 2:i) {
                if (A[j, j] == 0)
                    stop("All diagonal elements of matrix 'A' must be non-zero.")
                A[i, j-1] <- A[i, j-1] / A[j-1, j-1]
                k <- 1:(j-1)
                A[i, j] <- A[i, j] - A[i, k] %*% A[k, j]
            }
            if (i < n) {
                k <- 1:(i-1)
                for (j in (i+1):n) {
                    A[i, j] <- A[i, j] - A[i, k] %*% A[k, j]
                }
            }
        }
    }

    L <- eye(n) + tril(A, -1)
    U <- triu(A)
    return(list(L = L, U = U))
}

back to top