forked from JuliaLinearAlgebra/ArnoldiMethod.jl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
eigvals.jl
45 lines (38 loc) · 1.16 KB
/
eigvals.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
"""
copy_eigenvalues!(λs, A) -> λs
Puts the eigenvalues of a quasi-upper triangular matrix A in the λs vector.
"""
function copy_eigenvalues!(λs, A::AbstractMatrix{T}, tol = eps(real(T))) where {T}
n = size(A, 1)
i = 1
@inbounds while i < n
if is_offdiagonal_small(A, i, tol)
λs[i] = A[i, i]
i += 1
else
# Conjugate pair
d = A[i,i] * A[i+1,i+1] - A[i,i+1] * A[i+1,i]
x = (A[i,i] + A[i+1,i+1]) / 2
y = sqrt(complex(x*x - d))
λs[i] = x + y
λs[i + 1] = x - y
i += 2
end
end
@inbounds if i == n
λs[i] = A[n, n]
end
return λs
end
"""
eigenvalues(A::AbstractMatrix{T}) -> Vector{complex(T)}
Computes the eigenvalues of the matrix A. Assumes that A is quasi-upper triangular.
The eigenvalues are returned in complex arithmetic, even if their imaginary
part is 0.
"""
eigenvalues(A::AbstractMatrix{T}, tol = eps(real(T))) where {T} =
copy_eigenvalues!(Vector{complex(T)}(undef, size(A, 1)), A, tol)
function partialeigen(P::PartialSchur)
vals, vecs = eigen(P.R)
return vals, P.Q*vecs
end