LDLT分解Python

#!/usr/bin/env python3
'''
LDLT
'''
import numpy as np
 
def LDLT(amatrix):

    if len(np.shape(amatrix)) != 2 or np.shape(amatrix)[0]!=np.shape(amatrix)[1]:
        print("error shape")
        return
    
    for i in range(np.shape(amatrix)[0]):
        for j in range(np.shape(amatrix)[1]):
            if amatrix[i][j] != amatrix[j][i]:
                print("The input matrix should be symmetric")
                return 
        
    n = np.shape(amatrix)[0] #dimension of matrix
 
    l = np.eye(n)
 
    d = np.zeros((n,n))
 
    for k in range(n):
 
        if k == 0:
 
            d[k][k] = amatrix[k][k]

            if d[k][k] == 0:
    
                print('error matrix type with 0 sequential principal minor determinant')
    
                return
            
            for m in range(n):
 
                l[m][k] = amatrix[m][k]/d[k][k]
        else:
 
            temp_sum1 = 0
 
            for m in range(k):
 
                temp_sum1 += amatrix[k][m]*l[k][m]
 
            d[k][k] = amatrix[k][k] - temp_sum1
 
            for j in range(k+1,n):
 
                temp_sum2 = 0
 
                for m in range(k):
 
                    temp_sum2 += amatrix[j][m]*l[k][m]
 
 
                amatrix[j][k] = amatrix[j][k] - temp_sum2
 
 
                l[j][k] = amatrix[j][k]/d[k][k]
 
    print("l")
    print(l)
    print('d')
    print(d)
 
    return l,d
 
 
 
a=[[5,-4,1],[-4,6,-4],[1,-4,-6]]

l,d = LDLT(a)
 
print(np.dot(np.dot(l,d),l.T))
#print(np.shape(a))