# test_eigen.py   for each lambda   check det|A-lambda*I|=0
#                 check A*v=lambda*v   check max of rows and cols <= lambda
#                 check trace=sum diag|A|= sum(lambda)
#                 check eigen vectors orthogonal

import numpy              # numerical library for Python
from numpy import array   # shorthand, can now just use  array
import numpy.linalg
from numpy.linalg import det
from numpy.linalg import eig
from numpy.matlib import rand,ones,eye,zeros

print "test_eigen.py"

def eig_tests(a):
  e,v = eig(a)
  print "lambda="
  print e
  if len(e) <=4:
    print "vectors="
    print v
    print " "
  # check det|a - lambda*I|
  print "det|a - lambda*I| should be near zero"
  for k in range(4):
    for i in range(4):
      a[i,i] = a[i,i] - e[k]

    z = det(a1)
    print "check eig = ", e[k], ", det = ", z
    for i in range(4):
      a[i,i] = a[i,i] + e[k] # restore a1 for next time through loop

  print " "
  # check eigen vectors orthoganal
  print "sum of abs of dot product of all v should be near zero"
  gsum = 0.0
  for k in range(4-1):
    for k2 in range(k+1,4):
      sum = 0.0
      for j in range(4):
        sum = sum + v[k,j] * v[k2,j]

      gsum = gsum + abs(sum)


  print "orthogonal should be near zero =", gsum
  print " "

  # check a1*v = lambda*v  
  print "a1*v = lambda*v difference should be near zero"
  ve = array([0 for i in range(4)])
  print "len(ve)=", len(ve)
  for k in range(4):
    for i in range(4):
      ve[i] = v[k,i] * e[k]
      va = array([0 for i in range(4)])
      for j in range(4):
        va[i] = va[i] + v[k,j] * a[j,i] # ? order
      
    print "va=" , va
    print "ve= " , ve
    
  print " "
  # check trace = sum diag = sum eigenvalues
  trace = 0.0
  sum = 0.0
  for i in range(4):
    trace +=  a[i,i]
    sum   +=  e[i]

  print "trace = sum ?",trace, sum
  print " "
  # check max sum row or col >= largest eigenvalue
  maxsum = 0.0
  for i in range(4):
    sum = 0.0
    for j in range(4):     # rows
      sum += abs(a[i,j])

    if sum > maxsum:
      maxsum = sum

    sum = 0.0
    for j in range(4):    # cols
      sum += abs(a[j,i])

    if sum > maxsum:
      maxsum = sum # row or col

  print "max sum row col > max abs eigenvalue?", maxsum,  max(abs(e))
  print " "
  print " "
  

a1=array([[1.0/2,1.0/3,1.0/4,1.0/5],[1.0/3,1.0/4,1.0/5,1.0/6],[1.0/4,1.0/5,1.0/6,1.0/7],[1.0/5,1.0/6,1.0/7,1.0/8]])
print "a1="
print  a1

eig_tests(a1)

# check roots of polynomial = eigenvalues

# eigen values and eigen vectors, complex
print "test complex eigenvalues on diagonal"
a2 = zeros((4,4),dtype=complex,order='C')
print "a2="
for i in range(4):
  a2[i,i] = complex(i,i)
print a2

eig_tests(a2)

# using rand
print "generate random matrix, find and test eigenvalues"
a3 = rand(4,4)
print "a3="
print a3

eig_tests(a3)

# end test_eigen.py