Python - Tengo problemas con este codigo para calcular la factorizacion LDL^t, por favor ayuda.

n=4

A=creamatriz(n,n)

B=creamatriz(n,n)

D=creamatriz(n,n)

sumad=0

sumab=0

#B=np.identity(n)

#B=np.diag(n,n)

#b=creamatriz(n,1)

#x=creamatriz(n,1)

A=[[4,3,2,1],[3,3,2,1],[2,2,2,1],[1,1,1,1]]

for j in range(n):

    if j==0:

        for k1 in range (n):

            B[k1][k1]=1

            print("B",k1,B[k1][k1])

    print(np.array(B))

    for k in range(j):

        sumad+=(B[j][k]**2)*D[k][k]

        print("j",j)

        print("B",j+1,k+1,B[j][k])

        print("D",k+1,D[k][k])

    D[j][j]=A[j][j]-sumad

    print("Di",j+1,"-----",D[j][j])

    print(np.array(D))

    for i in range(j+1,n):

        B[j][i]=0

        for k in range(i):

            sumab+=B[i][k]*D[k][k]*B[j][k]

            print("i",i)

            print("Bik",i+1,k+1,B[i][k])

            print("Bjk",j+1,k+1,B[j][k])

        B[i][j]=1/D[j][j]*(A[i][j]-sumab)

        print("Bjk",i+1,j+1)

print(np.array(B))

print(np.array(D))

from matplotlib import pyplot as plt

import math

import numpy as np

import random

from scipy import sparse

from scipy.sparse import csr_matrix

import pandas as pd

##se crea una matriz de tamaño nxm

def creamatriz(n,m):

    matriz = []

    for i in range (n):

        a = [0]*m

        matriz.append(a)

    return matriz

n=4

A=creamatriz(n,n)

B=creamatriz(n,n)

D=creamatriz(n,n)

sumad=0

sumab=0

#B=np.identity(n)

#B=np.diag(n,n)

#b=creamatriz(n,1)

#x=creamatriz(n,1)

A=[[4,3,2,1],[3,3,2,1],[2,2,2,1],[1,1,1,1]]

for j in range(n):

    if j==0:

        for k1 in range (n):

            B[k1][k1]=1

            print("B",k1,B[k1][k1])

    print(np.array(B))

    for k in range(j):

        sumad+=(B[j][k]**2)*D[k][k]

        print("j",j)

        print("B",j+1,k+1,B[j][k])

        print("D",k+1,D[k][k])

    D[j][j]=A[j][j]-sumad

    print("Di",j+1,"-----",D[j][j])

    print(np.array(D))

    for i in range(j+1,n):

        B[j][i]=0

        for k in range(i):

            sumab+=B[i][k]*D[k][k]*B[j][k]

            print("i",i)

            print("Bik",i+1,k+1,B[i][k])

            print("Bjk",j+1,k+1,B[j][k])

        B[i][j]=1/D[j][j]*(A[i][j]-sumab)

        print("Bjk",i+1,j+1)

print(np.array(B))

print(np.array(D))

# Lo que estoy programando es el método de factorización LDL^t que se basa en las formulas para calculo de Lij y los valores de la diagonal D1-D2-D3D-4

# El codigo que intento implementar corresponde a las formulas que estan como imagen, y los resultados que debian ser la solucion, no me estan dando.

#-*- coding: utf -8 -*-

import math

import numpy as np

##se crea una matriz de tamaño nxm

def creamatriz(n,m):

    matriz = [[0]*n for x in range(m)]

    return matriz

n=4

# A=creamatriz(n,n) para que si después la reemplazás

B=creamatriz(n,n)

A=[[4,3,2,1],[3,3,2,1],[2,2,2,1],[1,1,1,1]]

#diagonal positiva en B

for j in range(n):

    for k1 in range (n):

        B[k1][k1]=1

print(np.array(B))