Hello I try to implement from basis (I need it for GPU kernel (just small 3 by 3 matrix one of multiple things that a kernel do - no point in using cublas …) so -can not using Linear algebra package for a moment , also for some unclear reason the julia linear algebra implementation is using inversion which purpose is hidden from me- and from my perspective much complicating thing)
Hence I try to adapt python code from
Cholesky Decomposition : Matrix Decomposition - GeeksforGeeks
and I know that it is very simple task, yet for some reason I am stuck for hours and I can not make it work
python working code
# Python3 program to decompose
# a matrix using Cholesky
# Decomposition
import math
MAX = 100;
def Cholesky_Decomposition(matrix, n):
lower = [[0 for x in range(n + 1)]
for y in range(n + 1)];
# Decomposing a matrix
# into Lower Triangular
for i in range(n):
for j in range(i + 1):
sum1 = 0;
# summation for diagonals
if (j == i):
for k in range(j):
sum1 += pow(lower[j][k], 2);
lower[j][j] = int(math.sqrt(matrix[j][j] - sum1));
else:
# Evaluating L(i, j)
# using L(j, j)
for k in range(j):
sum1 += (lower[i][k] *lower[j][k]);
if(lower[j][j] > 0):
lower[i][j] = int((matrix[i][j] - sum1) /
lower[j][j]);
# Displaying Lower Triangular
# and its Transpose
print("Lower Triangular\t\tTranspose");
for i in range(n):
# Lower Triangular
for j in range(n):
print(lower[i][j], end = "\t");
print("", end = "\t");
# Transpose of
# Lower Triangular
for j in range(n):
print(lower[j][i], end = "\t");
print("");
# Driver Code
n = 3;
matrix = [[4, 12, -16],
[12, 37, -43],
[-16, -43, 98]];
Cholesky_Decomposition(matrix, n);
# This code is contributed by mits
my Julia implementation that do not give correct results
using LinearAlgebra
A = [ 6 15 55;15 55 225;55 225 979 ]
L = zeros(3,3)
L = zeros(3,3)
lower =L
matrix = A
n=3
for i in 0:(n-1)
for j in 0:i
sum1 = 0;
# summation for diagonals
if (j == i)
for k in 0:j-1
sum1 += lower[j+1,k+1]*lower[j+1,k+1]
lower[j+1,j+1] = sqrt(matrix[j+1,j+1] - sum1)
end#for
else
# Evaluating L(i, j)
# using L(j, j)
for k in 0:j-1
sum1 += (lower[i+1,k+1] *lower[j+1,k+1]);
if(lower[j+1,j+1] > 0)
lower[i+1,j+1] = (matrix[i+1,j+1] - sum1) /lower[j+1,j+1]
end#if
end#for
end#if else
end#for
end#for
#this gives diffrent results
L
cholesky([ 6 15 55;15 55 225;55 225 979 ])
Where is my mistake?
