using LinearAlgebra
S=rand(3,3) .- 0.5
D=diagm([3,4,6])
A=S*D*inv(S)
eigvals(A)
x0=rand(3)
x=x0
for j=1:100
    global x,y
    y=A*x
    x=y/norm(y)
end
println("approximate eigenvector:")
display(x)
println("eigenvalue is ",x'*A*x/(x'*x))
B=A-6*I
x=x0
for j=1:100
    global x,y
    y=B*x 
    x=y/norm(y)
end    
println("\napproximate eigenvector:")
display(x)
beta=x'*B*x/(x'*x)
println("other eigenvalue ",beta+6)
B=inv(A-4.5*I)
x=x0
for j=1:100
    global x,y
    y=B*x
    x=y/norm(y)
end
println("\napproximate eigenvector:")
display(x)
beta=x'*B*x/(x'*x)
println("another eigenvalue ",1/beta+4.5)

println("\nShifted inverse iteration...")

function printeigen(lambda)
    println("\nInverse shift by $lambda")
    println("the eigenvalues of B are")
    println(1/(3-lambda))
    println(1/(4-lambda))
    println(1/(6-lambda))
end

B=inv(A-4.5*I)
printeigen(4.5)
y1=B*x0
x1=y1/norm(y1)
lambda=x1'*A*x1/(x1'*x1)
println("the new shift is ",lambda)

B=inv(A-lambda*I)
printeigen(lambda)
y2=B*x1
x2=y2/norm(y2)
lambda=x2'*A*x2/(x2'*x2)
println("the new shift is ",lambda)

B=inv(A-lambda*I)
printeigen(lambda)
y3=B*x2
x3=y3/norm(y3)
lambda=x3'*A*x3/(x3'*x3)
println("the new shift is ",lambda)