diff --git a/Project.toml b/Project.toml index 0551a6e..8ff53c4 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,8 @@ uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a" version = "0.2.4" [deps] +GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" +GenericSVD = "01680d73-4ee2-5a08-a1aa-533608c188bb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -16,3 +18,5 @@ test = ["Quaternions", "Test"] [compat] julia = "1.3" +GenericSchur = "0.4" +GenericSVD = "0.3" diff --git a/src/GenericLinearAlgebra.jl b/src/GenericLinearAlgebra.jl index 2ab778b..f796da0 100644 --- a/src/GenericLinearAlgebra.jl +++ b/src/GenericLinearAlgebra.jl @@ -1,6 +1,7 @@ module GenericLinearAlgebra import LinearAlgebra: mul!, ldiv! +import GenericSVD, GenericSchur include("juliaBLAS.jl") include("lapack.jl") diff --git a/src/eigenGeneral.jl b/src/eigenGeneral.jl index f1f7018..9f98c38 100644 --- a/src/eigenGeneral.jl +++ b/src/eigenGeneral.jl @@ -36,29 +36,6 @@ function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat) ldiv!(Triangular(Hd, :U), B) end -# Hessenberg factorization -struct HessenbergFactorization{T, S<:StridedMatrix,U} <: Factorization{T} - data::S - τ::Vector{U} -end - -function _hessenberg!(A::StridedMatrix{T}) where T - n = LinearAlgebra.checksquare(A) - τ = Vector{Householder{T}}(undef, n - 1) - for i = 1:n - 1 - xi = view(A, i + 1:n, i) - t = LinearAlgebra.reflector!(xi) - H = Householder{T,typeof(xi)}(view(xi, 2:n - i), t) - τ[i] = H - lmul!(H', view(A, i + 1:n, i + 1:n)) - rmul!(view(A, :, i + 1:n), H) - end - return HessenbergFactorization{T, typeof(A), eltype(τ)}(A, τ) -end -LinearAlgebra.hessenberg!(A::StridedMatrix) = _hessenberg!(A) - -Base.size(H::HessenbergFactorization, args...) = size(H.data, args...) - # Schur struct Schur{T,S<:StridedMatrix} <: Factorization{T} data::S @@ -81,7 +58,7 @@ end # We currently absorb extra unsupported keywords in kwargs. These could e.g. be scale and permute. Do we want to check that these are false? function _schur!( - H::HessenbergFactorization{T}; + H::GenericSchur.HessenbergArg{T}; tol = eps(real(T)), debug = false, shiftmethod = :Francis, @@ -91,7 +68,7 @@ function _schur!( n = size(H, 1) istart = 1 iend = n - HH = H.data + HH = GenericSchur._getdata(H) τ = Rotation(Givens{T}[]) # iteration count @@ -166,7 +143,7 @@ function _schur!( return Schur{T,typeof(HH)}(HH, τ) end -_schur!(A::StridedMatrix; kwargs...) = _schur!(_hessenberg!(A); kwargs...) +_schur!(A::StridedMatrix; kwargs...) = _schur!(hessenberg!(A); kwargs...) LinearAlgebra.schur!(A::StridedMatrix; kwargs...) = _schur!(A; kwargs...) function singleShiftQR!(HH::StridedMatrix, τ::Rotation, shift::Number, istart::Integer, iend::Integer) @@ -245,7 +222,7 @@ end _eigvals!(A::StridedMatrix; kwargs...) = _eigvals!(_schur!(A; kwargs...)) _eigvals!(H::HessenbergMatrix; kwargs...) = _eigvals!(_schur!(H; kwargs...)) -_eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(_schur!(H; kwargs...)) +_eigvals!(H::GenericSchur.HessenbergArg; kwargs...) = _eigvals!(_schur!(H; kwargs...)) # Overload methods from LinearAlgebra to make them work generically if VERSION > v"1.2.0-DEV.0" @@ -266,7 +243,7 @@ if VERSION > v"1.2.0-DEV.0" kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby) LinearAlgebra.eigvals!( - H::HessenbergFactorization; + H::GenericSchur.HessenbergArg; sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby) else @@ -280,7 +257,7 @@ else LinearAlgebra.eigvals!(H::HessenbergMatrix; kwargs...) = _eigvals!(H; kwargs...) - LinearAlgebra.eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(H; kwargs...) + LinearAlgebra.eigvals!(H::GenericSchur.HessenbergArg; kwargs...) = _eigvals!(H; kwargs...) end # To compute the eigenvalue of the pseudo triangular Schur matrix we just return diff --git a/src/qr.jl b/src/qr.jl index b26df03..d2fb3d2 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -13,28 +13,6 @@ QR2(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = QR2{T,typeof(factors) size(F::QR2, i::Integer...) = size(F.factors, i...) -# Similar to the definition in base but applies the reflector from the right -@inline function reflectorApply!(A::StridedMatrix, x::AbstractVector, τ::Number) # apply conjugate transpose reflector from right. - m, n = size(A) - if length(x) != n - throw(DimensionMismatch("reflector must have same length as second dimension of matrix")) - end - @inbounds begin - for i in 1:m - Aiv = A[i, 1] - for j in 2:n - Aiv += A[i, j]*x[j] - end - Aiv = Aiv*τ - A[i, 1] -= Aiv - for j in 2:n - A[i, j] -= Aiv*x[j]' - end - end - end - return A -end - # FixMe! Consider how to represent Q # immutable Q{T,S<:QR2} <: AbstractMatrix{T} diff --git a/src/svd.jl b/src/svd.jl index 68944df..e6d206c 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -2,9 +2,6 @@ using LinearAlgebra import LinearAlgebra: mul!, rmul! -lmul!(G::LinearAlgebra.Givens, ::Nothing) = nothing -rmul!(::Nothing, G::LinearAlgebra.Givens) = nothing - function svdvals2x2(d1, d2, e) d1sq = d1*d1 d2sq = d2*d2