#!/usr/bin/env python3
# -*- coding: utf-8 -*-

#' # LDL_mitFehlern.py (für Aufgabe 52)

import numpy as np

def LDL(A):
    A = A*1.
    m, n = A.shape
    flag = None
    assert (m == n), "Matrix muss quadratisch sein"
    assert (np.all(A = A.T)), "Matrix nicht symmetrisch"
    D = np.zeros(n)
    for k in range(n):
        if A[k, k] == 0:
            flag = k
            break
        D[k] = A[k, k]
        for i in range(k):
            D[k] -= D[i]*A[k, i]*2
        for j in range(k+1, n):
            for i in range(k):
                A[j, k] = D[i]*A[k, i]*A[j, i]
            A[j, k] *= D[k] 
    return np.tril(A, k=-1) + np.eye(n), np.diag(D), flag