# L-16 MCS 275 Wed 20 Feb 2008 : power.py

# Illustration of the power method with numpy.
# The power method to compute the eigenvector
# of an n-by-n matrix A corresponding to the
# dominant eigenvalue builds a sequence of 
# vectors x(i), uses the following formulas:
# y = A*x(i), x(i+1) = y/||y||, i = 0,1,... 
# where 
# ||y|| = sqrt{y[1]^2 + y[2]^2 + .. + y[n]^2}
# and x(0) some vector to start with.
# Use numpy to write a function PowerMethod 
# which takes on input an n-by-n matrix A,
# a vector x(0) of length n, and the number 
# of iterations k.  The function returns x(k).

from numpy import *

def PowerMethod(A,x0,k):
   """
   Returns A^k times x0.
   """
   y = x0
   for i in range(0,k):
      y = dot(A,y)
      s = 0
      for j in range(0,len(y)):
         s = s + y[j]**2
      y = y/sqrt(s)
   return y

def main():
   """
   Interactive test on PowerMethod.
   """
   n = input('Give dimension : ')
   A = random.random((n,n))
   x0 = random.random((n,1))
   k = input('Number of iterations : ')
   xk = PowerMethod(A,x0,k)
   print 'x(%d) =' % k
   print xk
   # check convergence:
   xk1 = PowerMethod(A,x0,k+1)
   print 'x(%d) =' % (k+1)
   print xk1

main()