jacobi.py

import numpy as np


EPS = 1e-7


def Givens(i, j, theta, n):
    res = DMatrix([1]*n)
    res[i, j] = np.sin(theta)
    res[j, i] = -np.sin(theta)
    res[i, i] = np.cos(theta)
    res[j, j] = np.cos(theta)
    return res

def DMatrix(diag):
    n = len(diag)
    res = np.zeros([n,n])
    for i in range(n):
        res[i][i] = diag[i]
    return res


def getDLU(mat):
    D = DMatrix(np.diag(mat))
    LU = mat - D
    return (D,LU)


def maxpos(mat):
    n = len(mat)
    p,q = 0,1
    for i in range(n):
        for j in range(n):
            if i == j:continue
            if abs(mat[p, q]) < abs(mat[i, j]):
                p,q = i,j
                
    return (p,q,abs(mat[p, q]))

def view(mat):
    x,y = mat.shape
    for i in range(x):
        tmp = []
        for j in range(y):
            tmp += [mat[i, j]]
            
        print(" ".join(list(map(lambda x:"{0: 10.5f}".format(x), tmp))))





def jacobi(A):
    B = A
    n = len(A)
    U = DMatrix([1] * n)


    itercnt = 0

    while True:
        D, LU = getDLU(B)
        (p,q,val) = maxpos(LU)
        if abs(val) < EPS:
            return (U, np.diag(B))

        theta = np.arctan(-2*B[p, q]/(B[p, p]-B[q, q]))/2
        G = Givens(p, q, theta, n)
        Gt = G.transpose()
        U = U.dot(G)
        B = Gt.dot(B).dot(G)
        itercnt += 1

        print("### iter{0}".format(itercnt))
        print("maxpos:({0}, {1})  maxval:{2}".format(p, q, abs(val)))
        print("theta:{0}".format(theta))
        view(B)
        print()
    
    

def main():
    A = np.matrix([[6, 5, 4, 3, 2, 1],
                   [5, 5, 4, 3, 2, 1],
                   [4, 4, 4, 3, 2, 1],
                   [3, 3, 3, 3, 2, 1],
                   [2, 2, 2, 2, 2, 1],
                   [1, 1, 1, 1, 1, 1]])

    print("###Original Matrix")
    view(A)
    print()
    print()
    
    (U, B) = jacobi(A)
    
    print()
    
    n = len(A)
    D = np.diag(B)

    print("###Eigenvalues/Eigenvectors")
    D = D*np.matrix([[1],[1],[1],[1],[1],[1]])
    eigen = np.concatenate((D, U), 1)

    view(eigen)
        

if __name__ == "__main__":
    main()