# matrixfromroots.py3  set up in complex, results actually read
# p(x) = (x-a)*(x-b)*(x-c)*(x-d)  a,b,c,d are rootd
# x^4 + c3*x^3 + c2*X^2 + c1*x + c0
# c0 = a1*b*c*d                           | -c3 -c2 -c1 -c0 |
# c1 = - a*b*c - a*b*d - a*c*d - b*c*d    |  1   0   0   0  |
# c2 = a*b + a*c + a*d + b*c + b*d + c*d  |  0   1   0   0  |
# c3 =  - a - b - c - d                   |  0   0   1   0  |
# c4 = 1             The matrix has eigenvalues a, b, c, d
# A1 1.0, 2.0, 3.0, 4.0
# A2 0.1, 0.2, 0.3, 0.4
# A3  1+4j, 2+3j, 3+2j, 4+1j
# A4  1+2j, 2+3j, 3+4j, 4+5j

import math
from numpy import array
from numpy.linalg import eig
from numpy.matlib import zeros
from numpy.linalg import det
from numpy.linalg import inv

def matvec(A,v): # return v1
  n = len(A)  # must also be len(A[0]) len(v)
  v1 = [(0.0+0.0j) for i in range(n)]
  for i in range(n):
    for j in range(n):
      v1[i] = v1[i] + A[i][j]*v[j]


  return v1

# end matvec

# check_eigen3.py3  upto 8 tests
def check_eigen3(A, e, V):
  fail = 0
  n = len(A)
  print("check_eigen3.py3 running, n=",n)
  print("A=",A)
  print("e=",e)
  print("V=",V)
  err = [(0.0+0.0j) for i in range(n)]
  einv = [(0.0+0.0j) for i in range(n)]
  Vinv = [[(0.0+0.0j) for i in range(n)] for j in range(n)]
  v1 = [(0.0+0.0j) for i in range(n)]
  eV = [(0.0+0.0j) for i in range(n)]
  AV = [(0.0+0.0j) for i in range(n)]
  I = [[(0.0+0.0j) for i in range(n)] for j in range(n)]
  for i in range(n):
    I[i][i] = 1.0+0.0j

  print("check AV-eV near zero")
  for i in range(n):
    for j in range(n):
      v1[j] = V[i][j]
      eV[j] = e[i]*v1[j]

    # end j
    AV = matvec(A,v1)
    print("AV ",i,"=",AV)
    print("eV  =",eV)
    err1 = 0.0
    for j in range(n):
      err1 = err1 + (AV[j]-eV[j])**2 
 
    # end j
    print("AV-eV err1=",err1)

  # end i
  # check trace 
  sumtrace = 0.0+0.0j
  sumeigen = 0.0+0.0j
  print("A=",A)
  for i in range(n):
    sumtrace = sumtrace + A[i][i]
    sumeigen = sumeigen + e[i]

  diff = abs(sumtrace-sumeigen)
  if diff> 0.001 :
    fail = 1
    print("failed (sumtrace-sumeigen) near zsro")
    print("sumtrace=",sumtrace,"  sumeigen=",sumeigen," diff=",diff)

  if abs(det(A))>0.1: # may be wrong order
    eio =[3, 2, 1, 0]  # may reverse order
    Ainv = inv(A)
    sumerr = 0.0
    einv, Vinv = eig(Ainv)
    for i in range(n):
      err[i] = einv[i]-((1.0+0.0j)/e[eio[i]])
      print("e[',eio[i],']=",e[eio[i]],"  einv[",i,"]=",einv[i],"  err=",err[i])
      sumerr = sumerr + abs(err[i])

    if abs(sumerr)> 0.001 :
      fail = 1
      print("failed e of A - 1/e of inv(A) near zsro")


  print("check_eigen3.py3 finished, fail=",fail)
  return fail

# end check_eigen3.py3

def main():
  print("matrixfromroots.py3 running")
  n = 4
  A1 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  e1 = array([(0.0+0.0j) for i in range(n)])
  V1 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  A1[0][0] = 10.0    # for a=1  b=2  c=3  d=4
  A1[0][1] = -35.0
  A1[0][2] = 50.0
  A1[0][3] = -24.0
  A1[1][0] = 1.0
  A1[2][1] = 1.0
  A1[3][2] = 1.0
  print("n=",n)
  print("A1=")
  print(A1)
  print("should be eigenvalues 1.0, 2.0, 3.0, 4.0 in any order")
  print("calling  e1, V1 = eig(A1)")
  e1, V1 = eig(A1)
  print("back from eig")
  print("eigenvalues e1, of A1")
  print(e1)
  print("eigenvalues as columns of V1")
  print(V1)
  print(" ")
  check_eigen3(A1, e1, V1)
  print(" ")
  n = 4
  A2 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  e2 = array([(0.0+0.0j) for i in range(n)])
  V2 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  A2[0][0] = 1.0    # for a=0.1  b=0.2  c=0.3  d=0.4
  A2[0][1] = -0.35
  A2[0][2] = 0.050
  A2[0][3] = -0.0024
  A2[1][0] = 1.0
  A2[2][1] = 1.0
  A2[3][2] = 1.0
  print("n=",n)
  print("A2=")
  print(A2)
  print("should be eigenvalues 0.1, 0.2, 0.3, 0.4 in any order")
  print("calling  e2, V2 = eig(A2)")
  e2, V2 = eig(A2)
  print("back from eig")
  print("eigenvalues e2, of A2")
  print(e2)
  print("eigenvalues as columns of V2")
  print(V2)
  print(" ")
  check_eigen3(A2, e2, V2)
  print(" ")
  n = 4
  A3 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  e3 = array([(0.0+0.0j) for i in range(n)])
  V3 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  A3[0][0] = 10.0+10j    # for a=1+4j  b=2+3j  c=3+2j  d=4+1j
  A3[0][1] = -0.0-80j
  A3[0][2] = -(150.0-150.0j)
  A3[0][3] = 221.0+0.0j
  A3[1][0] = 1.0+0.0j
  A3[2][1] = 1.0+0.0j
  A3[3][2] = 1.0+0.0j
  print("n=",n)
  print("A3=")
  print(A3)
  print("should be eigenvalues 1+4j, 2+3j, 3+2j, 4+1j in any order")
  print("calling  e3, V3 = eig(A3)")
  e3, V3 = eig(A3)
  print("back from eig")
  print("eigenvalues e3, of A3")
  print(e3)
  print("eigenvalues as columns of V3")
  print(V3)
  print(" ")
  check_eigen3(A3, e3, V3)
  print(" ")
  n = 4
  A4 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  e4 = array([(0.0+0.0j) for i in range(n)])
  V4 = array([[(0.0+0.0j) for j in range(n)] for i in range(n)])
  A4[0][0] = 10.0+14j    # for a=1+2j  b=2+3j  c=3+4j  d=4+5j
  A4[0][1] = 36.0-100j
  A4[0][2] = -270.0+66.0j
  A4[0][3] = 185.0+180.0j
  A4[1][0] = 1.0+0.0j
  A4[2][1] = 1.0+0.0j
  A4[3][2] = 1.0+0.0j
  print("n=",n)
  print("A4=")
  print(A4)
  print("should be eigenvalues 1+2j, 2+3j, 3+4j, 4+5j in any order")
  print("calling  e3, V4 = eig(A4)")
  e4, V4 = eig(A4)
  print("back from eig")
  print("eigenvalues e4, of A4")
  print(e4)
  print("eigenvalues as columns of V4")
  print(V4)
  print(" ")
  check_eigen3(A4, e4, V4)
  print(" ")

  print("matrixfromroots.py3 finished")

if __name__ == '__main__':
  main()