# L-39 MCS 275 Mon 19 Apr 2010 : powermethod.py

# The power method is a simple method for obtaining the
# eigenvector corresponding to the largest eigenvalue.
# We use it to illustrate a computational intensive procedure
# and to see how Cython could help us.

from numpy import *

def MatrixVectorMultiply(A,x):
   """
   Returns the product of A with x.
   """
   n = A.shape[0]
   y = zeros((n,1),double)
   for i in range(0,A.shape[0]):
      for j in range(0,A.shape[1]):
         y[i] = y[i] + A[i,j]*x[j]
   return y

def Norm(x):
   """
   Returns the 2-norm of the vector x.
   """
   z = 0
   for i in range(0,x.shape[0]):
      z = z + x[i]*x[i]
   return sqrt(z)

def NormalizedVector(x):
   """
   Returns x/Norm(x)
   """
   n = Norm(x)
   y = x
   for i in range(0,x.shape[0]):
      y[i] = y[i]/n
   return y

def PowerMethod(A,x,n):
   """
   Given a square numpy matrix in A,
   multiplies x n times with A.
   """
   y = x
   for i in range(0,n):
      y = MatrixVectorMultiply(A,y)
      y = NormalizedVector(y)
   return y

def main():
   """
   Test on the power method.
   """
   n = input('give the dimension : ')
   a = random.normal(0,1,(n,n))
   x = random.normal(0,1,(n,1))
   m = input('give number of iterations : ')
   y = PowerMethod(a,x,m)
   print y
   ans = raw_input('See all eigenvectors ? (y/n) ')
   if ans == 'y':
      [L,V] = linalg.eig(a)
      print V

main()