python code for solving eigenvalue problem by Jacobi’s method，eigenvaluejacobi,''' lam,x =

python code for solving eigenvalue problem by Jacobi’s method，eigenvaluejacobi,''' lam,x =

''' lam,x = jacobi(a,tol = 1.0e-9).    Solution of std. eigenvalue problem [a]{x} = lam{x}    by Jacobi's method. Returns eigenvalues in vector {lam}    and the eigenvectors as columns of matrix [x].'''from numpy import array,identity,diagonalfrom math import sqrtdef jacobi(a,tol = 1.0e-9): # Jacobi method    def maxElem(a): # Find largest off-diag. element a[k,l]        n = len(a)        aMax = 0.0        for i in range(n-1):            for j in range(i+1,n):                if abs(a[i,j]) >= aMax:                    aMax = abs(a[i,j])                    k = i; l = j        return aMax,k,l    def rotate(a,p,k,l): # Rotate to make a[k,l] = 0        n = len(a)        aDiff = a[l,l] - a[k,k]        if abs(a[k,l]) < abs(aDiff)*1.0e-36: t = a[k,l]/aDiff        else:            phi = aDiff/(2.0*a[k,l])            t = 1.0/(abs(phi) + sqrt(phi**2 + 1.0))            if phi < 0.0: t = -t        c = 1.0/sqrt(t**2 + 1.0); s = t*c        tau = s/(1.0 + c)        temp = a[k,l]        a[k,l] = 0.0        a[k,k] = a[k,k] - t*temp        a[l,l] = a[l,l] + t*temp        for i in range(k):      # Case of i < k            temp = a[i,k]            a[i,k] = temp - s*(a[i,l] + tau*temp)            a[i,l] = a[i,l] + s*(temp - tau*a[i,l])        for i in range(k+1,l):  # Case of k < i < l            temp = a[k,i]            a[k,i] = temp - s*(a[i,l] + tau*a[k,i])            a[i,l] = a[i,l] + s*(temp - tau*a[i,l])        for i in range(l+1,n):  # Case of i > l            temp = a[k,i]            a[k,i] = temp - s*(a[l,i] + tau*temp)            a[l,i] = a[l,i] + s*(temp - tau*a[l,i])        for i in range(n):      # Update transformation matrix            temp = p[i,k]            p[i,k] = temp - s*(p[i,l] + tau*p[i,k])            p[i,l] = p[i,l] + s*(temp - tau*p[i,l])    n = len(a)    maxRot = 5*(n**2)       # Set limit on number of rotations    p = identity(n)*1.0     # Initialize transformation matrix    for i in range(maxRot): # Jacobi rotation loop         aMax,k,l = maxElem(a)        if aMax < tol: return diagonal(a),p        rotate(a,p,k,l)    print 'Jacobi method did not converge'