# +----------------------------------------------------+
# |                  Lekce 10 - demo 1                 |
# +----------------------------------------------------+
#  Linearni algebra - numpy
#   diagonalizace pomoci numpy.linalg

import numpy as np
#import scipy.linalg as la
import matplotlib.pyplot as plt
  
#Testik na matici jejiz diagonalizaci zname
N=50
q=np.linalg.qr(np.random.rand(N,N))[0]  #konstrukce nahodne ortogonalni matice pomoci qr-faktorizace
d=np.random.rand(N)
d.sort()                                #zavolani tridici metody
a=np.dot(np.dot(q,np.diag(d)),q.transpose())

#testik: diagonalizace rutionou z numpy:
d1=np.linalg.eigvalsh(a)
print('max. rozdil vl. cisel: ',max(abs(d1-d)))

#Pomocna procedura: Givensova rotace - modifikuje matici
#  -> nesmi byt i=j (netestuji)
#  -> tohle je specialni pripad, kdy phi je voleno, aby a[i,j] vynulovalo
def Givens(a,p,q):
  n=a.shape[0]
  if (abs(a[p,q])>1e-20):
    theta=0.5*(a[q,q]-a[p,p])/a[p,q]
    t=1.0/(abs(theta)+np.sqrt(1.0+theta**2))
    if (theta<0.0):
      t=-t
    c=1.0/np.sqrt(1+t**2)
    s=t*c
    tau=s/(1+c)
    a[p,p], a[q,q] = a[p,p]-t*a[p,q], a[q,q]+t*a[p,q]
    a[p,q], a[q,p] = 0.0, 0.0
    for r in range(n):
      if ((r!=p)and(r!=q)):
        a[r,p], a[r,q] = a[r,p]-s*(a[r,q]+tau*a[r,p]), a[r,q]+s*(a[r,p]-tau*a[r,q]) 
        a[q,r]=a[r,q]
        a[p,r]=a[r,p]

#Testik, ze Givens nezmeni spektrum
Givens(a,2,5)
d2=np.linalg.eigvalsh(a)
print('max. rozdil vl. cisel: ',max(abs(d2-d)))