source('../Chapter3/utils.R')
source('utils.R')

A = matrix(data=c(4, -1, -1, -1, -1, 4, -1, -1, -1, -1, 4, -1, -1, -1, -1, 4), nrow=4, ncol=4, byrow=FALSE)
print('A=')
print(A)

res = householder_reduction(A, verbose=FALSE)
B = res$H
print('B= ')
print(round(B, 3))

eigenvalue_polynomial = function(x){
    res = evaluate_symmetric_tridiagonal_characteristic_polynomial(x, diag(B), diag(B[-nrow(B), -1]))
    res$ps[length(res$ps)] ## return the value of p_n(x)
}

## Lets get a feel for where the zeros of this polynomial i.e.:
##
xs = seq(0.5, +5.5, length.out=1000)
ys = sapply(xs, eigenvalue_polynomial)
plot(xs, ys, 'l', ylab='p(x)')
grid()

## Lets find all eigenvalues of B using Muller's method and the function "evaluate_symmetric_tridiagonal_characteristic_polynomial":
##
x0s = c(0.5, 2.5, 7.0)
res = mullers_iteration(x0s, eigenvalue_polynomial, 4, xtol=1.e-4, ftol=1.e-6, verbose=FALSE)
print('mullers roots=')
print(res$x)