Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
LU <- function(A){
# Check if matrix is square.
if(dim(A)[1]!=dim(A)[2]){
return('STOP! Matrix is not square')
}
n <- nrow(A) #number of rows in A
U <- A
L <- diag(n) #assign diagonal to L
for (j in c(1:n)){
for(i in c(2:n)){
if(i > j){
the_row <- U[j,]
multiplier <- U[i, j] / the_row[j]
U[i,] <- U[i,] - (multiplier * the_row)
L[i,j] <- multiplier
}
}
}
return(list(L=L, U=U))
}
# squre matrix A
A <- matrix(c(2,7,3,7,9,4,3,4,7), nrow=3, byrow=TRUE)
A
## [,1] [,2] [,3]
## [1,] 2 7 3
## [2,] 7 9 4
## [3,] 3 4 7
# using LU function to calculate lower and upper triangular matrices.
LU(A)
## $L
## [,1] [,2] [,3]
## [1,] 1.0 0.0000000 0
## [2,] 3.5 1.0000000 0
## [3,] 1.5 0.4193548 1
##
## $U
## [,1] [,2] [,3]
## [1,] 2 7.0 3.000000
## [2,] 0 -15.5 -6.500000
## [3,] 0 0.0 5.225806
#comaring with lu.decomposition function
library(matrixcalc)
A <- matrix(c(2,7,3,7,9,4,3,4,7), nrow=3, byrow=TRUE)
lu.decomposition(A)
## $L
## [,1] [,2] [,3]
## [1,] 1.0 0.0000000 0
## [2,] 3.5 1.0000000 0
## [3,] 1.5 0.4193548 1
##
## $U
## [,1] [,2] [,3]
## [1,] 2 7.0 3.000000
## [2,] 0 -15.5 -6.500000
## [3,] 0 0.0 5.225806