Problem #2 Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
I write a R factorize function for nxn matrix A into LU.
#check input error
checkinput <- function(v,num.dimantion){
total.element =num.dimantion^2
sum=0
for(i in 1: num.dimantion){
sum=sum+abs(v[1+(i-1)*num.dimantion])
}
if(length(v) != total.element ){
print ("It's not a square matrix. Place check your input numbers.")
return (-1)
}else if(sum == 0){
print("The first column is zero vector.")
return (-1)
}else {
#print ("square")
return (0)
}
}
#check pivor numvber of each row and swap next nonzero row to current row while pivor number is zero
checkPivorSwap <- function (M,jcol, num.dimantion) {
if (M[jcol,jcol] == 0) {
temp=M[jcol, ]
for(irow in (jcol+1):num.dimantion){
if(M[irow,jcol] != 0){
M[jcol, ]= M[irow, ]
M[irow, ]= temp
return(M)
}
}
}
return(M)
}
LUfunction<- function (a,num.dimantion){
#check input error
inputError = checkinput(a,num.dimantion)
if ( inputError == 0 ) {
#build original square matrix
A = matrix(a,nrow=num.dimantion,ncol=num.dimantion, byrow=TRUE)
#let lowerMatrix and uperMatrix be an n square matrix to store all factors:L and U
lowerMatrix = diag(num.dimantion)
uperMatrix = A
jcol=1
while(jcol < num.dimantion) {
uperMatrix = checkPivorSwap(uperMatrix,jcol,num.dimantion)
for(irow in (jcol+1):num.dimantion){
#built L matrix
if(uperMatrix[jcol, jcol] != 0) {
factorL= ((uperMatrix[irow, jcol]/uperMatrix[jcol, jcol]))
}else{
factorL=0
}
lowerMatrix[irow, jcol] = factorL
#builr U matrix
uperMatrix[irow,] = ((-1)*(factorL * uperMatrix[jcol, ]) + uperMatrix[irow, ])
}
jcol = jcol +1
}
print("L:lowerMatrix")
print(lowerMatrix)
print("U:uperMatrix")
print(uperMatrix)
print("LU")
print(lowerMatrix%*%uperMatrix)
print("A")
print(A)
}
}
#Sample1
num.dimantion=6
a1 = c(0,1,3,0,0,2,1,0,5,8,5,3,-2,4,8,6,0,8,1,0,1,4,5,11,8,10,-1,0,0,2,-10,0,3,9,0,0) #a1 square matrix A
LUfunction(a1,num.dimantion)## [1] "L:lowerMatrix"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0.0000000 0.00000 0.0000000 0
## [2,] 0 1 0.0000000 0.00000 0.0000000 0
## [3,] -2 4 1.0000000 0.00000 0.0000000 0
## [4,] 1 0 -0.6666667 1.00000 0.0000000 0
## [5,] 8 10 -11.8333333 18.40625 1.0000000 0
## [6,] -10 0 8.8333333 -9.87500 -0.6197183 1
## [1] "U:uperMatrix"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 5 8.00000 5.000000e+00 3.00000
## [2,] 0 1 3 0.00000 0.000000e+00 2.00000
## [3,] 0 0 6 22.00000 1.000000e+01 6.00000
## [4,] 0 0 0 10.66667 6.666667e+00 12.00000
## [5,] 0 0 0 0.00000 -4.437500e+01 -191.87500
## [6,] 0 0 0 0.00000 3.552714e-15 -23.40845
## [1] "LU"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 5 8 5 3
## [2,] 0 1 3 0 0 2
## [3,] -2 4 8 6 0 8
## [4,] 1 0 1 4 5 11
## [5,] 8 10 -1 0 0 2
## [6,] -10 0 3 9 0 0
## [1] "A"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 1 3 0 0 2
## [2,] 1 0 5 8 5 3
## [3,] -2 4 8 6 0 8
## [4,] 1 0 1 4 5 11
## [5,] 8 10 -1 0 0 2
## [6,] -10 0 3 9 0 0#Sample2
num.dimantion=6
a2 = c(0,1,3,0,0,2,1,0,5,8,5,3,-2,4,8,6,0,8,1,0,1,4,5,11,8,10,-1,0,0,2,0,0,0,0,0,0) #a2 square matrix A
LUfunction(a2,num.dimantion)## [1] "L:lowerMatrix"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0.0000000 0.00000 0 0
## [2,] 0 1 0.0000000 0.00000 0 0
## [3,] -2 4 1.0000000 0.00000 0 0
## [4,] 1 0 -0.6666667 1.00000 0 0
## [5,] 8 10 -11.8333333 18.40625 1 0
## [6,] 0 0 0.0000000 0.00000 0 1
## [1] "U:uperMatrix"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 5 8.00000 5.000000 3.000
## [2,] 0 1 3 0.00000 0.000000 2.000
## [3,] 0 0 6 22.00000 10.000000 6.000
## [4,] 0 0 0 10.66667 6.666667 12.000
## [5,] 0 0 0 0.00000 -44.375000 -191.875
## [6,] 0 0 0 0.00000 0.000000 0.000
## [1] "LU"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 5 8 5 3
## [2,] 0 1 3 0 0 2
## [3,] -2 4 8 6 0 8
## [4,] 1 0 1 4 5 11
## [5,] 8 10 -1 0 0 2
## [6,] 0 0 0 0 0 0
## [1] "A"
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 1 3 0 0 2
## [2,] 1 0 5 8 5 3
## [3,] -2 4 8 6 0 8
## [4,] 1 0 1 4 5 11
## [5,] 8 10 -1 0 0 2
## [6,] 0 0 0 0 0 0#Sample3
num.dimantion=2
a3 = c(2,3,4)
LUfunction(a3,num.dimantion)## [1] "It's not a square matrix. Place check your input numbers."#Sample4
num.dimantion=3
a4 = c(0,3,4,0,0,1,0,8,-1)
LUfunction(a4,num.dimantion)## [1] "The first column is zero vector."