diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl
index dc7152a1..1033e5b5 100644
--- a/src/BandedMatrices.jl
+++ b/src/BandedMatrices.jl
@@ -3,27 +3,30 @@ using Base, FillArrays, ArrayLayouts, LinearAlgebra, SparseArrays, Random
 using LinearAlgebra.LAPACK
 import Base: axes, axes1, getproperty, iterate, tail
 import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!,
-                        copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym,
+                        copy_oftype, checksquare, AdjOrTrans, HermOrSym,
                         _chol!, rot180
 import LinearAlgebra.BLAS: libblas
 import LinearAlgebra.LAPACK: liblapack, chkuplo, chktrans
-import LinearAlgebra: cholesky, cholesky!, cholcopy, norm, diag, eigvals!, eigvals, eigen!, eigen,
+import LinearAlgebra: cholesky, cholesky!, cholcopy, norm, diag,
+            eigvals!, eigvals, eigen!, eigen, eigmax, eigmin, eigvecs,
             qr, qr!, axpy!, ldiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular,
             chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, logabsdet,
-            svdvals, svdvals!, QRPackedQ, checknonsingular, ipiv2perm, tril!,
-            triu!, Givens, diagzero
+            svdvals, svdvals!, QRPackedQ, checknonsingular, ipiv2perm,
+            tril!, triu!, istril, istriu, isdiag,
+            Givens, diagzero
 import SparseArrays: sparse
 
 import Base: getindex, setindex!, *, +, -, ==, <, <=, >, isassigned,
-                >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds
+                >=, /, ^, \, adjoint, transpose, showerror, reindex, checkbounds, @propagate_inbounds
 
 import Base: convert, size, view, unsafe_indices,
                 first, last, size, length, unsafe_length, step,
                 to_indices, to_index, show, fill!, promote_op,
                 MultiplicativeInverses, OneTo, ReshapedArray,
-                               similar, copy, convert, promote_rule, rand,
-                            IndexStyle, real, imag, Slice, pointer, unsafe_convert, copyto!,
-                            hcat, vcat, hvcat
+                similar, copy, convert, promote_rule, rand,
+                IndexStyle, conj, real, imag, Slice, pointer, unsafe_convert, copyto!,
+                hcat, vcat, hvcat,
+                iszero, isone
 
 import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted, broadcasted,
                         materialize, materialize!
@@ -32,9 +35,9 @@ import ArrayLayouts: MemoryLayout, transposelayout, triangulardata,
                     conjlayout, symmetriclayout, symmetricdata,
                     triangularlayout, MatLdivVec, hermitianlayout, hermitiandata,
                     materialize!, BlasMatMulMatAdd, BlasMatMulVecAdd, BlasMatLmulVec, BlasMatLdivVec,
-                    colsupport, rowsupport, symmetricuplo, MatMulMatAdd, MatMulVecAdd, 
+                    colsupport, rowsupport, symmetricuplo, MatMulMatAdd, MatMulVecAdd,
                     sublayout, sub_materialize, _fill_lmul!,
-                    reflector!, reflectorApply!, _copyto!, 
+                    reflector!, reflectorApply!, _copyto!,
                     _qr!, _qr, _lu!, _lu, _factorize, TridiagonalLayout
 
 import FillArrays: AbstractFill, getindex_value
@@ -87,6 +90,8 @@ include("banded/bandedqr.jl")
 include("banded/gbmm.jl")
 include("banded/linalg.jl")
 
+include("hermtridiag.jl")
+
 include("symbanded/symbanded.jl")
 include("symbanded/ldlt.jl")
 include("symbanded/BandedCholesky.jl")
diff --git a/src/hermtridiag.jl b/src/hermtridiag.jl
new file mode 100644
index 00000000..73cdda35
--- /dev/null
+++ b/src/hermtridiag.jl
@@ -0,0 +1,273 @@
+## Hermitian tridiagonal matrices
+struct HermTridiagonal{S, T, U <: AbstractVector{S}, V <: AbstractVector{T}} <: AbstractMatrix{T}
+    dv::U                        # diagonal
+    ev::V                        # superdiagonal
+    function HermTridiagonal{S, T, U, V}(dv, ev) where {S, T, U <: AbstractVector{S}, V <: AbstractVector{T}}
+        require_one_based_indexing(dv, ev)
+        if !(length(dv) - 1 <= length(ev) <= length(dv))
+            throw(DimensionMismatch("superdiagonal has wrong length. Has length $(length(ev)), but should be either $(length(dv) - 1) or $(length(dv))."))
+        end
+        new{S,T,U,V}(dv,ev)
+    end
+end
+
+HermTridiagonal(dv::U, ev::V) where {S,T,U<:AbstractVector{S},V<:AbstractVector{T}} = HermTridiagonal{S,T}(dv, ev)
+HermTridiagonal{S,T}(dv::U, ev::V) where {S,T,U<:AbstractVector{S},V<:AbstractVector{T}} = HermTridiagonal{S,T,U,V}(dv, ev)
+function HermTridiagonal{S,T}(dv::AbstractVector, ev::AbstractVector) where {S,T}
+    HermTridiagonal(convert(AbstractVector{S}, dv)::AbstractVector{S},
+                   convert(AbstractVector{T}, ev)::AbstractVector{T})
+end
+
+function HermTridiagonal(A::AbstractMatrix)
+    if diag(A,1) == conj(diag(A,-1))
+        HermTridiagonal(diag(A,0), diag(A,1))
+    else
+        throw(ArgumentError("matrix is not Hermitian; cannot convert to HermTridiagonal"))
+    end
+end
+
+HermTridiagonal{S,T,U,V}(H::HermTridiagonal{S,T,U,V}) where {S,T,U<:AbstractVector{S},V<:AbstractVector{T}} = H
+HermTridiagonal{S,T,U,V}(H::HermTridiagonal) where {S,T,U<:AbstractVector{S},V<:AbstractVector{T}} =
+    HermTridiagonal(convert(U, H.dv)::U, convert(V, H.ev)::V)
+HermTridiagonal{S,T}(H::HermTridiagonal{S,T}) where {S,T} = H
+HermTridiagonal{S,T}(H::HermTridiagonal) where {S,T} =
+    HermTridiagonal(convert(AbstractVector{S}, H.dv)::AbstractVector{S},
+                   convert(AbstractVector{T}, H.ev)::AbstractVector{T})
+HermTridiagonal(H::HermTridiagonal) = H
+
+AbstractMatrix{T}(H::HermTridiagonal) where {T} =
+    HermTridiagonal(convert(AbstractVector{T}, H.dv)::AbstractVector{T},
+                   convert(AbstractVector{T}, H.ev)::AbstractVector{T})
+function Matrix{T}(M::HermTridiagonal) where T
+    n = size(M, 1)
+    Mf = zeros(T, n, n)
+    @inbounds begin
+        @simd for i = 1:n-1
+            Mf[i,i] = M.dv[i]
+            Mf[i+1,i] = conj(M.ev[i])
+            Mf[i,i+1] = M.ev[i]
+        end
+        Mf[n,n] = M.dv[n]
+    end
+    return Mf
+end
+Matrix(M::HermTridiagonal{S,T}) where {S,T} = Matrix{T}(M)
+Array(M::HermTridiagonal) = Matrix(M)
+
+size(A::HermTridiagonal) = (length(A.dv), length(A.dv))
+function size(A::HermTridiagonal, d::Integer)
+    if d < 1
+        throw(ArgumentError("dimension must be ≥ 1, got $d"))
+    elseif d<=2
+        return length(A.dv)
+    else
+        return 1
+    end
+end
+
+# For S<:HermTridiagonal, similar(S[, neweltype]) should yield a HermTridiagonal matrix.
+# On the other hand, similar(S, [neweltype,] shape...) should yield a sparse matrix.
+# The first method below effects the former, and the second the latter.
+#similar(S::HermTridiagonal, ::Type{T}) where {T} = HermTridiagonal(similar(S.dv, T), similar(S.ev, T))
+# The method below is moved to SparseArrays for now
+# similar(S::HermTridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)
+
+#Elementary operations
+for func in (:conj, :copy, :real, :imag)
+    @eval ($func)(M::HermTridiagonal) = HermTridiagonal(($func)(M.dv), ($func)(M.ev))
+end
+
+transpose(H::HermTridiagonal) = Transpose(H)
+adjoint(H::HermTridiagonal) = H
+Base.copy(S::Adjoint{<:Any,<:HermTridiagonal}) = HermTridiagonal(map(x -> copy.(adjoint.(x)), (S.parent.dv, S.parent.ev))...)
+Base.copy(S::Transpose{<:Any,<:HermTridiagonal}) = HermTridiagonal(map(x -> copy.(transpose.(x)), (S.parent.dv, S.parent.ev))...)
+
+#=
+function diag(M::HermTridiagonal, n::Integer=0)
+    # every branch call similar(..., ::Int) to make sure the
+    # same vector type is returned independent of n
+    absn = abs(n)
+    if absn == 0
+        return copyto!(similar(M.dv, length(M.dv)), M.dv)
+    elseif absn==1
+        return copyto!(similar(M.ev, length(M.ev)), M.ev)
+    elseif absn <= size(M,1)
+        return fill!(similar(M.dv, size(M,1)-absn), 0)
+    else
+        throw(ArgumentError(string("requested diagonal, $n, must be at least $(-size(M, 1)) ",
+            "and at most $(size(M, 2)) for an $(size(M, 1))-by-$(size(M, 2)) matrix")))
+    end
+end
+=#
+
++(A::HermTridiagonal, B::HermTridiagonal) = HermTridiagonal(A.dv+B.dv, A.ev+B.ev)
+-(A::HermTridiagonal, B::HermTridiagonal) = HermTridiagonal(A.dv-B.dv, A.ev-B.ev)
+*(A::HermTridiagonal, B::Number) = HermTridiagonal(A.dv*B, A.ev*B)
+*(B::Number, A::HermTridiagonal) = A*B
+/(A::HermTridiagonal, B::Number) = HermTridiagonal(A.dv/B, A.ev/B)
+==(A::HermTridiagonal, B::HermTridiagonal) = (A.dv==B.dv) && (A.ev==B.ev)
+
+#=
+@inline mul!(A::StridedVecOrMat, B::HermTridiagonal, C::StridedVecOrMat,
+             alpha::Number, beta::Number) =
+    _mul!(A, B, C, MulAddMul(alpha, beta))
+
+@inline function _mul!(C::StridedVecOrMat, S::HermTridiagonal, B::StridedVecOrMat,
+                          _add::MulAddMul)
+    m, n = size(B, 1), size(B, 2)
+    if !(m == size(S, 1) == size(C, 1))
+        throw(DimensionMismatch("A has first dimension $(size(S,1)), B has $(size(B,1)), C has $(size(C,1)) but all must match"))
+    end
+    if n != size(C, 2)
+        throw(DimensionMismatch("second dimension of B, $n, doesn't match second dimension of C, $(size(C,2))"))
+    end
+
+    if m == 0
+        return C
+    elseif iszero(_add.alpha)
+        return _rmul_or_fill!(C, _add.beta)
+    end
+
+    α = S.dv
+    β = S.ev
+    @inbounds begin
+        for j = 1:n
+            x₊ = B[1, j]
+            x₀ = zero(x₊)
+            # If m == 1 then β[1] is out of bounds
+            β₀ = m > 1 ? zero(β[1]) : zero(eltype(β))
+            for i = 1:m - 1
+                x₋, x₀, x₊ = x₀, x₊, B[i + 1, j]
+                β₋, β₀ = β₀, β[i]
+                _modify!(_add, β₋*x₋ + α[i]*x₀ + β₀*x₊, C, (i, j))
+            end
+            _modify!(_add, β₀*x₀ + α[m]*x₊, C, (m, j))
+        end
+    end
+
+    return C
+end
+
+(\)(T::HermTridiagonal, B::StridedVecOrMat) = ldlt(T)\B
+
+# division with optional shift for use in shifted-Hessenberg solvers (hessenberg.jl):
+ldiv!(A::HermTridiagonal, B::AbstractVecOrMat; shift::Number=false) = ldiv!(ldlt(A, shift=shift), B)
+rdiv!(B::AbstractVecOrMat, A::HermTridiagonal; shift::Number=false) = rdiv!(B, ldlt(A, shift=shift))
+=#
+
+
+function eigen!(A::HermTridiagonal{S,T}) where {S,T}
+    n = size(A, 1)
+    D = ones(T, n)
+    d = copy(A.dv)
+    e = zeros(S, n-1)
+    for i in 1:n-1
+        e[i] = abs(A.ev[i])
+        if e[i] != zero(S)
+            D[i+1] = A.ev[i]/e[i]
+        end
+        if i < n-1
+            A.ev[i+1] = A.ev[i+1]*D[i+1]
+        end
+    end
+    Λ, V = eigen(SymTridiagonal(d, e))
+    return Eigen(Λ, Diagonal(conj(D))*V)
+end
+eigen(A::HermTridiagonal{S,T}) where {S,T} = eigen!(copy(A))
+
+eigvals!(A::HermTridiagonal{S,T}, kwargs...) where {S,T} = eigvals!(SymTridiagonal(A.dv, map(abs, A.ev)), kwargs...)
+eigvals(A::HermTridiagonal{S,T}, kwargs...) where {S,T} = eigvals!(copy(A), kwargs...)
+
+#Computes largest and smallest eigenvalue
+eigmax(A::HermTridiagonal) = eigvals(A, size(A, 1):size(A, 1))[1]
+eigmin(A::HermTridiagonal) = eigvals(A, 1:1)[1]
+
+#Compute selected eigenvectors only corresponding to particular eigenvalues
+eigvecs(A::HermTridiagonal) = eigen(A).vectors
+
+function svdvals!(A::HermTridiagonal)
+    vals = eigvals!(A)
+    return sort!(map!(abs, vals, vals); rev=true)
+end
+
+#tril and triu
+
+istriu(M::HermTridiagonal) = iszero(M.ev)
+istril(M::HermTridiagonal) = iszero(M.ev)
+iszero(M::HermTridiagonal) = iszero(M.ev) && iszero(M.dv)
+isone(M::HermTridiagonal) = iszero(M.ev) && all(isone, M.dv)
+isdiag(M::HermTridiagonal) = iszero(M.ev)
+
+function tril!(M::HermTridiagonal, k::Integer=0)
+    n = length(M.dv)
+    if !(-n - 1 <= k <= n - 1)
+        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
+            "$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix")))
+    elseif k < -1
+        fill!(M.ev,0)
+        fill!(M.dv,0)
+        return Tridiagonal(conj(M.ev),M.dv,copy(M.ev))
+    elseif k == -1
+        fill!(M.dv,0)
+        return Tridiagonal(conj(M.ev),M.dv,zero(M.ev))
+    elseif k == 0
+        return Tridiagonal(conj(M.ev),M.dv,zero(M.ev))
+    elseif k >= 1
+        return Tridiagonal(conj(M.ev),M.dv,copy(M.ev))
+    end
+end
+
+function triu!(M::HermTridiagonal, k::Integer=0)
+    n = length(M.dv)
+    if !(-n + 1 <= k <= n + 1)
+        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
+            "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix")))
+    elseif k > 1
+        fill!(M.ev,0)
+        fill!(M.dv,0)
+        return Tridiagonal(conj(M.ev),M.dv,copy(M.ev))
+    elseif k == 1
+        fill!(M.dv,0)
+        return Tridiagonal(zero(M.ev),M.dv,M.ev)
+    elseif k == 0
+        return Tridiagonal(zero(M.ev),M.dv,M.ev)
+    elseif k <= -1
+        return Tridiagonal(conj(M.ev),M.dv,copy(M.ev))
+    end
+end
+
+###################
+# Generic methods #
+###################
+
+## structured matrix methods ##
+function Base.replace_in_print_matrix(A::HermTridiagonal, i::Integer, j::Integer, s::AbstractString)
+    i==j-1||i==j||i==j+1 ? s : Base.replace_with_centered_mark(s)
+end
+
+#logabsdet(A::HermTridiagonal; shift::Number=false) = logabsdet(ldlt(A; shift=shift))
+
+function getindex(A::HermTridiagonal{T}, i::Integer, j::Integer) where T
+    if !(1 <= i <= size(A,2) && 1 <= j <= size(A,2))
+        throw(BoundsError(A, (i,j)))
+    end
+    if i == j
+        return A.dv[i]
+    elseif i == j + 1
+        return conj(A.ev[j])
+    elseif i + 1 == j
+        return A.ev[i]
+    else
+        return zero(T)
+    end
+end
+
+function setindex!(A::HermTridiagonal, x, i::Integer, j::Integer)
+    @boundscheck checkbounds(A, i, j)
+    if i == j
+        @inbounds A.dv[i] = x
+    else
+        throw(ArgumentError("cannot set off-diagonal entry ($i, $j)"))
+    end
+    return x
+end
diff --git a/src/lapack.jl b/src/lapack.jl
index 5081edd0..66beaaad 100644
--- a/src/lapack.jl
+++ b/src/lapack.jl
@@ -16,8 +16,24 @@ for (fname, elty) in ((:dlargv_,:Float64),
         #       INTEGER INCC, INCX, INCY, N
         #       .. Array Arguments ..
         #       DOUBLE PRECISION C( * ), X( * ), Y( * )
-        function largv!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::Ptr{$elty}, incc::Int)
-            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, x, incx, y, incy, c, incc)
+        function largv!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::Ptr{real($elty)}, incc::Int)
+            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{real($elty)}, Ref{BlasInt}), n, x, incx, y, incy, c, incc)
+        end
+    end
+end
+
+function lacgv!(n::Int, x::Ptr{T}, incx::Int) where T <: Real end
+
+for (fname, elty) in ((:zlacgv_,:ComplexF64),
+                      (:clacgv_,:ComplexF32))
+    @eval begin
+        #       SUBROUTINE ZLACGV( N, X, INCX )
+        #       .. Scalar Arguments ..
+        #       INTEGER INCX, N
+        #       .. Array Arguments ..
+        #       COMPLEX*16 X( * )
+        function lacgv!(n::Int, x::Ptr{$elty}, incx::Int)
+            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, x, incx)
         end
     end
 end
@@ -32,8 +48,8 @@ for (fname, elty) in ((:dlartv_,:Float64),
         #       INTEGER INCC, INCX, INCY, N
         #       .. Array Arguments ..
         #       DOUBLE PRECISION C( * ), S( * ), X( * ), Y( * )
-        function lartv!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::Ptr{$elty}, s::Ptr{$elty}, incc::Int)
-            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}), n, x, incx, y, incy, c, s, incc)
+        function lartv!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::Ptr{real($elty)}, s::Ptr{$elty}, incc::Int)
+            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{real($elty)}, Ptr{$elty}, Ref{BlasInt}), n, x, incx, y, incy, c, s, incc)
         end
     end
 end
@@ -49,8 +65,8 @@ for (fname, elty) in ((:drot_,:Float64),
         #       INTEGER INCX,INCY,N
         #       .. Array Arguments ..
         #       DOUBLE PRECISION DX(*),DY(*)
-        function rot!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::$elty, s::$elty)
-            ccall((@blasfunc($fname), libblas), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ref{$elty}), n, x, incx, y, incy, c, s)
+        function rot!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::real($elty), s::$elty)
+            ccall((@blasfunc($fname), libblas), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{real($elty)}, Ref{$elty}), n, x, incx, y, incy, c, s)
         end
     end
 end
@@ -63,8 +79,8 @@ for (fname, elty) in ((:dlartg_,:Float64),
         #       SUBROUTINE DLARTG( F, G, CS, SN, R )
         #       .. Scalar Arguments ..
         #       DOUBLE PRECISION CS, F, G, R, SN
-        function lartg!(f::$elty, g::$elty, cs::Ref{$elty}, sn::Ref{$elty}, r::Ref{$elty})
-            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{$elty}, Ref{$elty}, Ref{$elty}, Ref{$elty}, Ref{$elty}), f, g, cs, sn, r)
+        function lartg!(f::$elty, g::$elty, cs::Ref{real($elty)}, sn::Ref{$elty}, r::Ref{$elty})
+            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{$elty}, Ref{$elty}, Ref{real($elty)}, Ref{$elty}, Ref{$elty}), f, g, cs, sn, r)
         end
     end
 end
@@ -79,8 +95,8 @@ for (fname, elty) in ((:dlar2v_,:Float64),
         #       INTEGER INCC, INCX, N
         #       .. Array Arguments ..
         #       DOUBLE PRECISION C( * ), S( * ), X( * ), Y( * ), Z( * )
-        function lar2v!(n::Int, x::Ptr{$elty}, y::Ptr{$elty}, z::Ptr{$elty}, incx::Int, c::Ptr{$elty}, s::Ptr{$elty}, incc::Int)
-            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}), n, x, y, z, incx, c, s, incc)
+        function lar2v!(n::Int, x::Ptr{$elty}, y::Ptr{$elty}, z::Ptr{$elty}, incx::Int, c::Ptr{real($elty)}, s::Ptr{$elty}, incc::Int)
+            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{real($elty)}, Ptr{$elty}, Ref{BlasInt}), n, x, y, z, incx, c, s, incc)
         end
     end
 end
diff --git a/src/symbanded/bandedeigen.jl b/src/symbanded/bandedeigen.jl
index b90edb7e..d4c7b127 100644
--- a/src/symbanded/bandedeigen.jl
+++ b/src/symbanded/bandedeigen.jl
@@ -17,16 +17,13 @@ end
 size(B::BandedGeneralizedEigenvectors) = size(B.W)
 getindex(B::BandedGeneralizedEigenvectors, i, j) = Matrix(B)[i, j]
 
-convert(::Type{Eigen{T, T, Matrix{T}, Vector{T}}}, F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T = Eigen(F.values, Matrix(F.vectors))
-convert(::Type{GeneralizedEigen{T, T, Matrix{T}, Vector{T}}}, F::GeneralizedEigen{T, T, BandedGeneralizedEigenvectors{T}, Vector{T}}) where T = GeneralizedEigen(F.values, Matrix(F.vectors))
+convert(::Type{Eigen{S, T, Matrix{S}, Vector{T}}}, F::Eigen{S, T, BandedEigenvectors{S}, Vector{T}}) where {S, T} = Eigen(F.values, Matrix(F.vectors))
+convert(::Type{GeneralizedEigen{S, T, Matrix{S}, Vector{T}}}, F::GeneralizedEigen{S, T, BandedGeneralizedEigenvectors{S}, Vector{T}}) where {S, T} = GeneralizedEigen(F.values, Matrix(F.vectors))
 
-compress(F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T = convert(Eigen{T, T, Matrix{T}, Vector{T}}, F)
-compress(F::GeneralizedEigen{T, T, BandedGeneralizedEigenvectors{T}, Vector{T}}) where T = convert(GeneralizedEigen{T, T, Matrix{T}, Vector{T}}, F)
+compress(F::Eigen{S, T, BandedEigenvectors{S}, Vector{T}}) where {S, T} = convert(Eigen{S, T, Matrix{S}, Vector{T}}, F)
+compress(F::GeneralizedEigen{S, T, BandedGeneralizedEigenvectors{S}, Vector{T}}) where {S, T} = convert(GeneralizedEigen{S, T, Matrix{S}, Vector{T}}, F)
 
-eigen(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigen!(copy(A))
-eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigen!(copy(A), copy(B))
-
-function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
+function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, args...) where T <: Real
     N = size(A, 1)
     KD = bandwidth(A)
     D = Vector{T}(undef, N)
@@ -34,12 +31,56 @@ function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
     G = Vector{Givens{T}}(undef, 0)
     WORK = Vector{T}(undef, N)
     AB = symbandeddata(A)
-    sbtrd!('V', A.uplo, N, KD, AB, D, E, G, WORK)
-    Λ, Q = eigen(SymTridiagonal(D, E))
+    sbtrd!('N', symmetricuplo(A), N, KD, AB, D, E, G, WORK)
+    return eigvals!(SymTridiagonal(D, E), args...)
+end
+
+function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, args...; kwds...) where T <: Real
+    N = size(A, 1)
+    KD = bandwidth(A)
+    NG = count_givens_sbtrd(N, KD)
+    D = Vector{T}(undef, N)
+    E = Vector{T}(undef, N-1)
+    G = Vector{Givens{T}}(undef, NG)
+    WORK = Vector{T}(undef, N)
+    AB = symbandeddata(A)
+    sbtrd!('V', symmetricuplo(A), N, KD, AB, D, E, G, WORK)
+    Λ, Q = eigen(SymTridiagonal(D, E), args...)
+    Eigen(Λ, BandedEigenvectors(G, Q))
+end
+
+function eigvals!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex
+    N = size(A, 1)
+    KD = bandwidth(A)
+    D = Vector{real(T)}(undef, N)
+    E = Vector{T}(undef, N-1)
+    G = Vector{Givens{T}}(undef, 0)
+    WORK = Vector{T}(undef, N)
+    AB = hermbandeddata(A)
+    sbtrd!('N', symmetricuplo(A), N, KD, AB, D, E, G, WORK)
+    return eigvals!(HermTridiagonal(D, E), args...)
+end
+eigvals!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T = eigvals!(Symmetric(A), args...)
+
+function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...; kwds...) where T <: Complex
+    N = size(A, 1)
+    KD = bandwidth(A)
+    NG = count_givens_sbtrd(N, KD)
+    D = Vector{real(T)}(undef, N)
+    E = Vector{T}(undef, N-1)
+    G = Vector{Givens{T}}(undef, NG)
+    WORK = Vector{T}(undef, N)
+    AB = hermbandeddata(A)
+    sbtrd!('V', symmetricuplo(A), N, KD, AB, D, E, G, WORK)
+    if symmetricuplo(A) == 'L'
+        conj!(E)
+    end
+    Λ, Q = eigen(HermTridiagonal(D, E), args...)
     Eigen(Λ, BandedEigenvectors(G, Q))
 end
+eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...; kwds...) where T = eigen!(Symmetric(A), args...; kwds...)
 
-function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
+function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}, args...) where T <: Real
     S = splitcholesky!(B)
     N = size(A, 1)
     KA = bandwidth(A)
@@ -48,10 +89,24 @@ function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix
     WORK = Vector{T}(undef, 2*N)
     AB = symbandeddata(A)
     BB = symbandeddata(B)
-    sbgst!('V', A.uplo, N, KA, KB, AB, BB, Q, WORK)
-    Λ, W = eigen!(A)
+    sbgst!('N', symmetricuplo(A), N, KA, KB, AB, BB, Q, WORK)
+    return eigvals!(A, args...)
+end
+
+function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}, args...; kwds...) where T <: Real
+    S = splitcholesky!(B)
+    N = size(A, 1)
+    KA = bandwidth(A)
+    KB = bandwidth(B)
+    Q = Vector{Givens{T}}(undef, 0)
+    WORK = Vector{T}(undef, 2*N)
+    AB = symbandeddata(A)
+    BB = symbandeddata(B)
+    sbgst!('V', symmetricuplo(A), N, KA, KB, AB, BB, Q, WORK)
+    Λ, W = eigen!(A, args...; kwds...)
     GeneralizedEigen(Λ, BandedGeneralizedEigenvectors(S, Q, W))
 end
+eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}, args...; kwds...) where T <: Real = eigen!(copy(A), copy(B), args...; kwds...)
 
 function Matrix(B::BandedEigenvectors)
     Q = copy(B.Q)
@@ -73,7 +128,7 @@ function Matrix(B::BandedGeneralizedEigenvectors)
 end
 
 
-function compress!(F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T
+function compress!(F::Eigen{S, T, BandedEigenvectors{S}, Vector{T}}) where {S, T}
     Q = F.vectors.Q
     G = F.vectors.G
     for k in length(G):-1:1
@@ -165,6 +220,29 @@ function mul!(y::Array{T,N}, x::Array{T,N}, B::BandedGeneralizedEigenvectors{T})
     y .= y'
 end
 
+Givens(i1::Int, i2::Int, c::S, s::T) where {S, T} = Givens(i1, i2, promote(c, s)...)
+
+function count_givens_sbtrd(N::Integer, KD::Integer)
+    KD1 = KD + 1
+    KDN = min(N-1, KD)
+    NG = 0
+    J1 = KDN + 2
+    J2 = 1
+    for I = 1:N-2
+        for K = KDN+1:-1:2
+            J1 = J1 + KDN
+            J2 = J2 + KDN
+            if K > 2
+                J1 = J1 - KDN - 1
+            end
+            NG += length(J1:KD1:J2)
+            if J2+KDN > N
+                J2 = J2 - KDN - 1
+            end
+        end
+    end
+    NG
+end
 
 
 #
@@ -187,8 +265,8 @@ end
 #       ..
 function sbtrd!(VECT::Char, UPLO::Char,
                 N::Int, KD::Int, AB::AbstractMatrix{T},
-                D::AbstractVector{T}, E::AbstractVector{T}, Q::Vector{Givens{T}},
-                WORK::AbstractVector{T}) where T
+                D::AbstractVector{S}, E::AbstractVector{T}, Q::Vector{Givens{T}},
+                WORK::AbstractVector{T}) where {S, T}
     require_one_based_indexing(AB)
     chkstride1(AB)
     chkuplo(UPLO)
@@ -206,7 +284,7 @@ function sbtrd!(VECT::Char, UPLO::Char,
     LDABF = max(1, stride(AB, 2))
 
     ZERO = zero(T)
-    TEMP1 = Ref{T}()
+    TEMP1 = Ref{S}()
     TEMP2 = Ref{T}()
     TEMP3 = Ref{T}()
     KD1 = KD + 1
@@ -217,7 +295,8 @@ function sbtrd!(VECT::Char, UPLO::Char,
     KDN = min(N-1, KD)
 
     UPPER = UPLO == 'U'
-    WANTQ = true
+    WANTQ = VECT != 'N'
+    NG = 1
 
     if UPPER
         if KD > 1
@@ -271,6 +350,7 @@ function sbtrd!(VECT::Char, UPLO::Char,
                     end
                     # apply plane rotations from the left
                     if NR > 0
+                        lacgv!(NR, pointer(WORK, J1), KD1)
                         if 2*KD-1 < NR
                             # Dependent on the the number of diagonals either
                             # DLARTV or DROT is used
@@ -291,7 +371,7 @@ function sbtrd!(VECT::Char, UPLO::Char,
                                     rot!(KDM1, pointer(AB, KDM1+LDAB*JIN), INCX, pointer(AB, KD+LDAB*JIN), INCX, D[JIN], WORK[JIN])
                                 end
                             end
-                            LEND = min( KDM1, N-J2 )
+                            LEND = min(KDM1, N-J2)
                             LAST = J1END + KD1
                             if LEND > 0
                                 rot!(LEND, pointer(AB, KDM1+LDAB*LAST), INCX, pointer(AB, KD+LDAB*LAST), INCX, D[LAST], WORK[LAST])
@@ -301,9 +381,11 @@ function sbtrd!(VECT::Char, UPLO::Char,
 
                     if WANTQ
                         for J = J1:KD1:J2
-                            push!(Q, Givens(J-1, J, D[J], -WORK[J]))
+                            Q[NG] = Givens(J-1, J, D[J], -WORK[J])
+                            NG += 1
                         end
                     end
+
                     if J2+KDN > N
                         # adjust J2 to keep within the bounds of the matrix
                         NR = NR - 1
@@ -333,7 +415,7 @@ function sbtrd!(VECT::Char, UPLO::Char,
         end
         # copy diagonal elements to D
         for I = 1:N
-            D[I] = AB[KD1, I]
+            D[I] = real(AB[KD1, I])
         end
     else # if UPPER
         if KD > 1
@@ -390,6 +472,7 @@ function sbtrd!(VECT::Char, UPLO::Char,
                     # DLARTV or DROT is used
 
                     if NR > 0
+                        lacgv!(NR, pointer(WORK, J1), KD1)
                         if NR > 2*KD-1
                             for L = 1:KD-1
                                 if J2+L > N
@@ -418,7 +501,8 @@ function sbtrd!(VECT::Char, UPLO::Char,
 
                     if WANTQ
                         for J = J1:KD1:J2
-                            push!(Q, Givens(J-1, J, D[J], -WORK[J]))
+                            Q[NG] = Givens(J-1, J, D[J], conj(-WORK[J]))
+                            NG += 1
                         end
                     end
 
@@ -437,6 +521,7 @@ function sbtrd!(VECT::Char, UPLO::Char,
                 end
             end
         end # if
+
         if KD > 0
             # copy off-diagonal elements to E
             for I = 1:N-1
@@ -450,7 +535,7 @@ function sbtrd!(VECT::Char, UPLO::Char,
         end
         # copy diagonal elements to D
         for I = 1:N
-            D[I] = AB[1, I]
+            D[I] = real(AB[1, I])
         end
     end # if UPPER
 end
@@ -496,7 +581,6 @@ function sbgst!(VECT::Char, UPLO::Char,
     LDAB = size(AB, 1)
     LDABF = max(1, stride(AB, 2))
 
-    ZERO = zero(T)
     TEMP1 = Ref{T}()
     TEMP2 = Ref{T}()
     TEMP3 = Ref{T}()
@@ -506,7 +590,7 @@ function sbgst!(VECT::Char, UPLO::Char,
     M = ( N+KB ) ÷ 2
 
     UPPER = UPLO == 'U'
-    WANTX = true
+    WANTX = VECT != 'N'
 
     UPDATE = true
     I = N + 1
diff --git a/src/symbanded/symbanded.jl b/src/symbanded/symbanded.jl
index c871a2f1..54e56719 100644
--- a/src/symbanded/symbanded.jl
+++ b/src/symbanded/symbanded.jl
@@ -101,43 +101,3 @@ function materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:BandedColumnMajor}
     l ≥ 0 || return lmul!(β, y)
     _banded_hbmv!(symmetricuplo(A), α, A, x, β, y)
 end
-
-
-## eigvals routine
-
-
-function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T
-    n=size(A, 1)
-    d = Vector{T}(undef,n)
-    e = Vector{T}(undef,n-1)
-    Q = Matrix{T}(undef,0,0)
-    work = Vector{T}(undef,n)
-    sbtrd!('N', symmetricuplo(A), size(A,1), bandwidth(A), symbandeddata(A), d, e, Q, work)
-    SymTridiagonal(d,e)
-end
-
-tridiagonalize!(A::AbstractMatrix) = _tridiagonalize!(A, MemoryLayout(typeof(A)))
-tridiagonalize(A::AbstractMatrix) = tridiagonalize!(copy(A))
-
-eigvals!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigvals!(tridiagonalize!(A))
-eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigvals!(copy(A))
-
-function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T<:Real
-    n = size(A, 1)
-    @assert n == size(B, 1)
-    @assert A.uplo == B.uplo
-    # compute split-Cholesky factorization of B.
-    kb = bandwidth(B)
-    B_data = symbandeddata(B)
-    pbstf!(B.uplo, n, kb, B_data)
-    # convert to a regular symmetric eigenvalue problem.
-    ka = bandwidth(A)
-    A_data = symbandeddata(A)
-    X = Array{T}(undef,0,0)
-    work = Vector{T}(undef,2n)
-    sbgst!('N', A.uplo, n, ka, kb, A_data, B_data, X, work)
-    # compute eigenvalues of symmetric eigenvalue problem.
-    eigvals!(A)
-end
-
-eigvals(A::Symmetric{<:Any,<:BandedMatrix}, B::Symmetric{<:Any,<:BandedMatrix}) = eigvals!(copy(A), copy(B))
diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl
index af573983..b943ffa4 100644
--- a/test/test_symbanded.jl
+++ b/test/test_symbanded.jl
@@ -38,15 +38,17 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol
     # (generalized) eigen & eigvals
     Random.seed!(0)
 
-    A = Symmetric(brand(Float64, 100, 100, 2, 4))
-    @test eigvals(A) ≈ eigvals(Symmetric(Matrix(A)))
-
-    F = eigen(A)
-    Λ, Q = F
-    @test Q'Matrix(A)*Q ≈ Diagonal(Λ)
-    FD = convert(Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}, F)
-    @test FD.vectors'Matrix(A)*FD.vectors ≈ Diagonal(F.values)
-
+    for T in (Float32, Float64)
+        for uplo in (:U, :L)
+            A = Symmetric(brand(T, 100, 100, 2, 4), uplo)
+            @test eigvals(A) ≈ eigvals(Symmetric(Matrix(A)))
+            F = eigen(A)
+            Λ, Q = F
+            @test Q'Matrix(A)*Q ≈ Diagonal(Λ)
+            ΛD, QD = convert(Eigen{T, T, Matrix{T}, Vector{T}}, F)
+            @test QD'Matrix(A)*QD ≈ Diagonal(ΛD)
+        end
+    end
 
     function An(::Type{T}, N::Int) where {T}
         A = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2)))
@@ -115,6 +117,33 @@ end
     @test all(A*x .=== muladd!(one(T),A,x,zero(T),copy(x)) .===
                 materialize!(MulAdd(one(T),A,x,zero(T),copy(x))) .===
                 BLAS.hbmv!('L', 1, one(T), view(parent(A).data,3:4,:), x, zero(T), similar(x)))
+
+    # (generalized) eigen & eigvals
+    Random.seed!(0)
+
+    for T in (Float32, Float64)
+        for uplo in (:U, :L)
+            A = Hermitian(brand(T, 100, 100, 2, 4), uplo)
+            @test eigvals(A) == eigvals(Symmetric(A.data, uplo))
+            F = eigen(A)
+            Λ, Q = F
+            @test Q'Matrix(A)*Q ≈ Diagonal(Λ)
+            ΛD, QD = convert(Eigen{T, T, Matrix{T}, Vector{T}}, F)
+            @test QD'Matrix(A)*QD ≈ Diagonal(ΛD)
+        end
+    end
+
+    for T in (ComplexF32, ComplexF64)
+        for uplo in (:U, :L)
+            A = Hermitian(brand(T, 100, 100, 2, 4), uplo)
+            @test eigvals(A) ≈ eigvals(Hermitian(Matrix(A), uplo))
+            F = eigen(A)
+            Λ, Q = F
+            @test Q'Matrix(A)*Q ≈ Diagonal(Λ)
+            ΛD, QD = convert(Eigen{T, real(T), Matrix{T}, Vector{real(T)}}, F)
+            @test QD'Matrix(A)*QD ≈ Diagonal(ΛD)
+        end
+    end
 end
 
 @testset "LDLᵀ" begin
@@ -139,7 +168,7 @@ end
         y = F\b
         @test x ≈ y
         @test_throws DimensionMismatch F\[b;b]
-        
+
         T ≠ Float16 && (@test det(F) ≈ det(SAL))
     end
     for T in (Int16, Int32, Int64, BigInt)