Creating a function called factorize which takes a square matrix and decomposes into L and U matrices.
factorize <- function(A)
{
if (nrow(A) == ncol(A))
{
size <- nrow(A)
# Constructing a diagonal matrix. This will become the matrix L after factorization.
L <- diag(size)
for (i in 1:(size - 1))
{
for (j in (i + 1):size)
{
L[j, i] <- A[j, i] / A[i, i] # building the L matrix
A[j, ] <- A[j, ] - L[j, i] * A[i, ] # building U matrix
}
}
print(L %*% A)
print('-----L--------')
print(L)
print('--------U--------')
print(A)
}
}
testing with a 3 X 3 Matrix.
A <- matrix(c(3,5,6,4,5,1,5,1,4), nrow = 3)
print(A)
## [,1] [,2] [,3]
## [1,] 3 4 5
## [2,] 5 5 1
## [3,] 6 1 4
factorize(A)
## [,1] [,2] [,3]
## [1,] 3 4 5
## [2,] 5 5 1
## [3,] 6 1 4
## [1] "-----L--------"
## [,1] [,2] [,3]
## [1,] 1.000000 0.0 0
## [2,] 1.666667 1.0 0
## [3,] 2.000000 4.2 1
## [1] "--------U--------"
## [,1] [,2] [,3]
## [1,] 3 4.000000 5.000000
## [2,] 0 -1.666667 -7.333333
## [3,] 0 0.000000 24.800000