# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Funcao para a solucao de sistemas lineares dado um sistema triangular 
# inferior
#==============================================================================
sistema.triang <- function(L,b){
    n    <- length(b)
    x    <- rep(NA,n)
    x[1] <- b[1]/L[1,1]
    for(i in 2:n)
        x[i] <- (b[i]-sum(L[i,1:(i-1)]*x[1:(i-1)]))/L[i,i]
    return(x)
}

#==============================================================================
# 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

#==============================================================================
# Utilizando as funcoes dos exemplos anteriores
#==============================================================================
L     <- t(decomp.cholesky(A)) # Note a transposicao
L.inv <- matrix(0,3,3) 
I     <- diag(c(1,1,1))        # Matriz identidade
for(i in 1:3){
	L.inv[,i] <- sistema.triang(L,I[,i])
} 
print("===============================")
print(L.inv)
print("===============================")
print(t(L.inv)%*%L.inv)

[1] "==============================="
           [,1]      [,2]     [,3]
[1,]  1,0000000  0,000000 0,000000
[2,] -0,9801961  1,400280 0,000000
[3,]  0,1980295 -1,485221 1,683251
[1] "==============================="
           [,1]      [,2]       [,3]
[1,]  2,0000000 -1,666667  0,3333333
[2,] -1,6666667  4,166667 -2,5000000
[3,]  0,3333333 -2,500000  2,8333333

#==============================================================================
# Utilizando a funcao do R
#==============================================================================
print("===============================")
print(solve(A)) # Diretamente
print("===============================")
print(solve(chol(A))%*%solve(t(chol(A))))

[1] "==============================="
           [,1]      [,2]       [,3]
[1,]  2,0000000 -1,666667  0,3333333
[2,] -1,6666667  4,166667 -2,5000000
[3,]  0,3333333 -2,500000  2,8333333
[1] "==============================="
           [,1]      [,2]       [,3]
[1,]  2,0000000 -1,666667  0,3333333
[2,] -1,6666667  4,166667 -2,5000000
[3,]  0,3333333 -2,500000  2,8333333

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