Inverse power method applied to a tridiagonal matrix in Python
''' lam,x = inversePower3(d,c,s,tol=1.0e-6).
Inverse power method applied to a tridiagonal matrix
[A] = [c\d\c]. Returns the eigenvalue closest to 's'
and the corresponding eigenvector.
'''
from numpy import dot,zeros
from LUdecomp3 import *
from math import sqrt
from random import random
def inversePower3(d,c,s,tol=1.0e-6):
n = len(d)
e = c.copy()
cc = c.copy() # Save original [c]
dStar = d - s # Form [A*] = [A] - s[I]
LUdecomp3(cc,dStar,e) # Decompose [A*]
x = zeros(n)
for i in range(n): # Seed [x] with random numbers
x[i] = random()
xMag = sqrt(dot(x,x)) # Normalize [x]
x =x/xMag
flag = 0
for i in range(30): # Begin iterations
xOld = x.copy() # Save current [x]
LUsolve3(cc,dStar,e,x) # Solve [A*][x] = [xOld]
xMag = sqrt(dot(x,x)) # Normalize [x]
x = x/xMag
if dot(xOld,x) < 0.0: # Detect change in sign of [x]
sign = -1.0
x = -x
else: sign = 1.0
if sqrt(dot(xOld - x,xOld - x)) < tol:
return s + sign/xMag,x
print 'Inverse power method did not converge'
