Skip to content

Commit d014c4f

Browse files
committed
Wrap LAPACK's non-symmetric eigensolver lahqr that only uses single
and double shifts.
1 parent db2b642 commit d014c4f

File tree

1 file changed

+57
-0
lines changed

1 file changed

+57
-0
lines changed

src/lapack.jl

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -210,6 +210,63 @@ module LAPACK2
210210
end
211211
end
212212

213+
# Non-symmetric eigenvalue problem
214+
function lahqr!(wantt::Bool,
215+
wantz::Bool,
216+
H::StridedMatrix{Float64},
217+
ilo::Integer,
218+
ihi::Integer,
219+
wr::StridedVector{Float64},
220+
wi::StridedVector{Float64},
221+
iloz::Integer,
222+
ihiz::Integer,
223+
Z::StridedMatrix{Float64})
224+
# SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
225+
# ILOZ, IHIZ, Z, LDZ, INFO )
226+
227+
# .. Scalar Arguments ..
228+
# INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
229+
# LOGICAL WANTT, WANTZ
230+
# ..
231+
# .. Array Arguments ..
232+
# DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
233+
234+
n = LinAlg.checksquare(H)
235+
ldh = stride(H, 2)
236+
ldz = stride(Z, 2)
237+
238+
info = Ref{BlasInt}(0)
239+
240+
ccall((@blasfunc("dlahqr_"), Base.liblapack_name), Void,
241+
(Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
242+
Ptr{BlasInt}, Ptr{Float64}, Ptr{BlasInt}, Ptr{Float64},
243+
Ptr{Float64}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{Float64},
244+
Ptr{BlasInt}, Ptr{BlasInt}),
245+
&wantt, &wantz, &n, &ilo,
246+
&ihi, H, &max(1, ldh), wr,
247+
wi, &iloz, &ihiz, Z,
248+
&max(1, ldz), info)
249+
250+
if info[] != 0
251+
throw(LAPACKException(info[]))
252+
end
253+
254+
return wr, wi, Z
255+
end
256+
function lahqr!(H::StridedMatrix{Float64})
257+
n = size(H, 1)
258+
return lahqr!(false,
259+
false,
260+
H,
261+
1,
262+
n,
263+
Vector{Float64}(n),
264+
Vector{Float64}(n),
265+
1,
266+
1,
267+
Matrix{Float64}(0, 0))
268+
end
269+
213270
## Cholesky + singular values
214271
function spteqr!(compz::Char, d::StridedVector{Float64}, e::StridedVector{Float64}, Z::StridedMatrix{Float64},
215272
work::StridedVector{Float64} = Vector{Float64}(4length(d)))

0 commit comments

Comments
 (0)