# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Funcao para a decomposicao LU:   A = L * U 
#==============================================================================
decomp.LU <- function(A){
	n <- dim(A)[1]
	U <- A
	L <- diag(rep(1,n)) # matriz identidade
	for(i in 1:(n-1)){
		for(j in (i+1):n){
			L[j,i]   <- U[j,i]/U[i,i]
			U[j,i:n] <- U[j,i:n] - L[j,i]*U[i,i:n]
		}
	}
	return(list(L=L,U=U))  
}

#==============================================================================
# Exemplo
#==============================================================================
aux  <- c(1.0,0.7,0.5,0.7,1.0,0.8,0.5,0.8,1.0)
A    <-  matrix(aux,3,3,byrow=T)
print("===============================")
print(A)
print("===============================")
LU.A <- decomp.LU(A)
print(LU.A)
print("===============================")
print(LU.A$L%*%LU.A$U)

[1] "==============================="
     [,1] [,2] [,3]
[1,]  1,0  0,7  0,5
[2,]  0,7  1,0  0,8
[3,]  0,5  0,8  1,0
[1] "==============================="
$L
     [,1]      [,2] [,3]
[1,]  1,0 0,0000000    0
[2,]  0,7 1,0000000    0
[3,]  0,5 0,8823529    1

$U
     [,1] [,2]      [,3]
[1,]    1 0,70 0,5000000
[2,]    0 0,51 0,4500000
[3,]    0 0,00 0,3529412

[1] "==============================="
     [,1] [,2] [,3]
[1,]  1,0  0,7  0,5
[2,]  0,7  1,0  0,8
[3,]  0,5  0,8  1,0

#==============================================================================
# Decomposicao LU utilizando uma funcao do R
#==============================================================================
library("Matrix")
A.lu <- expand(lu(A))
print("===============================")
print(A.lu) # 3 componentes: "L", "U", e "P", a permutacao
print("===============================")
B    <- A.lu$L %*% A.lu$U
print(B)

[1] "==============================="
$L
3 x 3 Matrix of class "dtrMatrix" (unitriangular)
     [,1]      [,2]      [,3]     
[1,] 1,0000000         .         .
[2,] 0,7000000 1,0000000         .
[3,] 0,5000000 0,8823529 1,0000000

$U
3 x 3 Matrix of class "dtrMatrix"
     [,1]      [,2]      [,3]     
[1,] 1,0000000 0,7000000 0,5000000
[2,]         . 0,5100000 0,4500000
[3,]         .         . 0,3529412

$P
3 x 3 sparse Matrix of class "pMatrix"
          
[1,] | . .
[2,] . | .
[3,] . . |

[1] "==============================="
3 x 3 Matrix of class "dgeMatrix"
     [,1] [,2] [,3]
[1,]  1,0  0,7  0,5
[2,]  0,7  1,0  0,8
[3,]  0,5  0,8  1,0

# rm(list=ls()) 
#==============================================================================
# Fim
#==============================================================================