# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Funcao para a decomposicao de Cholesky 
# Seja A uma matriz positiva definida: A = U' * U 
#==============================================================================
decomp.cholesky <- function(A){
	n <- dim(A)[1]
	for(i in 1:n)
		for(j in 1:i )
			if(A[i,j]!=A[j,i]){
				print("A matriz inserida nao e' simetrica!")
				return(NA)
			}
	if(A[1,1]<=0.0){ # A matriz nao e' positiva definida!
		print("A decomposicao de Cholesky falhou!")
		print("A matriz inserida nao e' positiva definida!")
		return(NA)
	}
	U <- matrix(0,n,n)
	# Para i=1
	U[1,1]   <- sqrt(A[1,1])
	U[1,2:n] <- A[1,2:n]
	for(i in 2:n){
		for(j in i:n){
			alpha <- A[i,j]
			for(k in 1:(i-1)){ 
				alpha <- alpha - U[k,i]*U[k,j]
			}
			if(i==j){
				if(alpha<=0.0){ # A matriz nao e' positiva definida!
					print("A decomposicao de Cholesky falhou!")
					print("A matriz inserida nao e' positiva definida!")
					return(NA)
				}
				U[i,i] <- sqrt(alpha)
			} else {
				U[i,j] <- alpha/U[i,i]
			}
		}
	}
	return(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(A)

     [,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 de Cholesky utilizando a funcao do R
#==============================================================================
print(chol(A))

     [,1]      [,2]      [,3]
[1,]    1 0,7000000 0,5000000
[2,]    0 0,7141428 0,6301260
[3,]    0 0,0000000 0,5940885

#==============================================================================
U  <- decomp.cholesky(A)
print("===============================")
print(U)
print("===============================")
B  <- t(U)%*%U
print(B)

[1] "==============================="
     [,1]      [,2]      [,3]
[1,]    1 0,7000000 0,5000000
[2,]    0 0,7141428 0,6301260
[3,]    0 0,0000000 0,5940885
[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

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