jordan normal form

J = matrix([[2 ,0, 0, 0, 0, 0],
            [0,-1, 0, 0, 0, 0],
            [0, 0, 3, 0, 0, 0],
            [0, 0, 0, 3, 1, 0],
            [0, 0, 0, 0, 3, 1],
            [0, 0, 0, 0, 0, 3]])

P = matrix([[1, 3, 5, 2,-3,-2],
            [1, 2, 1, 2, 1, 2],
            [1, 2, 3, 4, 5, 6],
            [1, 0, 0, 4, 2, 1],
            [6, 4, 2, 3,-1, 5],
            [0, 1, 2, 0, 1, 2]])

A = (P \ J * P).transpose()
#print(A)

I = identity_matrix(6)

def gauss(M):
    M = copy(M)
    print(M)
    for i in range(M.nrows()):
        r = range(i+1, M.nrows())
        if all(M[j,i] == 0 for j in r):
            break
#        try:
#            k = [M[j,i] != 0 for j in r].index(True)
#        except ValueError:
#            continue
#        if any(abs(M[j,i]) for j in r):
#            j =
        M[i,:] /= M[i,i]
        for j in r:
            M[j,:] -= M[j,i] * M[i,:]
    print(M)
    print(i)
    for i in range(i-1, -1, -1):
        for j in range(i-1, -1, -1):
            M[j,:] -= M[j,i] * M[i,:]
    return M #print(M)

B = 16*(A-3*I)
#gauss(B)

v1 = matrix([[1],[-1],[-3],[4],[2],[0]])
v2 = matrix([[-1],[3],[7],[-4],[0],[4]])
C = block_matrix([[B, 16*v1, 16*v2]])
gC = gauss(C)

w0 = matrix([[31/4],[-1/4],[-33/4],[10],[0],[0]])
D = block_matrix([[B, 16*w0, 16*v1, 16*v2]])
gD = gauss(D)

Q = matrix([[1, 1, 1, 0, 6, 8],
            [2, 3,-1, 1, 5, 0],
            [1, 5,-3, 2, 4,-7],
            [2, 2, 4, 0, 3,11],
            [1,-3, 2, 1, 0, 0],
            [2,-2, 0, 2, 7, 0]])

K = matrix([[-1,0, 0, 0, 0, 0],
            [0, 2, 0, 0, 0, 0],
            [0, 0, 3, 0, 0, 0],
            [0, 0, 0, 3, 1, 0],
            [0, 0, 0, 0, 3, 1],
            [0, 0, 0, 0, 0, 3]])

assert A*Q == Q*K