@staticmethod
def PDMA_LU(l2, l1, d, u1, u2):
    """LU decomposition of PDM (for testing only)"""
    n = d.shape[0]
    m = u1.shape[0]
    k = n - m
    for i in range(n - 2 * k):
        lam = l1[i] / d[i]
        d[i + k] -= lam * u1[i]
        u1[i + k] -= lam * u2[i]
        l1[i] = lam
        lam = l2[i] / d[i]
        l1[i + k] -= lam * u1[i]
        d[i + 2 * k] -= lam * u2[i]
        l2[i] = lam
    i = n - 4
    lam = l1[i] / d[i]
    d[i + k] -= lam * u1[i]
    l1[i] = lam
    i = n - 3
    lam = l1[i] / d[i]
    d[i + k] -= lam * u1[i]
    l1[i] = lam
