Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RFC: Rework copy_oftype a bit #44756

Merged
merged 2 commits into from
Mar 30, 2022
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 11 additions & 10 deletions stdlib/LinearAlgebra/src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -355,36 +355,37 @@ control over the factorization of `B`.
"""
rdiv!(A, B)


copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)

"""
copy_oftype(A, T)
copymutable_oftype(A, T)
dkarrasch marked this conversation as resolved.
Show resolved Hide resolved

Copy `A` to a mutable array with eltype `T` based on `similar(A, T)`.

The resulting matrix typically has similar algebraic structure as `A`. For
example, supplying a tridiagonal matrix results in another tridiagonal matrix.
In general, the type of the output corresponds to that of `similar(A, T)`.

There are three often used methods in LinearAlgebra to create a mutable copy
of an array with a given eltype. These copies can be passed to in-place
algorithms (such as `ldiv!`, `rdiv!`, `lu!` and so on). Which one to use in practice
depends on what is known (or assumed) about the structure of the array in that
algorithm.
In LinearAlgebra, mutable copies (of some desired eltype) are created to be passed
to in-place algorithms (such as `ldiv!`, `rdiv!`, `lu!` and so on). If the specific
algorithm is known to preserve the algebraic structure, use `copymutable_oftype`.
If the algorithm is known to return a dense matrix (or some wrapper backed by a dense
matrix), then use `copy_similar`.

See also: `copy_similar`.
"""
copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T), A)
copymutable_oftype(A::AbstractArray, ::Type{S}) where {S} = copyto!(similar(A, S), A)

"""
copy_similar(A, T)

Copy `A` to a mutable array with eltype `T` based on `similar(A, T, size(A))`.

Compared to `copy_oftype`, the result can be more flexible. In general, the type
Compared to `copymutable_oftype`, the result can be more flexible. In general, the type
of the output corresponds to that of the three-argument method `similar(A, T, size(A))`.

See also: `copy_oftype`.
See also: `copymutable_oftype`.
"""
copy_similar(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T, size(A)), A)

Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/bunchkaufman.jl
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@ julia> S.L*S.D*S.L' - A[S.p, S.p]
```
"""
bunchkaufman(A::AbstractMatrix{T}, rook::Bool=false; check::Bool = true) where {T} =
bunchkaufman!(copy_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
bunchkaufman!(copymutable_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)

BunchKaufman{T}(B::BunchKaufman) where {T} =
BunchKaufman(convert(Matrix{T}, B.LD), B.ipiv, B.uplo, B.symmetric, B.rook, B.info)
Expand Down
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/src/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -179,8 +179,8 @@ Base.iterate(C::CholeskyPivoted, ::Val{:done}) = nothing

# make a copy that allow inplace Cholesky factorization
@inline choltype(A) = promote_type(typeof(sqrt(oneunit(eltype(A)))), Float32)
@inline cholcopy(A::StridedMatrix) = copy_oftype(A, choltype(A))
@inline cholcopy(A::RealHermSymComplexHerm) = copy_oftype(A, choltype(A))
@inline cholcopy(A::StridedMatrix) = copymutable_oftype(A, choltype(A))
@inline cholcopy(A::RealHermSymComplexHerm) = copymutable_oftype(A, choltype(A))
@inline cholcopy(A::AbstractMatrix) = copy_similar(A, choltype(A))

# _chol!. Internal methods for calling unpivoted Cholesky
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -500,7 +500,7 @@ function (^)(A::AbstractMatrix{T}, p::Real) where T
# Quicker return if A is diagonal
if isdiag(A)
TT = promote_op(^, T, typeof(p))
retmat = copy_oftype(A, TT)
retmat = copymutable_oftype(A, TT)
for i in 1:n
retmat[i, i] = retmat[i, i] ^ p
end
Expand Down
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/src/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -763,8 +763,8 @@ end
@deprecate cholesky!(A::Diagonal, ::Val{false}; check::Bool = true) cholesky!(A::Diagonal, NoPivot(); check) false
@deprecate cholesky(A::Diagonal, ::Val{false}; check::Bool = true) cholesky(A::Diagonal, NoPivot(); check) false

@inline cholcopy(A::Diagonal) = copy_oftype(A, choltype(A))
@inline cholcopy(A::RealHermSymComplexHerm{<:Real,<:Diagonal}) = copy_oftype(A, choltype(A))
@inline cholcopy(A::Diagonal) = copymutable_oftype(A, choltype(A))
@inline cholcopy(A::RealHermSymComplexHerm{<:Real,<:Diagonal}) = copymutable_oftype(A, choltype(A))

function getproperty(C::Cholesky{<:Any,<:Diagonal}, d::Symbol)
Cfactors = getfield(C, :factors)
Expand Down
10 changes: 5 additions & 5 deletions stdlib/LinearAlgebra/src/eigen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -233,12 +233,12 @@ true
```
"""
function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
AA = copy_oftype(A, eigtype(T))
AA = copymutable_oftype(A, eigtype(T))
isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
return eigen!(AA; permute=permute, scale=scale, sortby=sortby)
end
function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where {T <: Union{Float16,Complex{Float16}}}
AA = copy_oftype(A, eigtype(T))
AA = copymutable_oftype(A, eigtype(T))
isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
A = eigen!(AA; permute, scale, sortby)
values = convert(AbstractVector{isreal(A.values) ? Float16 : Complex{Float16}}, A.values)
Expand Down Expand Up @@ -333,7 +333,7 @@ julia> eigvals(diag_matrix)
```
"""
eigvals(A::AbstractMatrix{T}; kws...) where T =
eigvals!(copy_oftype(A, eigtype(T)); kws...)
eigvals!(copymutable_oftype(A, eigtype(T)); kws...)

"""
For a scalar input, `eigvals` will return a scalar.
Expand Down Expand Up @@ -508,7 +508,7 @@ true
"""
function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
S = promote_type(eigtype(TA),TB)
eigen!(copy_oftype(A, S), copy_oftype(B, S); kws...)
eigen!(copymutable_oftype(A, S), copymutable_oftype(B, S); kws...)
end

eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1))
Expand Down Expand Up @@ -587,7 +587,7 @@ julia> eigvals(A,B)
"""
function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
S = promote_type(eigtype(TA),TB)
return eigvals!(copy_oftype(A, S), copy_oftype(B, S); kws...)
return eigvals!(copymutable_oftype(A, S), copymutable_oftype(B, S); kws...)
end

"""
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1780,7 +1780,7 @@ julia> normalize(a)
function normalize(a::AbstractArray, p::Real = 2)
nrm = norm(a, p)
if !isempty(a)
aa = copy_oftype(a, typeof(first(a)/nrm))
aa = copymutable_oftype(a, typeof(first(a)/nrm))
return __normalize!(aa, nrm)
else
T = typeof(zero(eltype(a))/nrm)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/givens.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ transpose(R::AbstractRotation) = error("transpose not implemented for $(typeof(R

function (*)(R::AbstractRotation{T}, A::AbstractVecOrMat{S}) where {T,S}
TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
lmul!(convert(AbstractRotation{TS}, R), copy_oftype(A, TS))
lmul!(convert(AbstractRotation{TS}, R), copy_similar(A, TS))
end
(*)(A::AbstractVector, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
(*)(A::AbstractMatrix, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/hessenberg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -502,7 +502,7 @@ true
```
"""
hessenberg(A::AbstractMatrix{T}) where T =
hessenberg!(copy_oftype(A, eigtype(T)))
hessenberg!(copymutable_oftype(A, eigtype(T)))

function show(io::IO, mime::MIME"text/plain", F::Hessenberg)
summary(io, F)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/ldlt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ julia> S \\ b
"""
function ldlt(M::SymTridiagonal{T}; shift::Number=false) where T
S = typeof((zero(T)+shift)/one(T))
Mₛ = SymTridiagonal{S}(copy_oftype(M.dv, S), copy_oftype(M.ev, S))
Mₛ = SymTridiagonal{S}(copymutable_oftype(M.dv, S), copymutable_oftype(M.ev, S))
if !iszero(shift)
Mₛ.dv .+= shift
end
Expand Down
12 changes: 6 additions & 6 deletions stdlib/LinearAlgebra/src/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ julia> l == S.L && q == S.Q
true
```
"""
lq(A::AbstractMatrix{T}) where {T} = lq!(copy_oftype(A, lq_eltype(T)))
lq(A::AbstractMatrix{T}) where {T} = lq!(copymutable_oftype(A, lq_eltype(T)))
lq(x::Number) = lq!(fill(convert(lq_eltype(typeof(x)), x), 1, 1))

lq_eltype(::Type{T}) where {T} = typeof(zero(T) / sqrt(abs2(one(T))))
Expand Down Expand Up @@ -197,15 +197,15 @@ function lmul!(A::LQ, B::StridedVecOrMat)
end
function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
TAB = promote_type(TA, TB)
_cut_B(lmul!(convert(Factorization{TAB}, A), copy_oftype(B, TAB)), 1:size(A,1))
_cut_B(lmul!(convert(Factorization{TAB}, A), copymutable_oftype(B, TAB)), 1:size(A,1))
end

## Multiplication by Q
### QB
lmul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B)
function (*)(A::LQPackedQ, B::StridedVecOrMat)
TAB = promote_type(eltype(A), eltype(B))
lmul!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB))
lmul!(AbstractMatrix{TAB}(A), copymutable_oftype(B, TAB))
end

### QcB
Expand All @@ -218,7 +218,7 @@ function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat)
A = adjA.parent
TAB = promote_type(eltype(A), eltype(B))
if size(B,1) == size(A.factors,2)
lmul!(adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB))
lmul!(adjoint(AbstractMatrix{TAB}(A)), copymutable_oftype(B, TAB))
elseif size(B,1) == size(A.factors,1)
lmul!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))])
else
Expand Down Expand Up @@ -269,7 +269,7 @@ rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasCo
function *(A::StridedVecOrMat, adjQ::Adjoint{<:Any,<:LQPackedQ})
Q = adjQ.parent
TR = promote_type(eltype(A), eltype(Q))
return rmul!(copy_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q)))
return rmul!(copymutable_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q)))
end
function *(adjA::Adjoint{<:Any,<:StridedMatrix}, adjQ::Adjoint{<:Any,<:LQPackedQ})
A, Q = adjA.parent, adjQ.parent
Expand All @@ -293,7 +293,7 @@ end
function *(A::StridedVecOrMat, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
if size(A, 2) == size(Q.factors, 2)
C = copy_oftype(A, TR)
C = copymutable_oftype(A, TR)
elseif size(A, 2) == size(Q.factors, 1)
C = zeros(TR, size(A, 1), size(Q.factors, 2))
copyto!(C, 1, A, 1, length(A))
Expand Down
10 changes: 5 additions & 5 deletions stdlib/LinearAlgebra/src/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -283,8 +283,8 @@ end
@deprecate lu(A::AbstractMatrix, ::Val{false}; check::Bool = true) lu(A, NoPivot(); check=check)

_lucopy(A::AbstractMatrix, T) = copy_similar(A, T)
_lucopy(A::HermOrSym, T) = copy_oftype(A, T)
_lucopy(A::Tridiagonal, T) = copy_oftype(A, T)
_lucopy(A::HermOrSym, T) = copymutable_oftype(A, T)
_lucopy(A::Tridiagonal, T) = copymutable_oftype(A, T)

lu(S::LU) = S
function lu(x::Number; check::Bool=true)
Expand Down Expand Up @@ -438,18 +438,18 @@ end

function (/)(A::AbstractMatrix, F::Adjoint{<:Any,<:LU})
T = promote_type(eltype(A), eltype(F))
return adjoint(ldiv!(F.parent, copy_oftype(adjoint(A), T)))
return adjoint(ldiv!(F.parent, copymutable_oftype(adjoint(A), T)))
end
# To avoid ambiguities with definitions in adjtrans.jl and factorizations.jl
(/)(adjA::Adjoint{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent)
(/)(adjA::Adjoint{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent)
function (/)(trA::Transpose{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU})
T = promote_type(eltype(trA), eltype(F))
return adjoint(ldiv!(F.parent, conj!(copy_oftype(trA.parent, T))))
return adjoint(ldiv!(F.parent, conj!(copymutable_oftype(trA.parent, T))))
end
function (/)(trA::Transpose{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU})
T = promote_type(eltype(trA), eltype(F))
return adjoint(ldiv!(F.parent, conj!(copy_oftype(trA.parent, T))))
return adjoint(ldiv!(F.parent, conj!(copymutable_oftype(trA.parent, T))))
end

function det(F::LU{T}) where T
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@ end
function (*)(A::AdjOrTransStridedMat{<:BlasComplex}, B::StridedMaybeAdjOrTransMat{<:BlasReal})
TS = promote_type(eltype(A), eltype(B))
mul!(similar(B, TS, (size(A, 1), size(B, 2))),
copy_oftype(A, TS), # remove AdjOrTrans to use reinterpret trick below
copymutable_oftype(A, TS), # remove AdjOrTrans to use reinterpret trick below
wrapperop(B)(convert(AbstractArray{real(TS)}, _parent(B))))
end
# the following case doesn't seem to benefit from the translation A*B = (B' * A')'
Expand Down
8 changes: 4 additions & 4 deletions stdlib/LinearAlgebra/src/qr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -656,7 +656,7 @@ function (*)(A::AbstractQ, b::StridedVector)
TAb = promote_type(eltype(A), eltype(b))
Anew = convert(AbstractMatrix{TAb}, A)
if size(A.factors, 1) == length(b)
bnew = copy_oftype(b, TAb)
bnew = copymutable_oftype(b, TAb)
elseif size(A.factors, 2) == length(b)
bnew = [b; zeros(TAb, size(A.factors, 1) - length(b))]
else
Expand All @@ -668,7 +668,7 @@ function (*)(A::AbstractQ, B::StridedMatrix)
TAB = promote_type(eltype(A), eltype(B))
Anew = convert(AbstractMatrix{TAB}, A)
if size(A.factors, 1) == size(B, 1)
Bnew = copy_oftype(B, TAB)
Bnew = copymutable_oftype(B, TAB)
elseif size(A.factors, 2) == size(B, 1)
Bnew = [B; zeros(TAB, size(A.factors, 1) - size(B,1), size(B, 2))]
else
Expand Down Expand Up @@ -722,7 +722,7 @@ end
function *(adjQ::Adjoint{<:Any,<:AbstractQ}, B::StridedVecOrMat)
Q = adjQ.parent
TQB = promote_type(eltype(Q), eltype(B))
return lmul!(adjoint(convert(AbstractMatrix{TQB}, Q)), copy_oftype(B, TQB))
return lmul!(adjoint(convert(AbstractMatrix{TQB}, Q)), copymutable_oftype(B, TQB))
end

### QBc/QcBc
Expand Down Expand Up @@ -774,7 +774,7 @@ end
function (*)(A::StridedMatrix, Q::AbstractQ)
TAQ = promote_type(eltype(A), eltype(Q))

return rmul!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
return rmul!(copymutable_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
end

function (*)(a::Number, B::AbstractQ)
Expand Down
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/src/special.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,8 +50,8 @@ Bidiagonal(A::AbstractTriangular) =
isbanded(A, -1, 0) ? Bidiagonal(diag(A, 0), diag(A, -1), :L) : # is lower bidiagonal
throw(ArgumentError("matrix cannot be represented as Bidiagonal"))

_lucopy(A::Bidiagonal, T) = copy_oftype(Tridiagonal(A), T)
_lucopy(A::Diagonal, T) = copy_oftype(Tridiagonal(A), T)
_lucopy(A::Bidiagonal, T) = copymutable_oftype(Tridiagonal(A), T)
_lucopy(A::Diagonal, T) = copymutable_oftype(Tridiagonal(A), T)
function _lucopy(A::SymTridiagonal, T)
du = copy_similar(_evview(A), T)
dl = copy.(transpose.(du))
Expand Down
10 changes: 5 additions & 5 deletions stdlib/LinearAlgebra/src/svd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -176,10 +176,10 @@ true
```
"""
function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T}
svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
svd!(copymutable_oftype(A, eigtype(T)), full = full, alg = alg)
end
function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T <: Union{Float16,Complex{Float16}}}
A = svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
A = svd!(copymutable_oftype(A, eigtype(T)), full = full, alg = alg)
return SVD{T}(A)
end
function svd(x::Number; full::Bool = false, alg::Algorithm = default_svd_alg(x))
Expand Down Expand Up @@ -240,7 +240,7 @@ julia> svdvals(A)
0.0
```
"""
svdvals(A::AbstractMatrix{T}) where {T} = svdvals!(copy_oftype(A, eigtype(T)))
svdvals(A::AbstractMatrix{T}) where {T} = svdvals!(copymutable_oftype(A, eigtype(T)))
svdvals(A::AbstractVector{T}) where {T} = [convert(eigtype(T), norm(A))]
svdvals(A::AbstractMatrix{<:BlasFloat}) = svdvals!(copy(A))
svdvals(A::AbstractVector{<:BlasFloat}) = [norm(A)]
Expand Down Expand Up @@ -459,7 +459,7 @@ true
"""
function svd(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB}
S = promote_type(eigtype(TA),TB)
return svd!(copy_oftype(A, S), copy_oftype(B, S))
return svd!(copymutable_oftype(A, S), copymutable_oftype(B, S))
end
# This method can be heavily optimized but it is probably not critical
# and might introduce bugs or inconsistencies relative to the 1x1 matrix
Expand Down Expand Up @@ -569,7 +569,7 @@ julia> svdvals(A, B)
"""
function svdvals(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB}
S = promote_type(eigtype(TA), TB)
return svdvals!(copy_oftype(A, S), copy_oftype(B, S))
return svdvals!(copymutable_oftype(A, S), copymutable_oftype(B, S))
end
svdvals(x::Number, y::Number) = abs(x/y)

Expand Down
12 changes: 6 additions & 6 deletions stdlib/LinearAlgebra/src/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -281,30 +281,30 @@ ldiv!(A::SymTridiagonal, B::AbstractVecOrMat; shift::Number=false) = ldiv!(ldlt(
rdiv!(B::AbstractVecOrMat, A::SymTridiagonal; shift::Number=false) = rdiv!(B, ldlt(A, shift=shift))

eigen!(A::SymTridiagonal{<:BlasReal}) = Eigen(LAPACK.stegr!('V', A.dv, A.ev)...)
eigen(A::SymTridiagonal{T}) where T = eigen!(copy_oftype(A, eigtype(T)))
eigen(A::SymTridiagonal{T}) where T = eigen!(copymutable_oftype(A, eigtype(T)))

eigen!(A::SymTridiagonal{<:BlasReal}, irange::UnitRange) =
Eigen(LAPACK.stegr!('V', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)...)
eigen(A::SymTridiagonal{T}, irange::UnitRange) where T =
eigen!(copy_oftype(A, eigtype(T)), irange)
eigen!(copymutable_oftype(A, eigtype(T)), irange)

eigen!(A::SymTridiagonal{<:BlasReal}, vl::Real, vu::Real) =
Eigen(LAPACK.stegr!('V', 'V', A.dv, A.ev, vl, vu, 0, 0)...)
eigen(A::SymTridiagonal{T}, vl::Real, vu::Real) where T =
eigen!(copy_oftype(A, eigtype(T)), vl, vu)
eigen!(copymutable_oftype(A, eigtype(T)), vl, vu)

eigvals!(A::SymTridiagonal{<:BlasReal}) = LAPACK.stev!('N', A.dv, A.ev)[1]
eigvals(A::SymTridiagonal{T}) where T = eigvals!(copy_oftype(A, eigtype(T)))
eigvals(A::SymTridiagonal{T}) where T = eigvals!(copymutable_oftype(A, eigtype(T)))

eigvals!(A::SymTridiagonal{<:BlasReal}, irange::UnitRange) =
LAPACK.stegr!('N', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)[1]
eigvals(A::SymTridiagonal{T}, irange::UnitRange) where T =
eigvals!(copy_oftype(A, eigtype(T)), irange)
eigvals!(copymutable_oftype(A, eigtype(T)), irange)

eigvals!(A::SymTridiagonal{<:BlasReal}, vl::Real, vu::Real) =
LAPACK.stegr!('N', 'V', A.dv, A.ev, vl, vu, 0, 0)[1]
eigvals(A::SymTridiagonal{T}, vl::Real, vu::Real) where T =
eigvals!(copy_oftype(A, eigtype(T)), vl, vu)
eigvals!(copymutable_oftype(A, eigtype(T)), vl, vu)

#Computes largest and smallest eigenvalue
eigmax(A::SymTridiagonal) = eigvals(A, size(A, 1):size(A, 1))[1]
Expand Down
Loading