From 10b5265b75e123f236ee50b03c845c642462c30e Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 7 Dec 2022 16:19:53 +0100 Subject: [PATCH 1/4] Bump test dep Quaternions.jl (#99) --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 79312ad..a2b58e0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "GenericLinearAlgebra" uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a" -version = "0.3.5" +version = "0.3.6-DEV" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -9,7 +9,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" libblastrampoline_jll = "8e850b90-86db-534c-a0d3-1478176c7d93" [compat] -Quaternions = "=0.6.0" +Quaternions = "0.7.0" julia = "1.6" [extras] From dd434ca3fe14ca073ab33bdf71393190810dac51 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Fri, 24 Feb 2023 11:40:44 +0100 Subject: [PATCH 2/4] Fix size of Householder (#101) Also improve error messages and extend Base's convert instead of defining a new one --- src/householder.jl | 14 +++++++------- test/eigengeneral.jl | 12 ++++++++++++ 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/householder.jl b/src/householder.jl index 65cfd51..34913c4 100644 --- a/src/householder.jl +++ b/src/householder.jl @@ -10,8 +10,8 @@ struct HouseholderBlock{T,S<:StridedMatrix,U<:StridedMatrix} T::UpperTriangular{T,U} end -size(H::Householder) = (length(H.v), length(H.v)) -size(H::Householder, i::Integer) = i <= 2 ? length(H.v) : 1 +size(H::Householder) = (length(H.v) + 1, length(H.v) + 1) +size(H::Householder, i::Integer) = i <= 2 ? length(H.v) + 1 : 1 eltype(H::Householder{T}) where T = T eltype(H::HouseholderBlock{T}) where T = T @@ -21,7 +21,7 @@ adjoint(H::HouseholderBlock{T}) where {T} = Adjoint{T,typeof(H)}(H) function lmul!(H::Householder, A::StridedMatrix) m, n = size(A) - length(H.v) == m - 1 || throw(DimensionMismatch("")) + length(H.v) == m - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))")) v = view(H.v, 1:m - 1) τ = H.τ for j = 1:n @@ -37,7 +37,7 @@ end function rmul!(A::StridedMatrix, H::Householder) m, n = size(A) - length(H.v) == n - 1 || throw(DimensionMismatch("")) + length(H.v) == n - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but second dimension of matrix is $(size(A, 2))")) v = view(H.v, :) τ = H.τ a1 = view(A, :, 1) @@ -52,7 +52,7 @@ end function lmul!(adjH::Adjoint{<:Any,<:Householder}, A::StridedMatrix) H = parent(adjH) m, n = size(A) - length(H.v) == m - 1 || throw(DimensionMismatch("")) + length(H.v) == m - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))")) v = view(H.v, 1:m - 1) τ = H.τ for j = 1:n @@ -142,5 +142,5 @@ end (*)(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}) where {T} = lmul!(adjH, copy(A), similar(A, (min(size(parent(adjH).V)...), size(A, 2)))) -convert(::Type{Matrix}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1))) -convert(::Type{Matrix{T}}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1))) +Base.convert(::Type{Matrix}, H::Householder{T}) where {T} = convert(Matrix{T}, H) +Base.convert(::Type{Matrix{T}}, H::Householder) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1))) diff --git a/test/eigengeneral.jl b/test/eigengeneral.jl index 5668029..12e4ab6 100644 --- a/test/eigengeneral.jl +++ b/test/eigengeneral.jl @@ -144,4 +144,16 @@ end @test sort(imag(vals)) ≈ sort(imag(λs)) atol=1e-25 end +@testset "_hessenberg! and Hessenberg" begin + n = 10 + A = randn(n, n) + HF = GenericLinearAlgebra._hessenberg!(copy(A)) + for i in 1:length(HF.τ) + HM = convert(Matrix, HF.τ[i]) + A[(i + 1):end, :] = HM * A[(i + 1):end, :] + A[:, (i + 1):end] = A[:, (i + 1):end] * HM' + end + @test tril(A, -2) ≈ zeros(n, n) atol = 1e-14 +end + end \ No newline at end of file From 06be73aa0024f8217ee4465f26d52ff56eaee0e1 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Fri, 24 Feb 2023 13:24:39 +0100 Subject: [PATCH 3/4] Format the source with JuliaFormatter (#102) --- src/cholesky.jl | 24 +- src/eigenGeneral.jl | 213 +++--- src/eigenSelfAdjoint.jl | 546 ++++++++------- src/householder.jl | 55 +- src/juliaBLAS.jl | 134 ++-- src/lapack.jl | 1375 +++++++++++++++++++++++++------------- src/qr.jl | 72 +- src/rectfullpacked.jl | 71 +- src/svd.jl | 317 +++++---- src/tridiag.jl | 14 +- test/cholesky.jl | 21 +- test/eigengeneral.jl | 343 ++++++---- test/eigenselfadjoint.jl | 52 +- test/juliaBLAS.jl | 16 +- test/lapack.jl | 4 +- test/qr.jl | 25 +- test/rectfullpacked.jl | 61 +- test/runtests.jl | 18 +- test/svd.jl | 90 +-- test/tridiag.jl | 5 +- 20 files changed, 2136 insertions(+), 1320 deletions(-) diff --git a/src/cholesky.jl b/src/cholesky.jl index 6cb85f4..9029fcf 100644 --- a/src/cholesky.jl +++ b/src/cholesky.jl @@ -1,11 +1,11 @@ using LinearAlgebra: rdiv! -function cholUnblocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}) where T<:Number +function cholUnblocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}) where {T<:Number} n = LinearAlgebra.checksquare(A) - A[1,1] = sqrt(A[1,1]) + A[1, 1] = sqrt(A[1, 1]) if n > 1 a21 = view(A, 2:n, 1) - rmul!(a21, inv(real(A[1,1]))) + rmul!(a21, inv(real(A[1, 1]))) A22 = view(A, 2:n, 2:n) rankUpdate!(Hermitian(A22, :L), a21, -1) @@ -14,36 +14,40 @@ function cholUnblocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}) where T<:Number A end -function cholBlocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}, blocksize::Integer) where T<:Number +function cholBlocked!( + A::AbstractMatrix{T}, + ::Type{Val{:L}}, + blocksize::Integer, +) where {T<:Number} n = LinearAlgebra.checksquare(A) mnb = min(n, blocksize) A11 = view(A, 1:mnb, 1:mnb) cholUnblocked!(A11, Val{:L}) if n > blocksize - A21 = view(A, (blocksize + 1):n, 1:blocksize) + A21 = view(A, (blocksize+1):n, 1:blocksize) rdiv!(A21, LowerTriangular(A11)') - A22 = view(A, (blocksize + 1):n, (blocksize + 1):n) + A22 = view(A, (blocksize+1):n, (blocksize+1):n) rankUpdate!(Hermitian(A22, :L), A21, -1) cholBlocked!(A22, Val{:L}, blocksize) end A end -function cholRecursive!(A::StridedMatrix{T}, ::Type{Val{:L}}, cutoff = 1) where T +function cholRecursive!(A::StridedMatrix{T}, ::Type{Val{:L}}, cutoff = 1) where {T} n = LinearAlgebra.checksquare(A) if n == 1 - A[1,1] = sqrt(A[1,1]) + A[1, 1] = sqrt(A[1, 1]) elseif n < cutoff cholUnblocked!(A, Val{:L}) else n2 = div(n, 2) A11 = view(A, 1:n2, 1:n2) cholRecursive!(A11, Val{:L}) - A21 = view(A, n2 + 1:n, 1:n2) + A21 = view(A, n2+1:n, 1:n2) rdiv!(A21, LowerTriangular(A11)') - A22 = view(A, n2 + 1:n, n2 + 1:n) + A22 = view(A, n2+1:n, n2+1:n) rankUpdate!(Hermitian(A22, :L), A21, -1) cholRecursive!(A22, Val{:L}, cutoff) end diff --git a/src/eigenGeneral.jl b/src/eigenGeneral.jl index ab8f71e..d6edf3f 100644 --- a/src/eigenGeneral.jl +++ b/src/eigenGeneral.jl @@ -6,8 +6,8 @@ using LinearAlgebra: Givens, Rotation function adiagmax(A::StridedMatrix) adm = zero(typeof(real(A[1]))) @inbounds begin - for i = size(A,1) - adm = max(adm, abs(A[i,i])) + for i in size(A, 1) + adm = max(adm, abs(A[i, i])) end end return adm @@ -20,7 +20,8 @@ end Base.copy(H::HessenbergMatrix{T,S}) where {T,S} = HessenbergMatrix{T,S}(copy(H.data)) -Base.getindex(H::HessenbergMatrix{T,S}, i::Integer, j::Integer) where {T,S} = i > j + 1 ? zero(T) : H.data[i,j] +Base.getindex(H::HessenbergMatrix{T,S}, i::Integer, j::Integer) where {T,S} = + i > j + 1 ? zero(T) : H.data[i, j] Base.size(H::HessenbergMatrix) = size(H.data) Base.size(H::HessenbergMatrix, i::Integer) = size(H.data, i) @@ -29,7 +30,7 @@ function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat) n = size(H, 1) Hd = H.data for i = 1:n-1 - G, _ = givens!(Hd, i, i+1, i) + G, _ = givens!(Hd, i, i + 1, i) lmul!(G, view(Hd, 1:n, i+1:n)) lmul!(G, B) end @@ -37,23 +38,23 @@ function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat) end # Hessenberg factorization -struct HessenbergFactorization{T, S<:StridedMatrix,U} <: Factorization{T} +struct HessenbergFactorization{T,S<:StridedMatrix,U} <: Factorization{T} data::S τ::Vector{U} end -function _hessenberg!(A::StridedMatrix{T}) where T +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) + 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) + 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, τ) + return HessenbergFactorization{T,typeof(A),eltype(τ)}(A, τ) end LinearAlgebra.hessenberg!(A::StridedMatrix) = _hessenberg!(A) @@ -74,8 +75,8 @@ function Base.getproperty(F::Schur, s::Symbol) end function wilkinson(Hmm, t, d) - λ1 = (t + sqrt(t*t - 4d))/2 - λ2 = (t - sqrt(t*t - 4d))/2 + λ1 = (t + sqrt(t * t - 4d)) / 2 + λ2 = (t - sqrt(t * t - 4d)) / 2 return ifelse(abs(Hmm - λ1) < abs(Hmm - λ2), λ1, λ2) end @@ -85,8 +86,9 @@ function _schur!( tol = eps(real(T)), debug = false, shiftmethod = :Francis, - maxiter = 30*size(H, 1), - kwargs...) where T + maxiter = 30 * size(H, 1), + kwargs..., +) where {T} n = size(H, 1) istart = 1 @@ -105,17 +107,22 @@ function _schur!( end # Determine if the matrix splits. Find lowest positioned subdiagonal "zero" - for _istart in iend - 1:-1:1 - if abs(HH[_istart + 1, _istart]) <= tol*(abs(HH[_istart, _istart]) + abs(HH[_istart + 1, _istart + 1])) + for _istart = iend-1:-1:1 + if abs(HH[_istart+1, _istart]) <= + tol * (abs(HH[_istart, _istart]) + abs(HH[_istart+1, _istart+1])) # Check if subdiagonal element H[i+1,i] is "zero" such that we can split the matrix - + # Set istart to the beginning of the second block istart = _istart + 1 - debug && @printf("Split! Subdiagonal element is: %10.3e and istart now %6d\n", HH[istart, istart - 1], istart) + debug && @printf( + "Split! Subdiagonal element is: %10.3e and istart now %6d\n", + HH[istart, istart-1], + istart + ) # Set the subdiagonal element to zero to signal that a split has taken place - HH[istart, istart - 1] = 0 + HH[istart, istart-1] = 0 break end @@ -125,43 +132,64 @@ function _schur!( # if block size is one we deflate if istart >= iend - debug && @printf("Bottom deflation! Block size is one. New iend is %6d\n", iend - 1) + debug && + @printf("Bottom deflation! Block size is one. New iend is %6d\n", iend - 1) iend -= 1 - # and the same for a 2x2 block + # and the same for a 2x2 block elseif istart + 1 == iend - debug && @printf("Bottom deflation! Block size is two. New iend is %6d\n", iend - 2) + debug && + @printf("Bottom deflation! Block size is two. New iend is %6d\n", iend - 2) iend -= 2 - # run a QR iteration - # shift method is specified with shiftmethod kw argument + # run a QR iteration + # shift method is specified with shiftmethod kw argument else Hmm = HH[iend, iend] - Hm1m1 = HH[iend - 1, iend - 1] + Hm1m1 = HH[iend-1, iend-1] if iszero(i % 10) # Use Eispack exceptional shift - β = abs(HH[iend, iend - 1]) + abs(HH[iend - 1, iend - 2]) - d = (Hmm + β)^2 - Hmm*β/2 - t = 2*Hmm + 3*β/2 + β = abs(HH[iend, iend-1]) + abs(HH[iend-1, iend-2]) + d = (Hmm + β)^2 - Hmm * β / 2 + t = 2 * Hmm + 3 * β / 2 else - d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend] + d = Hm1m1 * Hmm - HH[iend, iend-1] * HH[iend-1, iend] t = Hm1m1 + Hmm end - debug && @printf("block start is: %6d, block end is: %6d, d: %10.3e, t: %10.3e\n", istart, iend, d, t) + debug && @printf( + "block start is: %6d, block end is: %6d, d: %10.3e, t: %10.3e\n", + istart, + iend, + d, + t + ) if shiftmethod == :Francis - debug && @printf("Francis double shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2]) + debug && @printf( + "Francis double shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", + HH[iend, iend-1], + HH[iend-1, iend-2] + ) doubleShiftQR!(HH, τ, t, d, istart, iend) elseif shiftmethod == :Rayleigh - debug && @printf("Single shift with Rayleigh shift! Subdiagonal is: %10.3e\n", HH[iend, iend - 1]) + debug && @printf( + "Single shift with Rayleigh shift! Subdiagonal is: %10.3e\n", + HH[iend, iend-1] + ) # Run a bulge chase singleShiftQR!(HH, τ, Hmm, istart, iend) else - throw(ArgumentError("only support supported shift methods are :Francis (default) and :Rayleigh. You supplied $shiftmethod")) + throw( + ArgumentError( + "only support supported shift methods are :Francis (default) and :Rayleigh. You supplied $shiftmethod", + ), + ) end end - if iend <= 2 break end + if iend <= 2 + break + end end return Schur{T,typeof(HH)}(HH, τ) @@ -169,25 +197,31 @@ end _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) +function singleShiftQR!( + HH::StridedMatrix, + τ::Rotation, + shift::Number, + istart::Integer, + iend::Integer, +) m = size(HH, 1) H11 = HH[istart, istart] - H21 = HH[istart + 1, istart] + H21 = HH[istart+1, istart] if m > istart + 1 - Htmp = HH[istart + 2, istart] - HH[istart + 2, istart] = 0 + Htmp = HH[istart+2, istart] + HH[istart+2, istart] = 0 end G, _ = givens(H11 - shift, H21, istart, istart + 1) lmul!(G, view(HH, :, istart:m)) rmul!(view(HH, 1:min(istart + 2, iend), :), G') lmul!(G, τ) - for i = istart:iend - 2 - G, _ = givens(HH[i + 1, i], HH[i + 2, i], i + 1, i + 2) + for i = istart:iend-2 + G, _ = givens(HH[i+1, i], HH[i+2, i], i + 1, i + 2) lmul!(G, view(HH, :, i:m)) - HH[i + 2, i] = Htmp + HH[i+2, i] = Htmp if i < iend - 2 - Htmp = HH[i + 3, i + 1] - HH[i + 3, i + 1] = 0 + Htmp = HH[i+3, i+1] + HH[i+3, i+1] = 0 end rmul!(view(HH, 1:min(i + 3, iend), :), G') # mul!(G, τ) @@ -195,23 +229,35 @@ function singleShiftQR!(HH::StridedMatrix, τ::Rotation, shift::Number, istart:: return HH end -function doubleShiftQR!(HH::StridedMatrix, τ::Rotation, shiftTrace::Number, shiftDeterminant::Number, istart::Integer, iend::Integer) +function doubleShiftQR!( + HH::StridedMatrix, + τ::Rotation, + shiftTrace::Number, + shiftDeterminant::Number, + istart::Integer, + iend::Integer, +) m = size(HH, 1) H11 = HH[istart, istart] - H21 = HH[istart + 1, istart] - Htmp11 = HH[istart + 2, istart] - HH[istart + 2, istart] = 0 + H21 = HH[istart+1, istart] + Htmp11 = HH[istart+2, istart] + HH[istart+2, istart] = 0 if istart + 3 <= m - Htmp21 = HH[istart + 3, istart] - HH[istart + 3, istart] = 0 - Htmp22 = HH[istart + 3, istart + 1] - HH[istart + 3, istart + 1] = 0 + Htmp21 = HH[istart+3, istart] + HH[istart+3, istart] = 0 + Htmp22 = HH[istart+3, istart+1] + HH[istart+3, istart+1] = 0 else # values doen't matter in this case but variables should be initialized Htmp21 = Htmp22 = Htmp11 end - G1, r = givens(H11*H11 + HH[istart, istart + 1]*H21 - shiftTrace*H11 + shiftDeterminant, H21*(H11 + HH[istart + 1, istart + 1] - shiftTrace), istart, istart + 1) - G2, _ = givens(r, H21*HH[istart + 2, istart + 1], istart, istart + 2) + G1, r = givens( + H11 * H11 + HH[istart, istart+1] * H21 - shiftTrace * H11 + shiftDeterminant, + H21 * (H11 + HH[istart+1, istart+1] - shiftTrace), + istart, + istart + 1, + ) + G2, _ = givens(r, H21 * HH[istart+2, istart+1], istart, istart + 2) vHH = view(HH, :, istart:m) lmul!(G1, vHH) lmul!(G2, vHH) @@ -220,21 +266,23 @@ function doubleShiftQR!(HH::StridedMatrix, τ::Rotation, shiftTrace::Number, shi rmul!(vHH, G2') lmul!(G1, τ) lmul!(G2, τ) - for i = istart:iend - 2 + for i = istart:iend-2 for j = 1:2 - if i + j + 1 > iend break end + if i + j + 1 > iend + break + end # G, _ = givens(H.H,i+1,i+j+1,i) - G, _ = givens(HH[i + 1, i], HH[i + j + 1, i], i + 1, i + j + 1) + G, _ = givens(HH[i+1, i], HH[i+j+1, i], i + 1, i + j + 1) lmul!(G, view(HH, :, i:m)) - HH[i + j + 1, i] = Htmp11 + HH[i+j+1, i] = Htmp11 Htmp11 = Htmp21 # if i + j + 2 <= iend - # Htmp21 = HH[i + j + 2, i + 1] - # HH[i + j + 2, i + 1] = 0 + # Htmp21 = HH[i + j + 2, i + 1] + # HH[i + j + 2, i + 1] = 0 # end if i + 4 <= iend - Htmp22 = HH[i + 4, i + j] - HH[i + 4, i + j] = 0 + Htmp22 = HH[i+4, i+j] + HH[i+4, i+j] = 0 end rmul!(view(HH, 1:min(i + j + 2, iend), :), G') # mul!(G, τ) @@ -243,15 +291,16 @@ function doubleShiftQR!(HH::StridedMatrix, τ::Rotation, shiftTrace::Number, shi return HH end -_eigvals!(A::StridedMatrix; kwargs...) = _eigvals!(_schur!(A; kwargs...)) -_eigvals!(H::HessenbergMatrix; kwargs...) = _eigvals!(_schur!(H; kwargs...)) +_eigvals!(A::StridedMatrix; kwargs...) = _eigvals!(_schur!(A; kwargs...)) +_eigvals!(H::HessenbergMatrix; kwargs...) = _eigvals!(_schur!(H; kwargs...)) _eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(_schur!(H; kwargs...)) function LinearAlgebra.eigvals!( A::StridedMatrix; - sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, - kwargs...) - + sortby::Union{Function,Nothing} = LinearAlgebra.eigsortby, + kwargs..., +) + if ishermitian(A) return LinearAlgebra.sorteig!(eigvals!(Hermitian(A)), sortby) end @@ -260,38 +309,40 @@ end LinearAlgebra.eigvals!( H::HessenbergMatrix; - sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, - kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby) + sortby::Union{Function,Nothing} = LinearAlgebra.eigsortby, + kwargs..., +) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby) LinearAlgebra.eigvals!( H::HessenbergFactorization; - sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby, - kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby) + sortby::Union{Function,Nothing} = LinearAlgebra.eigsortby, + kwargs..., +) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby) # To compute the eigenvalue of the pseudo triangular Schur matrix we just return # the values of the 1x1 diagonal blocks and compute the eigenvalues of the 2x2 # blocks on the fly. We can check for exact zeros in the subdiagonal because # the QR algorithm sets the elements to zero when a split happens -function _eigvals!(S::Schur{T}; tol = eps(real(T))) where T +function _eigvals!(S::Schur{T}; tol = eps(real(T))) where {T} HH = S.data n = size(HH, 1) vals = Vector{complex(T)}(undef, n) i = 1 while i < n Hii = HH[i, i] - if iszero(HH[i + 1, i]) + if iszero(HH[i+1, i]) # 1x1 block vals[i] = Hii i += 1 else # 2x2 block - Hi1i1 = HH[i + 1, i + 1] - d = Hii*Hi1i1 - HH[i, i + 1]*HH[i + 1, i] + Hi1i1 = HH[i+1, i+1] + d = Hii * Hi1i1 - HH[i, i+1] * HH[i+1, i] t = Hii + Hi1i1 - x = 0.5*t - y = sqrt(complex(x*x - d)) - vals[i ] = x + y - vals[i + 1] = x - y + x = 0.5 * t + y = sqrt(complex(x * x - d)) + vals[i] = x + y + vals[i+1] = x - y i += 2 end end diff --git a/src/eigenSelfAdjoint.jl b/src/eigenSelfAdjoint.jl index 0b6dd3c..60627d7 100644 --- a/src/eigenSelfAdjoint.jl +++ b/src/eigenSelfAdjoint.jl @@ -42,32 +42,32 @@ function LinearAlgebra.lmul!(Q::EigenQ, B::StridedVecOrMat) throw(DimensionMismatch("")) end if Q.uplo == 'L' - for k in length(Q.τ):-1:1 - for j in 1:size(B, 2) + for k = length(Q.τ):-1:1 + for j = 1:size(B, 2) b = view(B, :, j) - vb = B[k + 1, j] - for i in (k + 2):m + vb = B[k+1, j] + for i = (k+2):m vb += Q.factors[i, k]'B[i, j] end τkvb = Q.τ[k]'vb - B[k + 1, j] -= τkvb - for i in (k + 2):m - B[i, j] -= Q.factors[i, k]*τkvb + B[k+1, j] -= τkvb + for i = (k+2):m + B[i, j] -= Q.factors[i, k] * τkvb end end end elseif Q.uplo == 'U' - for k in length(Q.τ):-1:1 - for j in 1:size(B, 2) + for k = length(Q.τ):-1:1 + for j = 1:size(B, 2) b = view(B, :, j) - vb = B[k + 1, j] - for i in (k + 2):m - vb += Q.factors[k, i]*B[i, j] + vb = B[k+1, j] + for i = (k+2):m + vb += Q.factors[k, i] * B[i, j] end τkvb = Q.τ[k]'vb - B[k + 1, j] -= τkvb - for i in (k + 2):m - B[i, j] -= Q.factors[k, i]'*τkvb + B[k+1, j] -= τkvb + for i = (k+2):m + B[i, j] -= Q.factors[k, i]' * τkvb end end end @@ -83,32 +83,32 @@ function LinearAlgebra.rmul!(A::StridedMatrix, Q::EigenQ) throw(DimensionMismatch("")) end if Q.uplo == 'L' - for k in 1:length(Q.τ) - for i in 1:size(A, 1) + for k = 1:length(Q.τ) + for i = 1:size(A, 1) a = view(A, i, :) - va = A[i, k + 1] - for j in (k + 2):n - va += A[i, j]*Q.factors[j, k] + va = A[i, k+1] + for j = (k+2):n + va += A[i, j] * Q.factors[j, k] end - τkva = va*Q.τ[k]' - A[i, k + 1] -= τkva - for j in (k + 2):n - A[i, j] -= τkva*Q.factors[j, k]' + τkva = va * Q.τ[k]' + A[i, k+1] -= τkva + for j = (k+2):n + A[i, j] -= τkva * Q.factors[j, k]' end end end elseif Q.uplo == 'U' # FixMe! Should consider reordering loops - for k in 1:length(Q.τ) - for i in 1:size(A, 1) + for k = 1:length(Q.τ) + for i = 1:size(A, 1) a = view(A, i, :) - va = A[i, k + 1] - for j in (k + 2):n - va += A[i, j]*Q.factors[k, j]' + va = A[i, k+1] + for j = (k+2):n + va += A[i, j] * Q.factors[k, j]' end - τkva = va*Q.τ[k]' - A[i, k + 1] -= τkva - for j in (k + 2):n - A[i, j] -= τkva*Q.factors[k, j] + τkva = va * Q.τ[k]' + A[i, k+1] -= τkva + for j = (k+2):n + A[i, j] -= τkva * Q.factors[k, j] end end end @@ -122,30 +122,30 @@ Base.Array(Q::EigenQ) = lmul!(Q, Matrix{eltype(Q)}(I, size(Q, 1), size(Q, 1))) function _updateVectors!(c, s, j, vectors) @inbounds for i = 1:size(vectors, 1) - v1 = vectors[i, j + 1] + v1 = vectors[i, j+1] v2 = vectors[i, j] - vectors[i, j + 1] = c*v1 + s*v2 - vectors[i, j] = c*v2 - s*v1 + vectors[i, j+1] = c * v1 + s * v2 + vectors[i, j] = c * v2 - s * v1 end end function eigvals2x2(d1::Number, d2::Number, e::Number) - r1 = (d1 + d2)/2 - r2 = hypot(d1 - d2, 2*e)/2 + r1 = (d1 + d2) / 2 + r2 = hypot(d1 - d2, 2 * e) / 2 return r1 + r2, r1 - r2 end function eigQL2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matrix) dj = d[j] - dj1 = d[j + 1] + dj1 = d[j+1] ej = e[j] - r1 = (dj + dj1)/2 - r2 = hypot(dj - dj1, 2*ej)/2 + r1 = (dj + dj1) / 2 + r2 = hypot(dj - dj1, 2 * ej) / 2 λ = r1 + r2 d[j] = λ - d[j + 1] = r1 - r2 + d[j+1] = r1 - r2 e[j] = 0 - c = ej/(dj - λ) + c = ej / (dj - λ) if isfinite(c) # FixMe! is this the right fix for overflow? h = hypot(one(c), c) c /= h @@ -160,7 +160,11 @@ function eigQL2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matr return c, s end -function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), debug::Bool=false) where T<:Real +function eigvalsPWK!( + S::SymTridiagonal{T}; + tol = eps(T), + debug::Bool = false, +) where {T<:Real} d = S.dv e = S.ev n = length(d) @@ -168,13 +172,13 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), debug::Bool=false) wher blockend = n iter = 0 @inbounds begin - for i = 1:n - 1 + for i = 1:n-1 e[i] = abs2(e[i]) end while true # Check for zero off diagonal elements - for be in blockstart + 1:n - if abs(e[be - 1]) <= tol*sqrt(abs(d[be - 1]))*sqrt(abs(d[be])) + for be = blockstart+1:n + if abs(e[be-1]) <= tol * sqrt(abs(d[be-1])) * sqrt(abs(d[be])) blockend = be - 1 break end @@ -187,20 +191,29 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), debug::Bool=false) wher blockstart += 1 elseif blockstart + 1 == blockend debug && println("Yes, but after sreolving 2x2 block") - d[blockstart], d[blockend] = eigvals2x2(d[blockstart], d[blockend], sqrt(e[blockstart])) + d[blockstart], d[blockend] = + eigvals2x2(d[blockstart], d[blockend], sqrt(e[blockstart])) e[blockstart] = zero(T) blockstart += 1 else # Calculate shift sqrte = sqrt(e[blockstart]) - μ = (d[blockstart + 1] - d[blockstart])/(2*sqrte) + μ = (d[blockstart+1] - d[blockstart]) / (2 * sqrte) r = hypot(μ, one(T)) - μ = d[blockstart] - sqrte/(μ + copysign(r, μ)) + μ = d[blockstart] - sqrte / (μ + copysign(r, μ)) # QL bulk chase singleShiftQLPWK!(S, μ, blockstart, blockend) - debug && @printf("QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %e, rotations: %d\n", blockstart, blockend, e[blockstart], e[blockend-1], μ, iter += blockend - blockstart) + debug && @printf( + "QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %e, rotations: %d\n", + blockstart, + blockend, + e[blockstart], + e[blockend-1], + μ, + iter += blockend - blockstart + ) end if blockstart >= n break @@ -210,19 +223,29 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), debug::Bool=false) wher sort!(d) end -function eigQL!(S::SymTridiagonal{T}; - vectors::Matrix = zeros(T, 0, size(S, 1)), - tol = eps(T), - debug::Bool=false) where T<:Real +function eigQL!( + S::SymTridiagonal{T}; + vectors::Matrix = zeros(T, 0, size(S, 1)), + tol = eps(T), + debug::Bool = false, +) where {T<:Real} d = S.dv e = S.ev n = length(d) if size(vectors, 2) != n - throw(DimensionMismatch("eigenvector matrix must have $(n) columns but had $(size(vectors, 2))")) + throw( + DimensionMismatch( + "eigenvector matrix must have $(n) columns but had $(size(vectors, 2))", + ), + ) end if size(vectors, 1) > n - throw(DimensionMismatch("eigenvector matrix must have at most $(n) rows but had $(size(vectors, 1))")) + throw( + DimensionMismatch( + "eigenvector matrix must have at most $(n) rows but had $(size(vectors, 1))", + ), + ) end blockstart = 1 @@ -231,7 +254,7 @@ function eigQL!(S::SymTridiagonal{T}; while true # Check for zero off diagonal elements for be = blockstart+1:n - if abs(e[be-1]) <= tol*sqrt(abs(d[be-1]))*sqrt(abs(d[be])) + if abs(e[be-1]) <= tol * sqrt(abs(d[be-1])) * sqrt(abs(d[be])) blockend = be - 1 break end @@ -247,13 +270,20 @@ function eigQL!(S::SymTridiagonal{T}; blockstart += 1 else # Calculate shift - μ = (d[blockstart + 1] - d[blockstart])/(2*e[blockstart]) + μ = (d[blockstart+1] - d[blockstart]) / (2 * e[blockstart]) r = hypot(μ, one(T)) - μ = d[blockstart] - (e[blockstart]/(μ + copysign(r, μ))) + μ = d[blockstart] - (e[blockstart] / (μ + copysign(r, μ))) # QL bulk chase singleShiftQL!(S, μ, blockstart, blockend, vectors) - debug && @printf("QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %f\n", blockstart, blockend, e[blockstart], e[blockend-1], μ) + debug && @printf( + "QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %f\n", + blockstart, + blockend, + e[blockstart], + e[blockend-1], + μ + ) end if blockstart >= n break @@ -261,13 +291,15 @@ function eigQL!(S::SymTridiagonal{T}; end end p = sortperm(d) - return d[p], vectors[:,p] + return d[p], vectors[:, p] end -function eigQR!(S::SymTridiagonal{T}; - vectors::Matrix = zeros(T, 0, size(S, 1)), - tol = eps(T), - debug::Bool=false) where T<:Real +function eigQR!( + S::SymTridiagonal{T}; + vectors::Matrix = zeros(T, 0, size(S, 1)), + tol = eps(T), + debug::Bool = false, +) where {T<:Real} d = S.dv e = S.ev n = length(d) @@ -276,8 +308,8 @@ function eigQR!(S::SymTridiagonal{T}; @inbounds begin while true # Check for zero off diagonal elements - for be in (blockstart + 1):n - if abs(e[be - 1]) <= tol*sqrt(abs(d[be - 1]))*sqrt(abs(d[be])) + for be = (blockstart+1):n + if abs(e[be-1]) <= tol * sqrt(abs(d[be-1])) * sqrt(abs(d[be])) blockend = be - 1 break end @@ -293,14 +325,22 @@ function eigQR!(S::SymTridiagonal{T}; blockstart += 1 else # Calculate shift - μ = (d[blockend - 1] - d[blockend])/(2*e[blockend - 1]) + μ = (d[blockend-1] - d[blockend]) / (2 * e[blockend-1]) r = hypot(μ, one(T)) - μ = d[blockend] - (e[blockend - 1]/(μ + copysign(r, μ))) + μ = d[blockend] - (e[blockend-1] / (μ + copysign(r, μ))) # QR bulk chase singleShiftQR!(S, μ, blockstart, blockend, vectors) - debug && @printf("QR, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, d[blockend]: %f, μ: %f\n", blockstart, blockend, e[blockstart], e[blockend-1], d[blockend], μ) + debug && @printf( + "QR, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, d[blockend]: %f, μ: %f\n", + blockstart, + blockend, + e[blockstart], + e[blockend-1], + d[blockend], + μ + ) end if blockstart >= n break @@ -308,14 +348,16 @@ function eigQR!(S::SymTridiagonal{T}; end end p = sortperm(d) - return d[p], vectors[:,p] + return d[p], vectors[:, p] end # Own implementation based on Parlett's book -function singleShiftQLPWK!(S::SymTridiagonal, - shift::Number, - istart::Integer = 1, - iend::Integer = length(S.dv)) +function singleShiftQLPWK!( + S::SymTridiagonal, + shift::Number, + istart::Integer = 1, + iend::Integer = length(S.dv), +) d = S.dv e = S.ev n = length(d) @@ -325,91 +367,95 @@ function singleShiftQLPWK!(S::SymTridiagonal, s = zero(eltype(S)) @inbounds for i = iend-1:-1:istart ei = e[i] - ζ = π + ei - if i < iend-1 - e[i+1] = s*ζ + ζ = π + ei + if i < iend - 1 + e[i+1] = s * ζ end - ci1 = ci - ci = π/ζ - s = ei/ζ - di = d[i] - γi1 = γi - γi = ci*(di - shift) - s*γi1 + ci1 = ci + ci = π / ζ + s = ei / ζ + di = d[i] + γi1 = γi + γi = ci * (di - shift) - s * γi1 d[i+1] = γi1 + di - γi - π = ci == 0 ? ei*ci1 : γi*γi/ci + π = ci == 0 ? ei * ci1 : γi * γi / ci end - e[istart] = s*π + e[istart] = s * π d[istart] = shift + γi S end # Own implementation based on Parlett's book -function singleShiftQL!(S::SymTridiagonal, - shift::Number, - istart::Integer = 1, - iend::Integer = length(S.dv), - vectors = zeros(eltype(S), 0, size(S, 1))) - d, e = S.dv, S.ev - n = length(d) - γi = d[iend] - shift - π = γi +function singleShiftQL!( + S::SymTridiagonal, + shift::Number, + istart::Integer = 1, + iend::Integer = length(S.dv), + vectors = zeros(eltype(S), 0, size(S, 1)), +) + d, e = S.dv, S.ev + n = length(d) + γi = d[iend] - shift + π = γi ci, si = one(eltype(S)), zero(eltype(S)) - @inbounds for i = (iend - 1):-1:istart - ei = e[i] - ci1 = ci - si1 = si + @inbounds for i = (iend-1):-1:istart + ei = e[i] + ci1 = ci + si1 = si ci, si, ζ = givensAlgorithm(π, ei) if i < iend - 1 - e[i + 1] = si1*ζ + e[i+1] = si1 * ζ end - di = d[i] - γi1 = γi - γi = ci*ci*(di - shift) - si*si*γi1 - d[i + 1] = γi1 + di - γi - π = ci == 0 ? -ei*ci1 : γi/ci + di = d[i] + γi1 = γi + γi = ci * ci * (di - shift) - si * si * γi1 + d[i+1] = γi1 + di - γi + π = ci == 0 ? -ei * ci1 : γi / ci # update eigen vectors _updateVectors!(ci, si, i, vectors) end - e[istart] = si*π + e[istart] = si * π d[istart] = shift + γi S end # Own implementation based on Parlett's book -function singleShiftQR!(S::SymTridiagonal, - shift::Number, - istart::Integer = 1, - iend::Integer = length(S.dv), - vectors = zeros(eltype(S), 0, size(S, 1))) - d, e = S.dv, S.ev - n = length(d) - γi = d[istart] - shift - π = γi +function singleShiftQR!( + S::SymTridiagonal, + shift::Number, + istart::Integer = 1, + iend::Integer = length(S.dv), + vectors = zeros(eltype(S), 0, size(S, 1)), +) + d, e = S.dv, S.ev + n = length(d) + γi = d[istart] - shift + π = γi ci, si = one(eltype(S)), zero(eltype(S)) - for i = (istart + 1):iend - ei = e[i - 1] - ci1 = ci - si1 = si + for i = (istart+1):iend + ei = e[i-1] + ci1 = ci + si1 = si ci, si, ζ = givensAlgorithm(π, ei) if i > istart + 1 - e[i - 2] = si1*ζ + e[i-2] = si1 * ζ end - di = d[i] - γi1 = γi - γi = ci*ci*(di - shift) - si*si*γi1 - d[i - 1] = γi1 + di - γi - π = ci == 0 ? -ei*ci1 : γi/ci + di = d[i] + γi1 = γi + γi = ci * ci * (di - shift) - si * si * γi1 + d[i-1] = γi1 + di - γi + π = ci == 0 ? -ei * ci1 : γi / ci # update eigen vectors _updateVectors!(ci, -si, i - 1, vectors) end - e[iend - 1] = si*π - d[iend] = shift + γi + e[iend-1] = si * π + d[iend] = shift + γi S end -function zeroshiftQR!(A::Bidiagonal{T}) where T +function zeroshiftQR!(A::Bidiagonal{T}) where {T} d = A.dv e = A.ev n = length(d) @@ -417,211 +463,225 @@ function zeroshiftQR!(A::Bidiagonal{T}) where T olds = oldc c = oldc for i = 1:n-1 - c, s, r = givensAlgorithm(d[i]*c,e[i]) + c, s, r = givensAlgorithm(d[i] * c, e[i]) if i > 1 - e[i-1] = olds*r + e[i-1] = olds * r end - oldc, olds, d[i] = givensAlgorithm(oldc*r,d[i+1]*s) + oldc, olds, d[i] = givensAlgorithm(oldc * r, d[i+1] * s) end - h = d[n]*c - e[n-1] = h*olds - d[n] = h*oldc + h = d[n] * c + e[n-1] = h * olds + d[n] = h * oldc return A end symtri!(A::Hermitian) = A.uplo == 'L' ? symtriLower!(A.data) : symtriUpper!(A.data) -symtri!(A::Symmetric{T}) where {T<:Real} = A.uplo == 'L' ? symtriLower!(A.data) : symtriUpper!(A.data) +symtri!(A::Symmetric{T}) where {T<:Real} = + A.uplo == 'L' ? symtriLower!(A.data) : symtriUpper!(A.data) # Assume that lower triangle stores the relevant part -function symtriLower!(AS::StridedMatrix{T}, - τ = zeros(T, size(AS, 1) - 1), - u = Vector{T}(undef, size(AS, 1))) where T +function symtriLower!( + AS::StridedMatrix{T}, + τ = zeros(T, size(AS, 1) - 1), + u = Vector{T}(undef, size(AS, 1)), +) where {T} n = size(AS, 1) # We ignore any non-real components of the diagonal - @inbounds for i in 1:n - AS[i,i] = real(AS[i,i]) + @inbounds for i = 1:n + AS[i, i] = real(AS[i, i]) end - @inbounds for k = 1:(n - 2 + !(T<:Real)) + @inbounds for k = 1:(n-2+!(T <: Real)) - τk = LinearAlgebra.reflector!(view(AS, (k + 1):n, k)) + τk = LinearAlgebra.reflector!(view(AS, (k+1):n, k)) τ[k] = τk - for i in (k + 1):n - u[i] = AS[i, k + 1] + for i = (k+1):n + u[i] = AS[i, k+1] end - for j in (k + 2):n + for j = (k+2):n ASjk = AS[j, k] - for i in j:n - u[i] += AS[i, j]*ASjk + for i = j:n + u[i] += AS[i, j] * ASjk end end - for j in (k + 1):(n - 1) + for j = (k+1):(n-1) tmp = zero(T) - for i in j+1:n - tmp += AS[i,j]'AS[i,k] + for i = j+1:n + tmp += AS[i, j]'AS[i, k] end u[j] += tmp end - vcAv = u[k + 1] - for i in (k + 2):n + vcAv = u[k+1] + for i = (k+2):n vcAv += AS[i, k]'u[i] end - ξτ2 = real(vcAv)*abs2(τk)/2 + ξτ2 = real(vcAv) * abs2(τk) / 2 - u[k + 1] = u[k + 1]*τk - ξτ2 - for i in (k + 2):n - u[i] = u[i]*τk - ξτ2*AS[i, k] + u[k+1] = u[k+1] * τk - ξτ2 + for i = (k+2):n + u[i] = u[i] * τk - ξτ2 * AS[i, k] end - AS[k + 1, k + 1] -= 2real(u[k + 1]) - for i in (k + 2):n - AS[i, k + 1] -= u[i] + AS[i, k]*u[k + 1]' + AS[k+1, k+1] -= 2real(u[k+1]) + for i = (k+2):n + AS[i, k+1] -= u[i] + AS[i, k] * u[k+1]' end - for j in (k + 2):n + for j = (k+2):n ASjk = AS[j, k] uj = u[j] - AS[j,j] -= 2real(uj*ASjk') - for i = (j + 1):n - AS[i, j] -= u[i]*ASjk' + AS[i, k]*uj' + AS[j, j] -= 2real(uj * ASjk') + for i = (j+1):n + AS[i, j] -= u[i] * ASjk' + AS[i, k] * uj' end end end - SymmetricTridiagonalFactorization('L', AS, τ, SymTridiagonal(real(diag(AS)), real(diag(AS, -1)))) + SymmetricTridiagonalFactorization( + 'L', + AS, + τ, + SymTridiagonal(real(diag(AS)), real(diag(AS, -1))), + ) end # Assume that upper triangle stores the relevant part -function symtriUpper!(AS::StridedMatrix{T}, - τ = zeros(T, size(AS, 1) - 1), - u = Vector{T}(undef, size(AS, 1))) where T +function symtriUpper!( + AS::StridedMatrix{T}, + τ = zeros(T, size(AS, 1) - 1), + u = Vector{T}(undef, size(AS, 1)), +) where {T} n = LinearAlgebra.checksquare(AS) # We ignore any non-real components of the diagonal - @inbounds for i in 1:n - AS[i,i] = real(AS[i,i]) + @inbounds for i = 1:n + AS[i, i] = real(AS[i, i]) end - @inbounds for k = 1:(n - 2 + !(T<:Real)) + @inbounds for k = 1:(n-2+!(T <: Real)) # This is a bit convoluted method to get the conjugated vector but conjugation is required for correctness of arrays of quaternions. Eventually, it should be sufficient to write vec(x') but it currently (July 10, 2018) hits a bug in LinearAlgebra - τk = LinearAlgebra.reflector!(vec(transpose(view(AS, k, (k + 1):n)'))) + τk = LinearAlgebra.reflector!(vec(transpose(view(AS, k, (k+1):n)'))) τ[k] = τk' - for j in (k + 1):n - tmp = AS[k + 1, j] - for i in (k + 2):j - tmp += AS[k, i]*AS[i, j] + for j = (k+1):n + tmp = AS[k+1, j] + for i = (k+2):j + tmp += AS[k, i] * AS[i, j] end - for i in (j + 1):n - tmp += AS[k, i]*AS[j, i]' + for i = (j+1):n + tmp += AS[k, i] * AS[j, i]' end - u[j] = τk'*tmp + u[j] = τk' * tmp end - vcAv = u[k + 1] - for i in (k + 2):n - vcAv += u[i]*AS[k, i]' + vcAv = u[k+1] + for i = (k+2):n + vcAv += u[i] * AS[k, i]' end - ξ = real(vcAv*τk) + ξ = real(vcAv * τk) - for j in (k + 1):n - ujt = u[j] - hjt = j > (k + 1) ? AS[k, j] : one(ujt) - ξhjt = ξ*hjt - for i in (k + 1):(j - 1) + for j = (k+1):n + ujt = u[j] + hjt = j > (k + 1) ? AS[k, j] : one(ujt) + ξhjt = ξ * hjt + for i = (k+1):(j-1) hit = i > (k + 1) ? AS[k, i] : one(ujt) - AS[i, j] -= hit'*ujt + u[i]'*hjt - hit'*ξhjt + AS[i, j] -= hit' * ujt + u[i]' * hjt - hit' * ξhjt end - AS[j, j] -= 2*real(hjt'*ujt) - abs2(hjt)*ξ + AS[j, j] -= 2 * real(hjt' * ujt) - abs2(hjt) * ξ end end - SymmetricTridiagonalFactorization('U', AS, τ, SymTridiagonal(real(diag(AS)), real(diag(AS, 1)))) + SymmetricTridiagonalFactorization( + 'U', + AS, + τ, + SymTridiagonal(real(diag(AS)), real(diag(AS, 1))), + ) end -_eigvals!(A::SymmetricTridiagonalFactorization; - tol = eps(real(eltype(A))), - debug = false) = eigvalsPWK!(A.diagonals, tol = tol, debug = debug) +_eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), debug = false) = + eigvalsPWK!(A.diagonals, tol = tol, debug = debug) -_eigvals!(A::SymTridiagonal; - tol = eps(real(eltype(A))), - debug = false) = eigvalsPWK!(A, tol = tol, debug = debug) +_eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), debug = false) = + eigvalsPWK!(A, tol = tol, debug = debug) -_eigvals!(A::Hermitian; - tol = eps(real(eltype(A))), - debug = false) = eigvals!(symtri!(A), tol = tol, debug = debug) +_eigvals!(A::Hermitian; tol = eps(real(eltype(A))), debug = false) = + eigvals!(symtri!(A), tol = tol, debug = debug) -LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization; - tol = eps(real(eltype(A))), - debug = false) = _eigvals!(A, tol = tol, debug = debug) +LinearAlgebra.eigvals!( + A::SymmetricTridiagonalFactorization; + tol = eps(real(eltype(A))), + debug = false, +) = _eigvals!(A, tol = tol, debug = debug) -LinearAlgebra.eigvals!(A::SymTridiagonal; - tol = eps(real(eltype(A))), - debug = false) = _eigvals!(A, tol = tol, debug = debug) +LinearAlgebra.eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), debug = false) = + _eigvals!(A, tol = tol, debug = debug) -LinearAlgebra.eigvals!(A::Hermitian; - tol = eps(real(eltype(A))), - debug = false) = _eigvals!(A, tol = tol, debug = debug) +LinearAlgebra.eigvals!(A::Hermitian; tol = eps(real(eltype(A))), debug = false) = + _eigvals!(A, tol = tol, debug = debug) -_eigen!(A::SymmetricTridiagonalFactorization; - tol = eps(real(eltype(A))), - debug = false) = - LinearAlgebra.Eigen(eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol, debug = debug)...) +_eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), debug = false) = + LinearAlgebra.Eigen( + eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol, debug = debug)..., + ) -_eigen!(A::SymTridiagonal; - tol = eps(real(eltype(A))), - debug = false) = - LinearAlgebra.Eigen(eigQL!(A, vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), tol = tol, debug = debug)...) +_eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), debug = false) = LinearAlgebra.Eigen( + eigQL!( + A, + vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), + tol = tol, + debug = debug, + )..., +) -_eigen!(A::Hermitian; - tol = eps(real(eltype(A))), - debug = false) = _eigen!(symtri!(A), tol = tol, debug = debug) +_eigen!(A::Hermitian; tol = eps(real(eltype(A))), debug = false) = + _eigen!(symtri!(A), tol = tol, debug = debug) -LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization; - tol = eps(real(eltype(A))), - debug = false) = _eigen!(A, tol = tol, debug = debug) +LinearAlgebra.eigen!( + A::SymmetricTridiagonalFactorization; + tol = eps(real(eltype(A))), + debug = false, +) = _eigen!(A, tol = tol, debug = debug) -LinearAlgebra.eigen!(A::SymTridiagonal; - tol = eps(real(eltype(A))), - debug = false) = _eigen!(A, tol = tol, debug = debug) +LinearAlgebra.eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), debug = false) = + _eigen!(A, tol = tol, debug = debug) -LinearAlgebra.eigen!(A::Hermitian; - tol = eps(real(eltype(A))), - debug = false) = _eigen!(A, tol = tol, debug = debug) +LinearAlgebra.eigen!(A::Hermitian; tol = eps(real(eltype(A))), debug = false) = + _eigen!(A, tol = tol, debug = debug) -function eigen2!(A::SymmetricTridiagonalFactorization; - tol = eps(real(float(one(eltype(A))))), - debug = false) +function eigen2!( + A::SymmetricTridiagonalFactorization; + tol = eps(real(float(one(eltype(A))))), + debug = false, +) V = zeros(eltype(A), 2, size(A, 1)) V[1] = 1 V[end] = 1 eigQL!(A.diagonals, vectors = rmul!(V, A.Q), tol = tol, debug = debug) end -function eigen2!(A::SymTridiagonal; - tol = eps(real(float(one(eltype(A))))), - debug = false) +function eigen2!(A::SymTridiagonal; tol = eps(real(float(one(eltype(A))))), debug = false) V = zeros(eltype(A), 2, size(A, 1)) V[1] = 1 V[end] = 1 eigQL!(A, vectors = V, tol = tol, debug = debug) end -eigen2!(A::Hermitian; - tol = eps(float(real(one(eltype(A))))), - debug = false) = eigen2!(symtri!(A), tol = tol, debug = debug) +eigen2!(A::Hermitian; tol = eps(float(real(one(eltype(A))))), debug = false) = + eigen2!(symtri!(A), tol = tol, debug = debug) -eigen2(A::SymTridiagonal; - tol = eps(float(real(one(eltype(A))))), - debug = false) = eigen2!(copy(A), tol = tol, debug = debug) +eigen2(A::SymTridiagonal; tol = eps(float(real(one(eltype(A))))), debug = false) = + eigen2!(copy(A), tol = tol, debug = debug) -eigen2(A::Hermitian , tol = eps(float(real(one(eltype(A))))), debug = false) = eigen2!(copy(A), tol = tol, debug = debug) +eigen2(A::Hermitian, tol = eps(float(real(one(eltype(A))))), debug = false) = + eigen2!(copy(A), tol = tol, debug = debug) # First method of each type here is identical to the method defined in # LinearAlgebra but is needed for disambiguation @@ -653,5 +713,5 @@ end # Aux (should go somewhere else at some point) function LinearAlgebra.givensAlgorithm(f::Real, g::Real) h = hypot(f, g) - return f/h, g/h, h + return f / h, g / h, h end diff --git a/src/householder.jl b/src/householder.jl index 34913c4..b34e1f1 100644 --- a/src/householder.jl +++ b/src/householder.jl @@ -13,23 +13,27 @@ end size(H::Householder) = (length(H.v) + 1, length(H.v) + 1) size(H::Householder, i::Integer) = i <= 2 ? length(H.v) + 1 : 1 -eltype(H::Householder{T}) where T = T -eltype(H::HouseholderBlock{T}) where T = T +eltype(H::Householder{T}) where {T} = T +eltype(H::HouseholderBlock{T}) where {T} = T -adjoint(H::Householder{T}) where {T} = Adjoint{T,typeof(H)}(H) +adjoint(H::Householder{T}) where {T} = Adjoint{T,typeof(H)}(H) adjoint(H::HouseholderBlock{T}) where {T} = Adjoint{T,typeof(H)}(H) function lmul!(H::Householder, A::StridedMatrix) m, n = size(A) - length(H.v) == m - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))")) - v = view(H.v, 1:m - 1) + length(H.v) == m - 1 || throw( + DimensionMismatch( + "size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))", + ), + ) + v = view(H.v, 1:m-1) τ = H.τ for j = 1:n - va = A[1,j] + va = A[1, j] Aj = view(A, 2:m, j) va += dot(v, Aj) - va = τ*va - A[1,j] -= va + va = τ * va + A[1, j] -= va axpy!(-va, v, Aj) end A @@ -37,12 +41,16 @@ end function rmul!(A::StridedMatrix, H::Householder) m, n = size(A) - length(H.v) == n - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but second dimension of matrix is $(size(A, 2))")) + length(H.v) == n - 1 || throw( + DimensionMismatch( + "size of reflector is $(length(H.v) + 1) but second dimension of matrix is $(size(A, 2))", + ), + ) v = view(H.v, :) τ = H.τ a1 = view(A, :, 1) A1 = view(A, :, 2:n) - x = A1*v + x = A1 * v axpy!(one(τ), a1, x) axpy!(-τ, x, a1) rankUpdate!(A1, x, v, -τ) @@ -52,22 +60,26 @@ end function lmul!(adjH::Adjoint{<:Any,<:Householder}, A::StridedMatrix) H = parent(adjH) m, n = size(A) - length(H.v) == m - 1 || throw(DimensionMismatch("size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))")) - v = view(H.v, 1:m - 1) + length(H.v) == m - 1 || throw( + DimensionMismatch( + "size of reflector is $(length(H.v) + 1) but first dimension of matrix is $(size(A, 1))", + ), + ) + v = view(H.v, 1:m-1) τ = H.τ for j = 1:n - va = A[1,j] + va = A[1, j] Aj = view(A, 2:m, j) va += dot(v, Aj) va = τ'va - A[1,j] -= va + A[1, j] -= va axpy!(-va, v, Aj) end A end # FixMe! This is a weird multiplication method and at least its behavior needs to be explained -function lmul!(H::HouseholderBlock{T}, A::StridedMatrix{T}, M::StridedMatrix{T}) where T +function lmul!(H::HouseholderBlock{T}, A::StridedMatrix{T}, M::StridedMatrix{T}) where {T} V = H.V mA, nA = size(A) nH = size(V, 1) @@ -102,9 +114,13 @@ function lmul!(H::HouseholderBlock{T}, A::StridedMatrix{T}, M::StridedMatrix{T}) return A end (*)(H::HouseholderBlock{T}, A::StridedMatrix{T}) where {T} = - lmul!(H, copy(A), similar(A, (min(size(H.V)...), size(A, 2)))) + lmul!(H, copy(A), similar(A, (min(size(H.V)...), size(A, 2)))) -function lmul!(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}, M::StridedMatrix) where T +function lmul!( + adjH::Adjoint{T,<:HouseholderBlock{T}}, + A::StridedMatrix{T}, + M::StridedMatrix, +) where {T} H = parent(adjH) V = H.V mA, nA = size(A) @@ -140,7 +156,8 @@ function lmul!(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}, M::S return A end (*)(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}) where {T} = - lmul!(adjH, copy(A), similar(A, (min(size(parent(adjH).V)...), size(A, 2)))) + lmul!(adjH, copy(A), similar(A, (min(size(parent(adjH).V)...), size(A, 2)))) Base.convert(::Type{Matrix}, H::Householder{T}) where {T} = convert(Matrix{T}, H) -Base.convert(::Type{Matrix{T}}, H::Householder) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1))) +Base.convert(::Type{Matrix{T}}, H::Householder) where {T} = + lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1))) diff --git a/src/juliaBLAS.jl b/src/juliaBLAS.jl index fcc188b..cb40e8b 100644 --- a/src/juliaBLAS.jl +++ b/src/juliaBLAS.jl @@ -9,7 +9,12 @@ export rankUpdate! ## General ### BLAS -rankUpdate!(A::StridedMatrix{T}, x::StridedVector{T}, y::StridedVector{T}, α::T) where {T<:BlasReal} = BLAS.ger!(α, x, y, A) +rankUpdate!( + A::StridedMatrix{T}, + x::StridedVector{T}, + y::StridedVector{T}, + α::T, +) where {T<:BlasReal} = BLAS.ger!(α, x, y, A) ### Generic function rankUpdate!(A::StridedMatrix, x::StridedVector, y::StridedVector, α::Number) @@ -19,30 +24,39 @@ function rankUpdate!(A::StridedMatrix, x::StridedVector, y::StridedVector, α::N for j = 1:n yjc = y[j]' for i = 1:m - A[i,j] += x[i]*α*yjc + A[i, j] += x[i] * α * yjc end end end ## Hermitian -function rankUpdate!(A::HermOrSym{T,S}, a::StridedVector{T}, α::T) where {T<:BlasReal,S<:StridedMatrix} +function rankUpdate!( + A::HermOrSym{T,S}, + a::StridedVector{T}, + α::T, +) where {T<:BlasReal,S<:StridedMatrix} BLAS.syr!(A.uplo, α, a, A.data) return A end -function rankUpdate!(A::Hermitian{Complex{T},S}, a::StridedVector{Complex{T}}, α::T) where {T<:BlasReal,S<:StridedMatrix} +function rankUpdate!( + A::Hermitian{Complex{T},S}, + a::StridedVector{Complex{T}}, + α::T, +) where {T<:BlasReal,S<:StridedMatrix} BLAS.her!(A.uplo, α, a, A.data) return A end -rankUpdate!(A::HermOrSym{T,S}, a::StridedVector{T}) where {T<:BlasFloat,S<:StridedMatrix} = rankUpdate!(A, a, one(real(T))) +rankUpdate!(A::HermOrSym{T,S}, a::StridedVector{T}) where {T<:BlasFloat,S<:StridedMatrix} = + rankUpdate!(A, a, one(real(T))) ### Generic function rankUpdate!(A::Hermitian, a::StridedVector, α::Real) n = size(A, 1) n == length(a) || throw(DimensionMismatch("a vector has wrong length")) - @inbounds for j in 1:n + @inbounds for j = 1:n ajc = a[j]' for i in ((A.uplo == 'L') ? (j:n) : (1:j)) - A.data[i,j] += a[i]*α*ajc + A.data[i, j] += a[i] * α * ajc end end return A @@ -50,13 +64,23 @@ end # Rank k update ## Real -function rankUpdate!(C::HermOrSym{T,S}, A::StridedMatrix{T}, α::T, β::T) where {T<:BlasReal,S<:StridedMatrix} +function rankUpdate!( + C::HermOrSym{T,S}, + A::StridedMatrix{T}, + α::T, + β::T, +) where {T<:BlasReal,S<:StridedMatrix} BLAS.syrk!(C.uplo, 'N', α, A, β, C.data) return C end ## Complex -function rankUpdate!(C::Hermitian{Complex{T},S}, A::StridedMatrix{Complex{T}}, α::T, β::T) where {T<:BlasReal,S<:StridedMatrix} +function rankUpdate!( + C::Hermitian{Complex{T},S}, + A::StridedMatrix{Complex{T}}, + α::T, + β::T, +) where {T<:BlasReal,S<:StridedMatrix} BLAS.herk!(C.uplo, 'N', α, A, β, C.data) return C end @@ -66,20 +90,20 @@ function rankUpdate!(C::Hermitian, A::StridedVecOrMat, α::Real) n = size(C, 1) n == size(A, 1) || throw(DimensionMismatch("first dimension of A has wrong size")) @inbounds if C.uplo == 'L' # branch outside the loop to have larger loop to optimize - for k in 1:size(A, 2) - for j in 1:n - ajkc = A[j,k]' - for i in j:n - C.data[i,j] += A[i,k]*α*ajkc + for k = 1:size(A, 2) + for j = 1:n + ajkc = A[j, k]' + for i = j:n + C.data[i, j] += A[i, k] * α * ajkc end end end else - for k in 1:size(A, 2) - for j in 1:n - ajkc = A[j,k]' - for i in 1:j - C.data[i,j] += A[i,k]*α*ajkc + for k = 1:size(A, 2) + for j = 1:n + ajkc = A[j, k]' + for i = 1:j + C.data[i, j] += A[i, k] * α * ajkc end end end @@ -93,9 +117,9 @@ function lmul!(A::UpperTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:Nu m, n = size(B) for i ∈ 1:m for j ∈ 1:n - B[i,j] = α*AA[i,i]*B[i,j] - for l ∈ (i + 1):m - B[i,j] += α*AA[i,l]*B[l,j] + B[i, j] = α * AA[i, i] * B[i, j] + for l ∈ (i+1):m + B[i, j] += α * AA[i, l] * B[l, j] end end end @@ -106,9 +130,9 @@ function lmul!(A::LowerTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:Nu m, n = size(B) for i ∈ m:-1:1 for j ∈ 1:n - B[i,j] = α*AA[i,i]*B[i,j] - for l ∈ 1:(i - 1) - B[i,j] += α*AA[i,l]*B[l,j] + B[i, j] = α * AA[i, i] * B[i, j] + for l ∈ 1:(i-1) + B[i, j] += α * AA[i, l] * B[l, j] end end end @@ -119,9 +143,9 @@ function lmul!(A::UnitUpperTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T m, n = size(B) for i ∈ 1:m for j ∈ 1:n - B[i,j] = α*B[i,j] - for l ∈ (i + 1):m - B[i,j] += α*AA[i,l]*B[l,j] + B[i, j] = α * B[i, j] + for l ∈ (i+1):m + B[i, j] += α * AA[i, l] * B[l, j] end end end @@ -132,61 +156,77 @@ function lmul!(A::UnitLowerTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T m, n = size(B) for i ∈ m:-1:1 for j ∈ 1:n - B[i,j] = α*B[i,j] - for l ∈ 1:(i - 1) - B[i,j] += α*AA[i,l]*B[l,j] + B[i, j] = α * B[i, j] + for l ∈ 1:(i-1) + B[i, j] += α * AA[i, l] * B[l, j] end end end return B end -function lmul!(A::LowerTriangular{T,Adjoint{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S} +function lmul!( + A::LowerTriangular{T,Adjoint{T,S}}, + B::StridedMatrix{T}, + α::T, +) where {T<:Number,S} AA = parent(A.data) m, n = size(B) for i ∈ m:-1:1 for j ∈ 1:n - B[i,j] = α*AA[i,i]'*B[i,j] - for l ∈ 1:(i - 1) - B[i,j] += α*AA[l,i]'*B[l,j] + B[i, j] = α * AA[i, i]' * B[i, j] + for l ∈ 1:(i-1) + B[i, j] += α * AA[l, i]' * B[l, j] end end end return B end -function lmul!(A::UpperTriangular{T,Adjoint{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S} +function lmul!( + A::UpperTriangular{T,Adjoint{T,S}}, + B::StridedMatrix{T}, + α::T, +) where {T<:Number,S} AA = parent(A.data) m, n = size(B) for i ∈ 1:m for j ∈ 1:n - B[i,j] = α*AA[i,i]'*B[i,j] - for l ∈ (i + 1):m - B[i,j] += α*AA[l,i]'*B[l,j] + B[i, j] = α * AA[i, i]' * B[i, j] + for l ∈ (i+1):m + B[i, j] += α * AA[l, i]' * B[l, j] end end end return B end -function lmul!(A::UnitLowerTriangular{T,Adjoint{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S} +function lmul!( + A::UnitLowerTriangular{T,Adjoint{T,S}}, + B::StridedMatrix{T}, + α::T, +) where {T<:Number,S} AA = parent(A.data) m, n = size(B) for i ∈ m:-1:1 for j ∈ 1:n - B[i,j] = α*B[i,j] - for l ∈ 1:(i - 1) - B[i,j] += α*AA[l,i]'*B[l,j] + B[i, j] = α * B[i, j] + for l ∈ 1:(i-1) + B[i, j] += α * AA[l, i]' * B[l, j] end end end return B end -function lmul!(A::UnitUpperTriangular{T,Adjoint{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S} +function lmul!( + A::UnitUpperTriangular{T,Adjoint{T,S}}, + B::StridedMatrix{T}, + α::T, +) where {T<:Number,S} AA = parent(A.data) m, n = size(B) for i ∈ 1:m for j ∈ 1:n - B[i,j] = α*B[i,j] - for l ∈ (i + 1):m - B[i,j] += α*AA[l,i]'*B[l,j] + B[i, j] = α * B[i, j] + for l ∈ (i+1):m + B[i, j] += α * AA[l, i]' * B[l, j] end end end diff --git a/src/lapack.jl b/src/lapack.jl index 7c5a130..571666f 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -1,402 +1,620 @@ module LAPACK2 - using libblastrampoline_jll - using LinearAlgebra - using LinearAlgebra: BlasInt, chkstride1, LAPACKException - using LinearAlgebra.BLAS: @blasfunc - using LinearAlgebra.LAPACK: chkdiag, chkside, chkuplo - - liblapack_name = libblastrampoline_jll.libblastrampoline - - ## Standard QR/QL - function steqr!(compz::Char, - d::StridedVector{Float64}, - e::StridedVector{Float64}, - Z::StridedMatrix{Float64} = compz == 'N' ? - Matrix{Float64}(undef, 0, 0) : - Matrix{Float64}(undef, length(d), length(d)), - work::StridedVector{Float64} = compz == 'N' ? - Vector{Float64}(undef, 0) : - Vector{Float64}(undef, max(1, 2*length(d) - 2))) - - # Extract sizes - n = length(d) - ldz = stride(Z, 2) - - # Checks - length(e) + 1 == n || throw(DimensionMismatch("")) - if compz == 'N' - elseif compz in ('V', 'I') - chkstride1(Z) - n <= size(Z, 1) || throw(DimensionMismatch("Z matrix has too few rows")) - n <= size(Z, 2) || throw(DimensionMismatch("Z matrix has too few colums")) - else - throw(ArgumentError("compz must be either 'N', 'V', or 'I'")) - end - - # Allocations - info = Vector{BlasInt}(undef, 1) - - ccall((@blasfunc("dsteqr_"), liblapack_name),Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64}, - Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ptr{BlasInt}), - compz, n, d, e, - Z, max(1, ldz), work, info) - info[1] == 0 || throw(LAPACKException(info[1])) - - return d, Z +using libblastrampoline_jll +using LinearAlgebra +using LinearAlgebra: BlasInt, chkstride1, LAPACKException +using LinearAlgebra.BLAS: @blasfunc +using LinearAlgebra.LAPACK: chkdiag, chkside, chkuplo + +liblapack_name = libblastrampoline_jll.libblastrampoline + +## Standard QR/QL +function steqr!( + compz::Char, + d::StridedVector{Float64}, + e::StridedVector{Float64}, + Z::StridedMatrix{Float64} = compz == 'N' ? Matrix{Float64}(undef, 0, 0) : + Matrix{Float64}(undef, length(d), length(d)), + work::StridedVector{Float64} = compz == 'N' ? Vector{Float64}(undef, 0) : + Vector{Float64}(undef, max(1, 2 * length(d) - 2)), +) + + # Extract sizes + n = length(d) + ldz = stride(Z, 2) + + # Checks + length(e) + 1 == n || throw(DimensionMismatch("")) + if compz == 'N' + elseif compz in ('V', 'I') + chkstride1(Z) + n <= size(Z, 1) || throw(DimensionMismatch("Z matrix has too few rows")) + n <= size(Z, 2) || throw(DimensionMismatch("Z matrix has too few colums")) + else + throw(ArgumentError("compz must be either 'N', 'V', or 'I'")) end - ## Square root free QR/QL (values only, but fast) - function sterf!(d::StridedVector{Float64}, e::StridedVector{Float64}) - - # Extract sizes - n = length(d) - - # Checks - length(e) >= n - 1 || throw(DimensionMismatch("subdiagonal is too short")) + # Allocations + info = Vector{BlasInt}(undef, 1) + + ccall( + (@blasfunc("dsteqr_"), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{BlasInt}, + Ptr{Float64}, + Ptr{Float64}, + Ptr{Float64}, + Ref{BlasInt}, + Ptr{Float64}, + Ptr{BlasInt}, + ), + compz, + n, + d, + e, + Z, + max(1, ldz), + work, + info, + ) + info[1] == 0 || throw(LAPACKException(info[1])) + + return d, Z +end - # Allocations - info = BlasInt[0] +## Square root free QR/QL (values only, but fast) +function sterf!(d::StridedVector{Float64}, e::StridedVector{Float64}) - ccall((@blasfunc("dsterf_"), liblapack_name), Cvoid, - (Ref{BlasInt}, Ptr{Float64}, Ptr{Float64}, Ptr{BlasInt}), - n, d, e, info) + # Extract sizes + n = length(d) - info[1] == 0 || throw(LAPACKException(info[1])) + # Checks + length(e) >= n - 1 || throw(DimensionMismatch("subdiagonal is too short")) - d - end + # Allocations + info = BlasInt[0] - ## Divide and Conquer - function stedc!(compz::Char, - d::StridedVector{Float64}, - e::StridedVector{Float64}, - Z::StridedMatrix{Float64}, - work::StridedVector{Float64}, - lwork::BlasInt, - iwork::StridedVector{BlasInt}, - liwork::BlasInt) - - # Extract sizes - n = length(d) - ldz = stride(Z, 2) - - # Checks - length(e) + 1 == n || throw(DimensionMismatch("")) - # length(work) - # length(lwork) - - # Allocations - info = BlasInt[0] - - ccall((@blasfunc("dstedc_"), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64}, - Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, - Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}), - compz, n, d, e, - Z, max(1, ldz), work, lwork, - iwork, liwork, info) - - info[1] == 0 || throw(LAPACKException(info[1])) - - return d, Z - end + ccall( + (@blasfunc("dsterf_"), liblapack_name), + Cvoid, + (Ref{BlasInt}, Ptr{Float64}, Ptr{Float64}, Ptr{BlasInt}), + n, + d, + e, + info, + ) - function stedc!(compz::Char, - d::StridedVector{Float64}, - e::StridedVector{Float64}, - Z::StridedMatrix{Float64} = compz == 'N' ? - Matrix{Float64}(undef, 0, 0) : - Matrix{Float64}(undef, length(d), length(d))) + info[1] == 0 || throw(LAPACKException(info[1])) - work::Vector{Float64} = Float64[0] - iwork::Vector{BlasInt} = BlasInt[0] - stedc!(compz, d, e, Z, work, -1, iwork, -1) + d +end - lwork::BlasInt = work[1] - work = Vector{Float64}(undef, lwork) - liwork = iwork[1] - iwork = Vector{BlasInt}(undef, liwork) +## Divide and Conquer +function stedc!( + compz::Char, + d::StridedVector{Float64}, + e::StridedVector{Float64}, + Z::StridedMatrix{Float64}, + work::StridedVector{Float64}, + lwork::BlasInt, + iwork::StridedVector{BlasInt}, + liwork::BlasInt, +) + + # Extract sizes + n = length(d) + ldz = stride(Z, 2) + + # Checks + length(e) + 1 == n || throw(DimensionMismatch("")) + # length(work) + # length(lwork) + + # Allocations + info = BlasInt[0] + + ccall( + (@blasfunc("dstedc_"), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{BlasInt}, + Ptr{Float64}, + Ptr{Float64}, + Ptr{Float64}, + Ref{BlasInt}, + Ptr{Float64}, + Ref{BlasInt}, + Ptr{BlasInt}, + Ref{BlasInt}, + Ptr{BlasInt}, + ), + compz, + n, + d, + e, + Z, + max(1, ldz), + work, + lwork, + iwork, + liwork, + info, + ) + + info[1] == 0 || throw(LAPACKException(info[1])) + + return d, Z +end - return stedc!(compz, d, e, Z, work, lwork, iwork, liwork) - end +function stedc!( + compz::Char, + d::StridedVector{Float64}, + e::StridedVector{Float64}, + Z::StridedMatrix{Float64} = compz == 'N' ? Matrix{Float64}(undef, 0, 0) : + Matrix{Float64}(undef, length(d), length(d)), +) - ## RRR - for (lsymb, elty) in ((:dstemr_, :Float64), (:sstemr_, :Float32)) - @eval begin - function stemr!(jobz::Char, - range::Char, - dv::StridedVector{$elty}, - ev::StridedVector{$elty}, - vl::$elty, - vu::$elty, - il::BlasInt, - iu::BlasInt, - w::StridedVector{$elty}, - Z::StridedMatrix{$elty}, - nzc::BlasInt, - isuppz::StridedVector{BlasInt}, - work::StridedVector{$elty}, - lwork::BlasInt, - iwork::StridedVector{BlasInt}, - liwork::BlasInt) - - # Extract sizes - n = length(dv) - ldz = stride(Z, 2) - - # Checks - length(ev) >= n - 1 || throw(DimensionMismatch("subdiagonal is too short")) - - # Allocations - eev::Vector{$elty} = length(ev) == n - 1 ? [ev; zero($elty)] : copy(ev) - abstol = Vector{$elty}(undef, 1) - m = Vector{BlasInt}(undef, 1) - tryrac = BlasInt[1] - info = Vector{BlasInt}(undef, 1) - - ccall((@blasfunc($lsymb), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ref{$elty}, Ref{$elty}, Ref{BlasInt}, - Ref{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, - Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, - Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}, - Ptr{BlasInt}), - jobz, range, n, dv, - eev, vl, vu, il, - iu, m, w, Z, - max(1, ldz), nzc, isuppz, tryrac, - work, lwork, iwork, liwork, - info) - - info[1] == 0 || throw(LAPACKException(info[1])) - - w, Z, tryrac[1] - end + work::Vector{Float64} = Float64[0] + iwork::Vector{BlasInt} = BlasInt[0] + stedc!(compz, d, e, Z, work, -1, iwork, -1) - function stemr!(jobz::Char, - range::Char, - dv::StridedVector{$elty}, - ev::StridedVector{$elty}, - vl::$elty = typemin($elty), - vu::$elty = typemax($elty), - il::BlasInt = 1, - iu::BlasInt = length(dv)) - - n = length(dv) - w = Vector{$elty}(undef, n) - if jobz == 'N' - # Z = Matrix{$elty}(0, 0) - nzc::BlasInt = 0 - isuppz = Vector{BlasInt}(undef, 0) - elseif jobz == 'V' - nzc = -1 - isuppz = similar(dv, BlasInt, 2n) - else - throw(ArgumentError("jobz must be either 'N' or 'V'")) - end + lwork::BlasInt = work[1] + work = Vector{Float64}(undef, lwork) + liwork = iwork[1] + iwork = Vector{BlasInt}(undef, liwork) - work = Vector{$elty}(undef, 1) - lwork::BlasInt = -1 - iwork = BlasInt[0] - liwork::BlasInt = -1 - Z = Matrix{$elty}(undef, jobz == 'N' ? 1 : n, 1) - nzc = -1 + return stedc!(compz, d, e, Z, work, lwork, iwork, liwork) +end - stemr!(jobz, range, dv, ev, vl, vu, il, iu, w, Z, nzc, isuppz, work, lwork, iwork, liwork) +## RRR +for (lsymb, elty) in ((:dstemr_, :Float64), (:sstemr_, :Float32)) + @eval begin + function stemr!( + jobz::Char, + range::Char, + dv::StridedVector{$elty}, + ev::StridedVector{$elty}, + vl::$elty, + vu::$elty, + il::BlasInt, + iu::BlasInt, + w::StridedVector{$elty}, + Z::StridedMatrix{$elty}, + nzc::BlasInt, + isuppz::StridedVector{BlasInt}, + work::StridedVector{$elty}, + lwork::BlasInt, + iwork::StridedVector{BlasInt}, + liwork::BlasInt, + ) + + # Extract sizes + n = length(dv) + ldz = stride(Z, 2) + + # Checks + length(ev) >= n - 1 || throw(DimensionMismatch("subdiagonal is too short")) + + # Allocations + eev::Vector{$elty} = length(ev) == n - 1 ? [ev; zero($elty)] : copy(ev) + abstol = Vector{$elty}(undef, 1) + m = Vector{BlasInt}(undef, 1) + tryrac = BlasInt[1] + info = Vector{BlasInt}(undef, 1) - lwork = work[1] - work = Vector{$elty}(undef, lwork) - liwork = iwork[1] - iwork = Vector{BlasInt}(undef, liwork) - nzc = Z[1] - Z = similar(dv, $elty, n, nzc) + ccall( + (@blasfunc($lsymb), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ptr{$elty}, + Ptr{$elty}, + Ref{$elty}, + Ref{$elty}, + Ref{BlasInt}, + Ref{BlasInt}, + Ptr{BlasInt}, + Ptr{$elty}, + Ptr{$elty}, + Ref{BlasInt}, + Ref{BlasInt}, + Ptr{BlasInt}, + Ptr{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{BlasInt}, + Ref{BlasInt}, + Ptr{BlasInt}, + ), + jobz, + range, + n, + dv, + eev, + vl, + vu, + il, + iu, + m, + w, + Z, + max(1, ldz), + nzc, + isuppz, + tryrac, + work, + lwork, + iwork, + liwork, + info, + ) + + info[1] == 0 || throw(LAPACKException(info[1])) + + w, Z, tryrac[1] + end - return stemr!(jobz, range, dv, ev, vl, vu, il, iu, w, Z, nzc, isuppz, work, lwork, iwork, liwork) + function stemr!( + jobz::Char, + range::Char, + dv::StridedVector{$elty}, + ev::StridedVector{$elty}, + vl::$elty = typemin($elty), + vu::$elty = typemax($elty), + il::BlasInt = 1, + iu::BlasInt = length(dv), + ) + + n = length(dv) + w = Vector{$elty}(undef, n) + if jobz == 'N' + # Z = Matrix{$elty}(0, 0) + nzc::BlasInt = 0 + isuppz = Vector{BlasInt}(undef, 0) + elseif jobz == 'V' + nzc = -1 + isuppz = similar(dv, BlasInt, 2n) + else + throw(ArgumentError("jobz must be either 'N' or 'V'")) end - end - end - # Non-symmetric eigenvalue problem - function lahqr!(wantt::Bool, - wantz::Bool, - H::StridedMatrix{Float64}, - ilo::Integer, - ihi::Integer, - wr::StridedVector{Float64}, - wi::StridedVector{Float64}, - iloz::Integer, - ihiz::Integer, - Z::StridedMatrix{Float64}) - # SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, - # ILOZ, IHIZ, Z, LDZ, INFO ) - - # .. Scalar Arguments .. - # INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N - # LOGICAL WANTT, WANTZ - # .. - # .. Array Arguments .. - # DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) - - n = LinearAlgebra.checksquare(H) - ldh = stride(H, 2) - ldz = stride(Z, 2) - - info = Ref{BlasInt}(0) - - ccall((@blasfunc("dlahqr_"), liblapack_name), Cvoid, - (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, - Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, - Ptr{Float64}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, - Ref{BlasInt}, Ptr{BlasInt}), - wantt, wantz, n, ilo, - ihi, H, max(1, ldh), wr, - wi, iloz, ihiz, Z, - max(1, ldz), info) - - if info[] != 0 - throw(LAPACKException(info[])) + work = Vector{$elty}(undef, 1) + lwork::BlasInt = -1 + iwork = BlasInt[0] + liwork::BlasInt = -1 + Z = Matrix{$elty}(undef, jobz == 'N' ? 1 : n, 1) + nzc = -1 + + stemr!( + jobz, + range, + dv, + ev, + vl, + vu, + il, + iu, + w, + Z, + nzc, + isuppz, + work, + lwork, + iwork, + liwork, + ) + + lwork = work[1] + work = Vector{$elty}(undef, lwork) + liwork = iwork[1] + iwork = Vector{BlasInt}(undef, liwork) + nzc = Z[1] + Z = similar(dv, $elty, n, nzc) + + return stemr!( + jobz, + range, + dv, + ev, + vl, + vu, + il, + iu, + w, + Z, + nzc, + isuppz, + work, + lwork, + iwork, + liwork, + ) end - - return wr, wi, Z - end - function lahqr!(H::StridedMatrix{Float64}) - n = size(H, 1) - return lahqr!(false, - false, - H, - 1, - n, - Vector{Float64}(undef, n), - Vector{Float64}(undef, n), - 1, - 1, - Matrix{Float64}(undef, 0, 0)) end +end - ## Cholesky + singular values - function spteqr!( - compz::Char, - d::StridedVector{Float64}, - e::StridedVector{Float64}, - Z::StridedMatrix{Float64}, - work::StridedVector{Float64} = Vector{Float64}(undef, 4length(d))) - - n = length(d) - ldz = stride(Z, 2) - - # Checks - length(e) >= n - 1 || throw(DimensionMismatch("subdiagonal vector is too short")) - if compz == 'N' - elseif compz == in('V', 'I') - size(Z, 1) >= n || throw(DimensionMismatch("Z does not have enough rows")) - size(Z, 2) >= ldz || throw(DimensionMismatch("Z does not have enough columns")) - else - throw(ArgumentError("compz must be either 'N', 'V', or 'I'")) - end +# Non-symmetric eigenvalue problem +function lahqr!( + wantt::Bool, + wantz::Bool, + H::StridedMatrix{Float64}, + ilo::Integer, + ihi::Integer, + wr::StridedVector{Float64}, + wi::StridedVector{Float64}, + iloz::Integer, + ihiz::Integer, + Z::StridedMatrix{Float64}, +) + # SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, + # ILOZ, IHIZ, Z, LDZ, INFO ) + + # .. Scalar Arguments .. + # INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N + # LOGICAL WANTT, WANTZ + # .. + # .. Array Arguments .. + # DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) + + n = LinearAlgebra.checksquare(H) + ldh = stride(H, 2) + ldz = stride(Z, 2) + + info = Ref{BlasInt}(0) + + ccall( + (@blasfunc("dlahqr_"), liblapack_name), + Cvoid, + ( + Ref{BlasInt}, + Ref{BlasInt}, + Ref{BlasInt}, + Ref{BlasInt}, + Ref{BlasInt}, + Ptr{Float64}, + Ref{BlasInt}, + Ptr{Float64}, + Ptr{Float64}, + Ref{BlasInt}, + Ref{BlasInt}, + Ptr{Float64}, + Ref{BlasInt}, + Ptr{BlasInt}, + ), + wantt, + wantz, + n, + ilo, + ihi, + H, + max(1, ldh), + wr, + wi, + iloz, + ihiz, + Z, + max(1, ldz), + info, + ) + + if info[] != 0 + throw(LAPACKException(info[])) + end - info = Vector{BlasInt}(undef, 1) + return wr, wi, Z +end +function lahqr!(H::StridedMatrix{Float64}) + n = size(H, 1) + return lahqr!( + false, + false, + H, + 1, + n, + Vector{Float64}(undef, n), + Vector{Float64}(undef, n), + 1, + 1, + Matrix{Float64}(undef, 0, 0), + ) +end - ccall((@blasfunc(:dpteqr_), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64}, - Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ptr{BlasInt}), - compz, n, d, e, - Z, max(1, ldz), work, info) +## Cholesky + singular values +function spteqr!( + compz::Char, + d::StridedVector{Float64}, + e::StridedVector{Float64}, + Z::StridedMatrix{Float64}, + work::StridedVector{Float64} = Vector{Float64}(undef, 4length(d)), +) + + n = length(d) + ldz = stride(Z, 2) + + # Checks + length(e) >= n - 1 || throw(DimensionMismatch("subdiagonal vector is too short")) + if compz == 'N' + elseif compz == in('V', 'I') + size(Z, 1) >= n || throw(DimensionMismatch("Z does not have enough rows")) + size(Z, 2) >= ldz || throw(DimensionMismatch("Z does not have enough columns")) + else + throw(ArgumentError("compz must be either 'N', 'V', or 'I'")) + end - info[1] == 0 || throw(LAPACKException(info[1])) + info = Vector{BlasInt}(undef, 1) + + ccall( + (@blasfunc(:dpteqr_), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{BlasInt}, + Ptr{Float64}, + Ptr{Float64}, + Ptr{Float64}, + Ref{BlasInt}, + Ptr{Float64}, + Ptr{BlasInt}, + ), + compz, + n, + d, + e, + Z, + max(1, ldz), + work, + info, + ) + + info[1] == 0 || throw(LAPACKException(info[1])) + + d, Z +end - d, Z - end +# Gu's dnc eigensolver +for (f, elty) in ((:dsyevd_, :Float64), (:ssyevd_, :Float32)) - # Gu's dnc eigensolver - for (f, elty) in ((:dsyevd_, :Float64), - (:ssyevd_, :Float32)) - - @eval begin - function syevd!(jobz::Char, uplo::Char, A::StridedMatrix{$elty}) - n = LinearAlgebra.checksquare(A) - lda = stride(A, 2) - w = Vector{$elty}(n) - work = Vector{$elty}(undef, 1) - lwork = BlasInt(-1) - iwork = Vector{BlasInt}(undef, 1) - liwork = BlasInt(-1) - info = BlasInt[0] - for i = 1:2 - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, - Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, - Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}), - jobz, uplo, n, A, - lda, w, work, lwork, - iwork, liwork, info) - - if info[1] != 0 - return LinearAlgebra.LAPACKException(info[1]) - end - if i == 1 - lwork = BlasInt(work[1]) - resize!(work, lwork) - liwork = iwork[1] - resize!(iwork, liwork) - end + @eval begin + function syevd!(jobz::Char, uplo::Char, A::StridedMatrix{$elty}) + n = LinearAlgebra.checksquare(A) + lda = stride(A, 2) + w = Vector{$elty}(n) + work = Vector{$elty}(undef, 1) + lwork = BlasInt(-1) + iwork = Vector{BlasInt}(undef, 1) + liwork = BlasInt(-1) + info = BlasInt[0] + for i = 1:2 + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{$elty}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{BlasInt}, + Ref{BlasInt}, + Ptr{BlasInt}, + ), + jobz, + uplo, + n, + A, + lda, + w, + work, + lwork, + iwork, + liwork, + info, + ) + + if info[1] != 0 + return LinearAlgebra.LAPACKException(info[1]) + end + if i == 1 + lwork = BlasInt(work[1]) + resize!(work, lwork) + liwork = iwork[1] + resize!(iwork, liwork) end - return w, A end + return w, A end end - for (f, elty, relty) in ((:zheevd_, :ComplexF64, :Float64), - (:cheevd_, :ComplexF32, :Float32)) - - @eval begin - function heevd!(jobz::Char, uplo::Char, A::StridedMatrix{$elty}) - n = LinearAlgebra.checksquare(A) - lda = stride(A, 2) - w = Vector{$relty}(n) - work = Vector{$elty}(undef, 1) - lwork = BlasInt(-1) - rwork = Vector{$relty}(undef, 1) - lrwork = BlasInt(-1) - iwork = Vector{BlasInt}(undef, 1) - liwork = BlasInt(-1) - info = BlasInt[0] - for i = 1:2 - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, - Ref{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ref{BlasInt}, - Ptr{$relty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}, - Ptr{BlasInt}), - jobz, uplo, n, A, - lda, w, work, lwork, - rwork, lrwork, iwork, liwork, - info) - - if info[1] != 0 - return LinearAlgebra.LAPACKException(info[1]) - end - - if i == 1 - lwork = BlasInt(work[1]) - resize!(work, lwork) - lrwork = BlasInt(rwork[1]) - resize!(rwork, lrwork) - liwork = iwork[1] - resize!(iwork, liwork) - end +end +for (f, elty, relty) in + ((:zheevd_, :ComplexF64, :Float64), (:cheevd_, :ComplexF32, :Float32)) + + @eval begin + function heevd!(jobz::Char, uplo::Char, A::StridedMatrix{$elty}) + n = LinearAlgebra.checksquare(A) + lda = stride(A, 2) + w = Vector{$relty}(n) + work = Vector{$elty}(undef, 1) + lwork = BlasInt(-1) + rwork = Vector{$relty}(undef, 1) + lrwork = BlasInt(-1) + iwork = Vector{BlasInt}(undef, 1) + liwork = BlasInt(-1) + info = BlasInt[0] + for i = 1:2 + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{$relty}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{$relty}, + Ref{BlasInt}, + Ptr{BlasInt}, + Ref{BlasInt}, + Ptr{BlasInt}, + ), + jobz, + uplo, + n, + A, + lda, + w, + work, + lwork, + rwork, + lrwork, + iwork, + liwork, + info, + ) + + if info[1] != 0 + return LinearAlgebra.LAPACKException(info[1]) + end + + if i == 1 + lwork = BlasInt(work[1]) + resize!(work, lwork) + lrwork = BlasInt(rwork[1]) + resize!(rwork, lrwork) + liwork = iwork[1] + resize!(iwork, liwork) end - return w, A end + return w, A end end +end # Generalized Eigen ## Eigenvectors of (pseudo) triangular matrices -for (f, elty) in ((:dtgevc_, :Float64), - (:stgevc_, :Float32)) +for (f, elty) in ((:dtgevc_, :Float64), (:stgevc_, :Float32)) @eval begin - function tgevc!(side::Char, howmny::Char, select::Vector{BlasInt}, S::StridedMatrix{$elty}, P::StridedMatrix{$elty}, VL::StridedMatrix{$elty}, VR::StridedMatrix{$elty}, work::Vector{$elty}) + function tgevc!( + side::Char, + howmny::Char, + select::Vector{BlasInt}, + S::StridedMatrix{$elty}, + P::StridedMatrix{$elty}, + VL::StridedMatrix{$elty}, + VR::StridedMatrix{$elty}, + work::Vector{$elty}, + ) n = LinearAlgebra.checksquare(S) if LinearAlgebra.checksquare(P) != n @@ -407,29 +625,54 @@ for (f, elty) in ((:dtgevc_, :Float64), if side == 'L' mm = size(VL, 2) if size(VL, 1) != n - throw(DimensionMismatch("first dimension of output is $size(VL, 1) but should be $n")) + throw( + DimensionMismatch( + "first dimension of output is $size(VL, 1) but should be $n", + ), + ) end elseif side == 'R' mm = size(VR, 2) if size(VR, 1) != n - throw(DimensionMismatch("first dimension of output is $size(VR, 1) but should be $n")) + throw( + DimensionMismatch( + "first dimension of output is $size(VR, 1) but should be $n", + ), + ) end elseif side == 'B' mm = size(VL, 2) if mm != size(VR, 2) - throw(DimensionMismatch("the matrices for storing the left and right eigenvectors must have the same number of columns when side == 'B'")) + throw( + DimensionMismatch( + "the matrices for storing the left and right eigenvectors must have the same number of columns when side == 'B'", + ), + ) end if size(VL, 1) != n - throw(DimensionMismatch("first dimension of output is $size(VL, 1) but should be $n")) - end else + throw( + DimensionMismatch( + "first dimension of output is $size(VL, 1) but should be $n", + ), + ) + end + else if size(VR, 1) != n - throw(DimensionMismatch("first dimension of output is $size(VR, 1) but should be $n")) + throw( + DimensionMismatch( + "first dimension of output is $size(VR, 1) but should be $n", + ), + ) end throw(ArgumentError("side must be either 'L', 'R', or 'B'")) end if howmny == 'A' || howmny == 'B' if mm != n - throw(DimensionMismatch("the number of columne in the output matrices are $mm, but should be $n")) + throw( + DimensionMismatch( + "the number of columne in the output matrices are $mm, but should be $n", + ), + ) end elseif howmny == 'S' mx, mn = extrama(select) @@ -437,7 +680,11 @@ for (f, elty) in ((:dtgevc_, :Float64), throw(ArgumentError("the elements of select must be either 0 or 1")) end if sum(howmny) != mm - throw(DimensionMismatch("the number of columns in the output arrays is $mm, but you have selected $(sum(howmny)) vectors")) + throw( + DimensionMismatch( + "the number of columns in the output arrays is $mm, but you have selected $(sum(howmny)) vectors", + ), + ) end else throw(ArgumentError("howmny must be either A, B, or S")) @@ -446,15 +693,44 @@ for (f, elty) in ((:dtgevc_, :Float64), m = BlasInt[0] info = BlasInt[0] - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt}, - Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, - Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, - Ref{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), - side, howmny, select, n, - S, lds, P, ldp, - VL, ldvl, VR, ldvr, - mm, m, work, info) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ptr{BlasInt}, + Ref{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ref{BlasInt}, + Ptr{BlasInt}, + Ptr{$elty}, + Ptr{BlasInt}, + ), + side, + howmny, + select, + n, + S, + lds, + P, + ldp, + VL, + ldvl, + VR, + ldvr, + mm, + m, + work, + info, + ) if info[1] != 0 throw(LAPACKException(info[1])) @@ -463,7 +739,15 @@ for (f, elty) in ((:dtgevc_, :Float64), return VL, VR, m[1] end - function tgevc!(side::Char, howmny::Char, select::Vector{BlasInt}, S::StridedMatrix{$elty}, P::StridedMatrix{$elty}, VL::StridedMatrix{$elty}, VR::StridedMatrix{$elty}) + function tgevc!( + side::Char, + howmny::Char, + select::Vector{BlasInt}, + S::StridedMatrix{$elty}, + P::StridedMatrix{$elty}, + VL::StridedMatrix{$elty}, + VR::StridedMatrix{$elty}, + ) n = LinearAlgebra.checksquare(S) work = Vector{$elty}(6n) @@ -471,7 +755,13 @@ for (f, elty) in ((:dtgevc_, :Float64), return tgevc!(side, howmny, select, S, P, VL, VR, work) end - function tgevc!(side::Char, howmny::Char, select::Vector{BlasInt}, S::StridedMatrix{$elty}, P::StridedMatrix{$elty}) + function tgevc!( + side::Char, + howmny::Char, + select::Vector{BlasInt}, + S::StridedMatrix{$elty}, + P::StridedMatrix{$elty}, + ) # No checks here as they are done in method above n = LinearAlgebra.checksquare(S) if side == 'L' @@ -507,13 +797,23 @@ end # Rectangular full packed format ## Symmetric rank-k operation for matrix in RFP format. -for (f, elty, relty) in ((:dsfrk_, :Float64, :Float64), - (:ssfrk_, :Float32, :Float32), - (:zhfrk_, :ComplexF64, :Float64), - (:chfrk_, :ComplexF32, :Float32)) +for (f, elty, relty) in ( + (:dsfrk_, :Float64, :Float64), + (:ssfrk_, :Float32, :Float32), + (:zhfrk_, :ComplexF64, :Float64), + (:chfrk_, :ComplexF32, :Float32), +) @eval begin - function sfrk!(transr::Char, uplo::Char, trans::Char, alpha::Real, A::StridedMatrix{$elty}, beta::Real, C::StridedVector{$elty}) + function sfrk!( + transr::Char, + uplo::Char, + trans::Char, + alpha::Real, + A::StridedMatrix{$elty}, + beta::Real, + C::StridedVector{$elty}, + ) chkuplo(uplo) chkstride1(A) if trans in ('N', 'n') @@ -523,57 +823,90 @@ for (f, elty, relty) in ((:dsfrk_, :Float64, :Float64), end lda = max(1, stride(A, 2)) - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, - Ref{BlasInt}, Ref{$relty}, Ptr{$elty}, Ref{BlasInt}, - Ref{$relty}, Ptr{$elty}), - transr, uplo, trans, n, - k, alpha, A, lda, - beta, C) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ref{BlasInt}, + Ref{$relty}, + Ptr{$elty}, + Ref{BlasInt}, + Ref{$relty}, + Ptr{$elty}, + ), + transr, + uplo, + trans, + n, + k, + alpha, + A, + lda, + beta, + C, + ) C end end end # Cholesky factorization of a real symmetric positive definite matrix A -for (f, elty) in ((:dpftrf_, :Float64), - (:spftrf_, :Float32), - (:zpftrf_, :ComplexF64), - (:cpftrf_, :ComplexF32)) +for (f, elty) in ( + (:dpftrf_, :Float64), + (:spftrf_, :Float32), + (:zpftrf_, :ComplexF64), + (:cpftrf_, :ComplexF32), +) @eval begin function pftrf!(transr::Char, uplo::Char, A::StridedVector{$elty}) chkuplo(uplo) - n = round(Int,div(sqrt(8length(A)), 2)) + n = round(Int, div(sqrt(8length(A)), 2)) info = Vector{BlasInt}(undef, 1) - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, - Ptr{BlasInt}), - transr, uplo, n, A, - info) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), + transr, + uplo, + n, + A, + info, + ) A end end end # Computes the inverse of a (real) symmetric positive definite matrix A using the Cholesky factorization -for (f, elty) in ((:dpftri_, :Float64), - (:spftri_, :Float32), - (:zpftri_, :ComplexF64), - (:cpftri_, :ComplexF32)) +for (f, elty) in ( + (:dpftri_, :Float64), + (:spftri_, :Float32), + (:zpftri_, :ComplexF64), + (:cpftri_, :ComplexF32), +) @eval begin function pftri!(transr::Char, uplo::Char, A::StridedVector{$elty}) chkuplo(uplo) - n = round(Int,div(sqrt(8length(A)), 2)) + n = round(Int, div(sqrt(8length(A)), 2)) info = Vector{BlasInt}(undef, 1) - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, - Ptr{BlasInt}), - transr, uplo, n, A, - info) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), + transr, + uplo, + n, + A, + info, + ) A end @@ -581,16 +914,23 @@ for (f, elty) in ((:dpftri_, :Float64), end # DPFTRS solves a system of linear equations A*X = B with a symmetric positive definite matrix A using the Cholesky factorization -for (f, elty) in ((:dpftrs_, :Float64), - (:spftrs_, :Float32), - (:zpftrs_, :ComplexF64), - (:cpftrs_, :ComplexF32)) +for (f, elty) in ( + (:dpftrs_, :Float64), + (:spftrs_, :Float32), + (:zpftrs_, :ComplexF64), + (:cpftrs_, :ComplexF32), +) @eval begin - function pftrs!(transr::Char, uplo::Char, A::StridedVector{$elty}, B::StridedVecOrMat{$elty}) + function pftrs!( + transr::Char, + uplo::Char, + A::StridedVector{$elty}, + B::StridedVecOrMat{$elty}, + ) chkuplo(uplo) chkstride1(B) - n = round(Int,div(sqrt(8length(A)), 2)) + n = round(Int, div(sqrt(8length(A)), 2)) if n != size(B, 1) throw(DimensionMismatch("B has first dimension $(size(B,1)) but needs $n")) end @@ -598,11 +938,28 @@ for (f, elty) in ((:dpftrs_, :Float64), ldb = max(1, stride(B, 2)) info = Vector{BlasInt}(undef, 1) - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, - Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}), - transr, uplo, n, nhrs, - A, B, ldb, info) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ref{BlasInt}, + Ptr{$elty}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{BlasInt}, + ), + transr, + uplo, + n, + nhrs, + A, + B, + ldb, + info, + ) B end @@ -610,37 +967,66 @@ for (f, elty) in ((:dpftrs_, :Float64), end # Solves a matrix equation (one operand is a triangular matrix in RFP format) -for (f, elty) in ((:dtfsm_, :Float64), - (:stfsm_, :Float32), - (:ztfsm_, :ComplexF64), - (:ctfsm_, :ComplexF32)) +for (f, elty) in ( + (:dtfsm_, :Float64), + (:stfsm_, :Float32), + (:ztfsm_, :ComplexF64), + (:ctfsm_, :ComplexF32), +) @eval begin - function tfsm!(transr::Char, - side::Char, - uplo::Char, - trans::Char, - diag::Char, - alpha::$elty, - A::StridedVector{$elty}, - B::StridedVecOrMat{$elty}) + function tfsm!( + transr::Char, + side::Char, + uplo::Char, + trans::Char, + diag::Char, + alpha::$elty, + A::StridedVector{$elty}, + B::StridedVecOrMat{$elty}, + ) chkuplo(uplo) chkside(side) chkdiag(diag) chkstride1(B) m, n = size(B, 1), size(B, 2) if round(Int, div(sqrt(8length(A)), 2)) != m - throw(DimensionMismatch("First dimension of B must equal $(round(Int, div(sqrt(8length(A)), 2))), got $m")) + throw( + DimensionMismatch( + "First dimension of B must equal $(round(Int, div(sqrt(8length(A)), 2))), got $m", + ), + ) end ldb = max(1, stride(B, 2)) - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, - Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, - Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}), - transr, side, uplo, trans, - diag, m, n, alpha, - A, B, ldb) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{UInt8}, + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ref{BlasInt}, + Ref{$elty}, + Ptr{$elty}, + Ptr{$elty}, + Ref{BlasInt}, + ), + transr, + side, + uplo, + trans, + diag, + m, + n, + alpha, + A, + B, + ldb, + ) return B end @@ -648,23 +1034,38 @@ for (f, elty) in ((:dtfsm_, :Float64), end # Computes the inverse of a triangular matrix A stored in RFP format. -for (f, elty) in ((:dtftri_, :Float64), - (:stftri_, :Float32), - (:ztftri_, :ComplexF64), - (:ctftri_, :ComplexF32)) +for (f, elty) in ( + (:dtftri_, :Float64), + (:stftri_, :Float32), + (:ztftri_, :ComplexF64), + (:ctftri_, :ComplexF32), +) @eval begin function tftri!(transr::Char, uplo::Char, diag::Char, A::StridedVector{$elty}) chkuplo(uplo) chkdiag(diag) - n = round(Int,div(sqrt(8length(A)), 2)) + n = round(Int, div(sqrt(8length(A)), 2)) info = Vector{BlasInt}(undef, 1) - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, - Ptr{$elty}, Ptr{BlasInt}), - transr, uplo, diag, n, - A, info) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ptr{$elty}, + Ptr{BlasInt}, + ), + transr, + uplo, + diag, + n, + A, + info, + ) A end @@ -672,23 +1073,40 @@ for (f, elty) in ((:dtftri_, :Float64), end # Copies a triangular matrix from the rectangular full packed format (TF) to the standard full format (TR) -for (f, elty) in ((:dtfttr_, :Float64), - (:stfttr_, :Float32), - (:ztfttr_, :ComplexF64), - (:ctfttr_, :ComplexF32)) +for (f, elty) in ( + (:dtfttr_, :Float64), + (:stfttr_, :Float32), + (:ztfttr_, :ComplexF64), + (:ctfttr_, :ComplexF32), +) @eval begin function tfttr!(transr::Char, uplo::Char, Arf::StridedVector{$elty}) chkuplo(uplo) - n = round(Int,div(sqrt(8length(Arf)), 2)) + n = round(Int, div(sqrt(8length(Arf)), 2)) info = Vector{BlasInt}(undef, 1) A = similar(Arf, $elty, n, n) - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}), - transr, uplo, n, Arf, - A, n, info) + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ptr{$elty}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{BlasInt}, + ), + transr, + uplo, + n, + Arf, + A, + n, + info, + ) A end @@ -696,10 +1114,12 @@ for (f, elty) in ((:dtfttr_, :Float64), end # Copies a triangular matrix from the standard full format (TR) to the rectangular full packed format (TF). -for (f, elty) in ((:dtrttf_, :Float64), - (:strttf_, :Float32), - (:ztrttf_, :ComplexF64), - (:ctrttf_, :ComplexF32)) +for (f, elty) in ( + (:dtrttf_, :Float64), + (:strttf_, :Float32), + (:ztrttf_, :ComplexF64), + (:ctrttf_, :ComplexF32), +) @eval begin function trttf!(transr::Char, uplo::Char, A::StridedMatrix{$elty}) @@ -708,13 +1128,28 @@ for (f, elty) in ((:dtrttf_, :Float64), n = size(A, 1) lda = max(1, stride(A, 2)) info = Vector{BlasInt}(undef, 1) - Arf = similar(A, $elty, div(n*(n+1), 2)) - - ccall((@blasfunc($f), liblapack_name), Cvoid, - (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, - Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), - transr, uplo, n, A, - lda, Arf, info) + Arf = similar(A, $elty, div(n * (n + 1), 2)) + + ccall( + (@blasfunc($f), liblapack_name), + Cvoid, + ( + Ref{UInt8}, + Ref{UInt8}, + Ref{BlasInt}, + Ptr{$elty}, + Ref{BlasInt}, + Ptr{$elty}, + Ptr{BlasInt}, + ), + transr, + uplo, + n, + A, + lda, + Arf, + info, + ) Arf end diff --git a/src/qr.jl b/src/qr.jl index b26df03..323a5d0 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -6,8 +6,10 @@ import LinearAlgebra: reflectorApply! struct QR2{T,S<:AbstractMatrix{T},V<:AbstractVector{T}} <: Factorization{T} factors::S τ::V - QR2{T,S,V}(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T,S<:AbstractMatrix,V<:AbstractVector} = - new(factors, τ) + QR2{T,S,V}( + factors::AbstractMatrix{T}, + τ::AbstractVector{T}, + ) where {T,S<:AbstractMatrix,V<:AbstractVector} = new(factors, τ) end QR2(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = QR2{T,typeof(factors)}(factors, τ) @@ -17,18 +19,22 @@ size(F::QR2, i::Integer...) = size(F.factors, i...) @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")) + throw( + DimensionMismatch( + "reflector must have same length as second dimension of matrix", + ), + ) end @inbounds begin - for i in 1:m + for i = 1:m Aiv = A[i, 1] - for j in 2:n - Aiv += A[i, j]*x[j] + for j = 2:n + Aiv += A[i, j] * x[j] end - Aiv = Aiv*τ + Aiv = Aiv * τ A[i, 1] -= Aiv - for j in 2:n - A[i, j] -= Aiv*x[j]' + for j = 2:n + A[i, j] -= Aiv * x[j]' end end end @@ -46,7 +52,7 @@ end # getindex{T}(A::QR2{T}, ::Type{Tuple{:Q}}) = Q{T,typeof(A)}(A) -function getindex(A::QR2{T}, ::Type{Tuple{:R}}) where T +function getindex(A::QR2{T}, ::Type{Tuple{:R}}) where {T} m, n = size(A) if m >= n UpperTriangular(view(A.factors, 1:n, 1:n)) @@ -55,30 +61,32 @@ function getindex(A::QR2{T}, ::Type{Tuple{:R}}) where T end end -function getindex(A::QR2{T}, ::Type{Tuple{:QBlocked}}) where T +function getindex(A::QR2{T}, ::Type{Tuple{:QBlocked}}) where {T} m, n = size(A) - mmn = min(m,n) + mmn = min(m, n) F = A.factors τ = A.τ Tmat = zeros(T, mmn, mmn) for j = 1:mmn - for i = 1:j - 1 - Tmat[i,j] = τ[i]*(F[j,i] + dot(view(F, j + 1:m, i), view(F, j + 1:m, j))) + for i = 1:j-1 + Tmat[i, j] = τ[i] * (F[j, i] + dot(view(F, j+1:m, i), view(F, j+1:m, j))) end end LinearAlgebra.inv!(LinearAlgebra.UnitUpperTriangular(Tmat)) - for j = 1:min(m,n) - Tmat[j,j] = τ[j] - for i = 1:j - 1 - Tmat[i,j] = Tmat[i,j]*τ[j] + for j = 1:min(m, n) + Tmat[j, j] = τ[j] + for i = 1:j-1 + Tmat[i, j] = Tmat[i, j] * τ[j] end end HouseholderBlock{T,typeof(F),Matrix{T}}(F, UpperTriangular(Tmat)) end # qrUnblocked!(A::StridedMatrix) = LinearAlgebra.qrfactUnblocked!(A) -function qrUnblocked!(A::StridedMatrix{T}, - τ::StridedVector{T} = fill(zero(T), min(size(A)...))) where {T} +function qrUnblocked!( + A::StridedMatrix{T}, + τ::StridedVector{T} = fill(zero(T), min(size(A)...)), +) where {T} m, n = size(A) @@ -102,10 +110,12 @@ function qrUnblocked!(A::StridedMatrix{T}, return QR2{T,typeof(A),typeof(τ)}(A, τ) end -function qrBlocked!(A::StridedMatrix{T}, - blocksize::Integer = 12, - τ::StridedVector{T} = fill(zero(T), min(size(A)...)), - work = Matrix{T}(undef, blocksize, size(A, 2))) where T +function qrBlocked!( + A::StridedMatrix{T}, + blocksize::Integer = 12, + τ::StridedVector{T} = fill(zero(T), min(size(A)...)), + work = Matrix{T}(undef, blocksize, size(A, 2)), +) where {T} m, n = size(A) @@ -117,17 +127,19 @@ function qrBlocked!(A::StridedMatrix{T}, # Apply the QR factorization to the trailing columns (if any) and recurse if n > blocksize # Make a view of the trailing columns - A2 = view(A, :, blocksize + 1:n) + A2 = view(A, :, blocksize+1:n) # Apply Q' to the trailing columns - lmul!(F[Tuple{:QBlocked}]', A2, view(work, 1:min(m, blocksize), 1:(n - blocksize))) + lmul!(F[Tuple{:QBlocked}]', A2, view(work, 1:min(m, blocksize), 1:(n-blocksize))) end # Compute QR factorization of trailing block if m > blocksize && n > blocksize - qrBlocked!(view(A, (blocksize + 1):m, (blocksize + 1):n), - blocksize, - view(τ, (blocksize + 1):min(m, n)), - work) + qrBlocked!( + view(A, (blocksize+1):m, (blocksize+1):n), + blocksize, + view(τ, (blocksize+1):min(m, n)), + work, + ) end return QR2{T,typeof(A),typeof(τ)}(A, τ) diff --git a/src/rectfullpacked.jl b/src/rectfullpacked.jl index 80cd1e5..e6f17ae 100644 --- a/src/rectfullpacked.jl +++ b/src/rectfullpacked.jl @@ -14,7 +14,7 @@ end function Base.size(A::HermitianRFP, i::Integer) if i == 1 || i == 2 - return (isqrt(8*length(A.data) + 1) - 1) >> 1 + return (isqrt(8 * length(A.data) + 1) - 1) >> 1 elseif i > 2 return 1 else @@ -30,56 +30,72 @@ function Base.getindex(A::HermitianRFP, i::Integer, j::Integer) n = size(A, 1) n2 = n >> 1 if i < 1 || i > n || j < 0 || j > n - throw(BoundsError(A, (i,j))) + throw(BoundsError(A, (i, j))) end if A.uplo == 'L' if i < j - return conj(A[j,i]) + return conj(A[j, i]) end if j <= n2 + isodd(n) if A.transr == 'N' - return A.data[i + iseven(n) + (j - 1)*(n + iseven(n))] + return A.data[i+iseven(n)+(j-1)*(n+iseven(n))] else - return conj(A.data[(i - 1)*(n + iseven(n)) + j + iseven(n)]) + return conj(A.data[(i-1)*(n+iseven(n))+j+iseven(n)]) end else if A.transr == 'N' - return conj(A.data[(i - n2 - 1)*(n + iseven(n)) + j - n2 - isodd(n)]) + return conj(A.data[(i-n2-1)*(n+iseven(n))+j-n2-isodd(n)]) else - return A.data[i - n2 - isodd(n) + (j - n2 - 1)*(n + iseven(n))] + return A.data[i-n2-isodd(n)+(j-n2-1)*(n+iseven(n))] end end else if i > j - return conj(A[j,i]) + return conj(A[j, i]) end if j > n2 if A.transr == 'N' - return A.data[i + (j - n2 - 1)*(n + iseven(n))] + return A.data[i+(j-n2-1)*(n+iseven(n))] else - return conj(A.data[(i - n2 - 1)*(n + iseven(n)) + j]) + return conj(A.data[(i-n2-1)*(n+iseven(n))+j]) end else if A.transr == 'N' - return conj(A.data[(i - 1)*(n + iseven(n)) + j + n2 + 1]) + return conj(A.data[(i-1)*(n+iseven(n))+j+n2+1]) else - return A.data[i - n2 - isodd(n) + (j - n2 - 1)*(n + iseven(n))] + return A.data[i-n2-isodd(n)+(j-n2-1)*(n+iseven(n))] end end end end -function Ac_mul_A_RFP(A::Matrix{T}, uplo = :U) where T<:BlasFloat +function Ac_mul_A_RFP(A::Matrix{T}, uplo = :U) where {T<:BlasFloat} n = size(A, 2) if uplo == :U - C = LAPACK2.sfrk!('N', 'U', T <: Complex ? 'C' : 'T', 1.0, A, 0.0, Vector{T}(undef, (n*(n + 1)) >> 1)) + C = LAPACK2.sfrk!( + 'N', + 'U', + T <: Complex ? 'C' : 'T', + 1.0, + A, + 0.0, + Vector{T}(undef, (n * (n + 1)) >> 1), + ) return HermitianRFP(C, 'N', 'U') elseif uplo == :L - C = LAPACK2.sfrk!('N', 'L', T <: Complex ? 'C' : 'T', 1.0, A, 0.0, Vector{T}(undef, (n*(n + 1)) >> 1)) - return HermitianRFP(C, 'N', 'L') - else + C = LAPACK2.sfrk!( + 'N', + 'L', + T <: Complex ? 'C' : 'T', + 1.0, + A, + 0.0, + Vector{T}(undef, (n * (n + 1)) >> 1), + ) + return HermitianRFP(C, 'N', 'L') + else throw(ArgumentError("uplo must be either :L or :U")) end end @@ -104,7 +120,7 @@ end function Base.size(A::TriangularRFP, i::Integer) if i == 1 || i == 2 - return (isqrt(8*length(A.data) + 1) - 1) >> 1 + return (isqrt(8 * length(A.data) + 1) - 1) >> 1 elseif i > 2 return 1 else @@ -127,10 +143,11 @@ function Base.Array(A::TriangularRFP) end end -LinearAlgebra.inv!(A::TriangularRFP) = TriangularRFP(LAPACK2.tftri!(A.transr, A.uplo, 'N', A.data), A.transr, A.uplo) -LinearAlgebra.inv(A::TriangularRFP) = LinearAlgebra.inv!(copy(A)) +LinearAlgebra.inv!(A::TriangularRFP) = + TriangularRFP(LAPACK2.tftri!(A.transr, A.uplo, 'N', A.data), A.transr, A.uplo) +LinearAlgebra.inv(A::TriangularRFP) = LinearAlgebra.inv!(copy(A)) -ldiv!(A::TriangularRFP{T}, B::StridedVecOrMat{T}) where T = +ldiv!(A::TriangularRFP{T}, B::StridedVecOrMat{T}) where {T} = LAPACK2.tfsm!(A.transr, 'L', A.uplo, 'N', 'N', one(T), A.data, B) (\)(A::TriangularRFP, B::StridedVecOrMat) = ldiv!(A, copy(B)) @@ -140,16 +157,18 @@ struct CholeskyRFP{T<:BlasFloat} <: Factorization{T} uplo::Char end -LinearAlgebra.cholesky!(A::HermitianRFP{T}) where {T<:BlasFloat} = CholeskyRFP(LAPACK2.pftrf!(A.transr, A.uplo, copy(A.data)), A.transr, A.uplo) +LinearAlgebra.cholesky!(A::HermitianRFP{T}) where {T<:BlasFloat} = + CholeskyRFP(LAPACK2.pftrf!(A.transr, A.uplo, copy(A.data)), A.transr, A.uplo) LinearAlgebra.cholesky(A::HermitianRFP{T}) where {T<:BlasFloat} = cholesky!(copy(A)) LinearAlgebra.factorize(A::HermitianRFP) = cholesky(A) -Base.copy(F::CholeskyRFP{T}) where T = CholeskyRFP{T}(copy(F.data), F.transr, F.uplo) +Base.copy(F::CholeskyRFP{T}) where {T} = CholeskyRFP{T}(copy(F.data), F.transr, F.uplo) # Solve (\)(A::CholeskyRFP, B::StridedVecOrMat) = LAPACK2.pftrs!(A.transr, A.uplo, A.data, copy(B)) -(\)(A::HermitianRFP, B::StridedVecOrMat) = cholesky(A)\B +(\)(A::HermitianRFP, B::StridedVecOrMat) = cholesky(A) \ B -LinearAlgebra.inv!(A::CholeskyRFP) = HermitianRFP(LAPACK2.pftri!(A.transr, A.uplo, A.data), A.transr, A.uplo) -LinearAlgebra.inv(A::CholeskyRFP) = LinearAlgebra.inv!(copy(A)) +LinearAlgebra.inv!(A::CholeskyRFP) = + HermitianRFP(LAPACK2.pftri!(A.transr, A.uplo, A.data), A.transr, A.uplo) +LinearAlgebra.inv(A::CholeskyRFP) = LinearAlgebra.inv!(copy(A)) LinearAlgebra.inv(A::HermitianRFP) = LinearAlgebra.inv!(cholesky(A)) diff --git a/src/svd.jl b/src/svd.jl index 119cdc5..80c3f3a 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -8,60 +8,67 @@ lmul!(G::LinearAlgebra.Givens, ::Nothing) = nothing rmul!(::Nothing, G::LinearAlgebra.Givens) = nothing function svdvals2x2(d1, d2, e) - d1sq = d1*d1 - d2sq = d2*d2 - esq = e*e + d1sq = d1 * d1 + d2sq = d2 * d2 + esq = e * e b = d1sq + d2sq + esq - D = sqrt(abs2((d1 + d2)*(d1 - d2)) + esq*(2*(d1sq + d2sq) + esq)) + D = sqrt(abs2((d1 + d2) * (d1 - d2)) + esq * (2 * (d1sq + d2sq) + esq)) D2 = b + D - λ1 = 2*d1sq*d2sq/D2 - λ2 = D2/2 + λ1 = 2 * d1sq * d2sq / D2 + λ2 = D2 / 2 return minmax(sqrt(λ1), sqrt(λ2)) end -function svdIter!(B::Bidiagonal{T}, n1, n2, shift, U = nothing, Vᴴ = nothing) where T<:Real +function svdIter!( + B::Bidiagonal{T}, + n1, + n2, + shift, + U = nothing, + Vᴴ = nothing, +) where {T<:Real} if B.uplo === 'U' d = B.dv e = B.ev - G, r = givens(d[n1] - abs2(shift)/d[n1], e[n1], n1, n1 + 1) + G, r = givens(d[n1] - abs2(shift) / d[n1], e[n1], n1, n1 + 1) lmul!(G, Vᴴ) - ditmp = d[n1] - ei1 = e[n1] - di = ditmp*G.c + ei1*G.s - ei1 = -ditmp*G.s + ei1*G.c - di1 = d[n1 + 1] - bulge = di1*G.s - di1 *= G.c + ditmp = d[n1] + ei1 = e[n1] + di = ditmp * G.c + ei1 * G.s + ei1 = -ditmp * G.s + ei1 * G.c + di1 = d[n1+1] + bulge = di1 * G.s + di1 *= G.c - for i = n1:n2 - 2 - G, r = givens(di, bulge, i, i + 1) + for i = n1:n2-2 + G, r = givens(di, bulge, i, i + 1) rmul!(U, G') - d[i] = G.c*di + G.s*bulge - ei = G.c*ei1 + G.s*di1 - di1 = -G.s*ei1 + G.c*di1 - ei1 = e[i + 1] - bulge = G.s*ei1 - ei1 *= G.c - - G, r = givens(ei, bulge, i + 1, i + 2) + d[i] = G.c * di + G.s * bulge + ei = G.c * ei1 + G.s * di1 + di1 = -G.s * ei1 + G.c * di1 + ei1 = e[i+1] + bulge = G.s * ei1 + ei1 *= G.c + + G, r = givens(ei, bulge, i + 1, i + 2) lmul!(G, Vᴴ) - e[i] = ei*G.c + bulge*G.s - di = di1*G.c + ei1*G.s - ei1 = -di1*G.s + ei1*G.c - di2 = d[i + 2] - bulge = di2*G.s - di1 = di2*G.c + e[i] = ei * G.c + bulge * G.s + di = di1 * G.c + ei1 * G.s + ei1 = -di1 * G.s + ei1 * G.c + di2 = d[i+2] + bulge = di2 * G.s + di1 = di2 * G.c end - G, r = givens(di, bulge, n2 - 1, n2) + G, r = givens(di, bulge, n2 - 1, n2) rmul!(U, G') - d[n2 - 1] = G.c*di + G.s*bulge - e[n2 - 1] = G.c*ei1 + G.s*di1 - d[n2] = -G.s*ei1 + G.c*di1 + d[n2-1] = G.c * di + G.s * bulge + e[n2-1] = G.c * ei1 + G.s * di1 + d[n2] = -G.s * ei1 + G.c * di1 else throw(ArgumentError("lower bidiagonal version not implemented yet")) @@ -71,7 +78,13 @@ function svdIter!(B::Bidiagonal{T}, n1, n2, shift, U = nothing, Vᴴ = nothing) end # See LAWN 3 -function svdDemmelKahan!(B::Bidiagonal{T}, n1, n2, U = nothing, Vᴴ = nothing) where T<:Real +function svdDemmelKahan!( + B::Bidiagonal{T}, + n1, + n2, + U = nothing, + Vᴴ = nothing, +) where {T<:Real} if B.uplo === 'U' @@ -79,23 +92,23 @@ function svdDemmelKahan!(B::Bidiagonal{T}, n1, n2, U = nothing, Vᴴ = nothing) e = B.ev oldcs = one(T) - G = LinearAlgebra.Givens(1, 2, one(T), zero(T)) - Gold = G + G = LinearAlgebra.Givens(1, 2, one(T), zero(T)) + Gold = G - for i = n1:n2 - 1 - G, r = givens(d[i] * G.c, e[i], i, i + 1) + for i = n1:n2-1 + G, r = givens(d[i] * G.c, e[i], i, i + 1) lmul!(G, Vᴴ) if i != n1 - e[i - 1] = Gold.s * r + e[i-1] = Gold.s * r end - Gold, d[i] = givens(Gold.c * r, d[i + 1] * G.s, i, i + 1) + Gold, d[i] = givens(Gold.c * r, d[i+1] * G.s, i, i + 1) rmul!(U, Gold') end - h = d[n2] * G.c - e[n2 - 1] = h * Gold.s - d[n2] = h * Gold.c + h = d[n2] * G.c + e[n2-1] = h * Gold.s + d[n2] = h * Gold.c else throw(ArgumentError("lower bidiagonal version not implemented yet")) @@ -106,17 +119,17 @@ end # Recurrence to estimate smallest singular value from LAWN3 Lemma 1 function estimate_σ⁻(dv::AbstractVector, ev::AbstractVector, n1::Integer, n2::Integer) - λ = abs(dv[n2]) + λ = abs(dv[n2]) B∞ = λ - for j in (n2 - 1):-1:n1 - λ = abs(dv[j])*(λ/(λ + abs(ev[j]))) + for j = (n2-1):-1:n1 + λ = abs(dv[j]) * (λ / (λ + abs(ev[j]))) B∞ = min(B∞, λ) end - μ = abs(dv[n1]) + μ = abs(dv[n1]) B1 = μ - for j in n1:(n2 - 1) - μ = abs(dv[j + 1])*(μ/(μ + abs(ev[j]))) + for j = n1:(n2-1) + μ = abs(dv[j+1]) * (μ / (μ + abs(ev[j]))) B1 = min(B1, μ) end @@ -125,7 +138,13 @@ end # The actual SVD solver routine # Notice that this routine doesn't adjust the sign and sorts the values -function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T), debug = false) where T<:Real +function __svd!( + B::Bidiagonal{T}, + U = nothing, + Vᴴ = nothing; + tol = 100eps(T), + debug = false, +) where {T<:Real} n = size(B, 1) n1 = 1 @@ -137,7 +156,7 @@ function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T), # See LAWN3 page 6 and 22 σ⁻ = estimate_σ⁻(d, e, n1, n2) fudge = n - thresh = tol*σ⁻ + thresh = tol * σ⁻ if B.uplo === 'U' while true @@ -148,7 +167,7 @@ function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T), # We are done! return nothing else - if abs(e[n2i - 1]) > thresh + if abs(e[n2i-1]) > thresh n2 = n2i break end @@ -156,17 +175,32 @@ function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T), end # Search for largest sub-bidiagonal matrix ending at n2 - for _n1 = (n2 - 1):-1:1 + for _n1 = (n2-1):-1:1 if _n1 == 1 n1 = _n1 break - elseif abs(e[_n1 - 1]) < thresh + elseif abs(e[_n1-1]) < thresh n1 = _n1 break end end - debug && println("n1=", n1, ", n2=", n2, ", d[n1]=", d[n1], ", d[n2]=", d[n2], ", e[n1]=", e[n1], " e[n2-1]=", e[n2-1], " thresh=", thresh) + debug && println( + "n1=", + n1, + ", n2=", + n2, + ", d[n1]=", + d[n1], + ", d[n2]=", + d[n2], + ", e[n1]=", + e[n1], + " e[n2-1]=", + e[n2-1], + " thresh=", + thresh, + ) # See LAWN 3 p 18 and @@ -175,12 +209,12 @@ function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T), # the zero shift algorithm is required to maintain high relative # accuracy σ⁻ = estimate_σ⁻(d, e, n1, n2) - σ⁺ = max(maximum(view(d, n1:n2)), maximum(view(e, n1:(n2 - 1)))) + σ⁺ = max(maximum(view(d, n1:n2)), maximum(view(e, n1:(n2-1)))) - if fudge*tol*σ⁻ <= eps(σ⁺) + if fudge * tol * σ⁻ <= eps(σ⁺) svdDemmelKahan!(B, n1, n2, U, Vᴴ) else - shift = svdvals2x2(d[n2 - 1], d[n2], e[n2 - 1])[1] + shift = svdvals2x2(d[n2-1], d[n2], e[n2-1])[1] if shift^2 < eps(B.dv[n1]^2) # Shift is too small to affect the iteration so just use # zero-shift algorithm anyway @@ -199,7 +233,7 @@ function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T), end end -function _svdvals!(B::Bidiagonal{T}; tol = eps(T), debug = false) where T<:Real +function _svdvals!(B::Bidiagonal{T}; tol = eps(T), debug = false) where {T<:Real} __svd!(B, tol = tol, debug = debug) return sort(abs.(diag(B)), rev = true) end @@ -217,18 +251,18 @@ function _sort_and_adjust!(U, s, Vᴴ) # FixMe! Try to do this without (much) allocation p = sortperm(s, by = abs, rev = true) if p != 1:n - U[:] = U[:,p] - s[:] = s[p] - Vᴴ[:] = Vᴴ[p,:] + U[:] = U[:, p] + s[:] = s[p] + Vᴴ[:] = Vᴴ[p, :] end return nothing end -function _svd!(B::Bidiagonal{T}; tol = eps(T), debug = false) where T<:Real - n = size(B, 1) +function _svd!(B::Bidiagonal{T}; tol = eps(T), debug = false) where {T<:Real} + n = size(B, 1) - U = Matrix{T}(I, n, n) + U = Matrix{T}(I, n, n) Vᴴ = Matrix{T}(I, n, n) __svd!(B, U, Vᴴ, tol = tol, debug = debug) @@ -254,42 +288,52 @@ function bidiagonalize!(A::AbstractMatrix) if m >= n # tall case: upper bidiagonal - for i in 1:min(m, n) + for i = 1:min(m, n) x = view(A, i:m, i) τi = LinearAlgebra.reflector!(x) push!(τl, τi) - LinearAlgebra.reflectorApply!(x, τi, view(A, i:m, (i + 1):n)) + LinearAlgebra.reflectorApply!(x, τi, view(A, i:m, (i+1):n)) if i < n - x = view(A, i, i + 1:n) + x = view(A, i, i+1:n) conj!(x) τi = LinearAlgebra.reflector!(x) push!(τr, τi) - LinearAlgebra.reflectorApply!(view(A, (i + 1):m, (i + 1):n), x, τi) + LinearAlgebra.reflectorApply!(view(A, (i+1):m, (i+1):n), x, τi) end end bd = Bidiagonal(real(diag(A)), real(diag(A, 1)), :U) - return BidiagonalFactorization{eltype(bd),typeof(bd.dv),typeof(A),typeof(τl)}(bd, A, τl, τr) + return BidiagonalFactorization{eltype(bd),typeof(bd.dv),typeof(A),typeof(τl)}( + bd, + A, + τl, + τr, + ) else # wide case: lower bidiagonal - for i in 1:min(m, n) + for i = 1:min(m, n) x = view(A, i, i:n) conj!(x) τi = LinearAlgebra.reflector!(x) push!(τr, τi) - LinearAlgebra.reflectorApply!(view(A, (i + 1):m, i:n), x, τi) + LinearAlgebra.reflectorApply!(view(A, (i+1):m, i:n), x, τi) if i < m - x = view(A, i + 1:m, i) + x = view(A, i+1:m, i) τi = LinearAlgebra.reflector!(x) push!(τl, τi) - LinearAlgebra.reflectorApply!(x, τi, view(A, (i + 1):m, (i + 1):n)) + LinearAlgebra.reflectorApply!(x, τi, view(A, (i+1):m, (i+1):n)) end end bd = Bidiagonal(real(diag(A)), real(diag(A, -1)), :L) - return BidiagonalFactorization{eltype(bd),typeof(bd.dv),typeof(A),typeof(τl)}(bd, A, τl, τr) + return BidiagonalFactorization{eltype(bd),typeof(bd.dv),typeof(A),typeof(τl)}( + bd, + A, + τl, + τr, + ) end end @@ -298,22 +342,40 @@ end # such that multiplication methods are available for free. function Base.getproperty(F::BidiagonalFactorization, s::Symbol) BD = getfield(F, :bidiagonal) - R = getfield(F, :reflectors) + R = getfield(F, :reflectors) τl = getfield(F, :τl) τr = getfield(F, :τr) if BD.uplo === 'U' if s === :leftQ return LinearAlgebra.QRPackedQ(R, τl) elseif s === :rightQ - factors = copy(transpose(R[1:size(R,2),:])) - return LinearAlgebra.HessenbergQ{eltype(factors),typeof(factors),typeof(τr),false}('U', factors, τr) + factors = copy(transpose(R[1:size(R, 2), :])) + return LinearAlgebra.HessenbergQ{ + eltype(factors), + typeof(factors), + typeof(τr), + false, + }( + 'U', + factors, + τr, + ) else return getfield(F, s) end else if s === :leftQ - factors = R[:,1:size(R,1)] - return LinearAlgebra.HessenbergQ{eltype(factors),typeof(factors),typeof(τr),false}('U', factors, τl) + factors = R[:, 1:size(R, 1)] + return LinearAlgebra.HessenbergQ{ + eltype(factors), + typeof(factors), + typeof(τr), + false, + }( + 'U', + factors, + τl, + ) elseif s === :rightQ # return transpose(LinearAlgebra.LQPackedQ(R, τr)) # FixMe! check that this shouldn't be adjoint LinearAlgebra.QRPackedQ(copy(transpose(R)), τr) @@ -332,16 +394,16 @@ function lmul!(Q::LinearAlgebra.HessenbergQ, B::AbstractVecOrMat) throw(DimensionMismatch("")) end - for j in 1:n - for l in length(Q.τ):-1:1 + for j = 1:n + for l = length(Q.τ):-1:1 τl = Q.τ[l] - ṽ = view(Q.factors, (l + 2):m, l) - b = view(B, (l + 2):m, j) - vᴴb = B[l + 1, j] + ṽ = view(Q.factors, (l+2):m, l) + b = view(B, (l+2):m, j) + vᴴb = B[l+1, j] if length(ṽ) > 0 vᴴb += ṽ'b end - B[l + 1, j] -= τl*vᴴb + B[l+1, j] -= τl * vᴴb b .-= ṽ .* τl .* vᴴb end end @@ -358,16 +420,16 @@ function lmul!(adjQ::AdjointQtype{<:Any,<:LinearAlgebra.HessenbergQ}, B::Abstrac throw(DimensionMismatch("")) end - for j in 1:n - for l in 1:length(Q.τ) + for j = 1:n + for l = 1:length(Q.τ) τl = Q.τ[l] - ṽ = view(Q.factors, (l + 2):m, l) - b = view(B, (l + 2):m, j) - vᴴb = B[l + 1, j] + ṽ = view(Q.factors, (l+2):m, l) + b = view(B, (l+2):m, j) + vᴴb = B[l+1, j] if length(ṽ) > 0 vᴴb += ṽ'b end - B[l + 1, j] -= τl'*vᴴb + B[l+1, j] -= τl' * vᴴb b .-= ṽ .* τl' .* vᴴb end end @@ -383,16 +445,16 @@ function rmul!(A::AbstractMatrix, Q::LinearAlgebra.HessenbergQ) throw(DimensionMismatch("")) end - for i in 1:m - for l in 1:(n - 1) + for i = 1:m + for l = 1:(n-1) τl = Q.τ[l] - ṽ = view(Q.factors, (l + 2):n, l) - aᵀ = transpose(view(A, i, (l + 2):n)) - aᵀv = A[i, l + 1] + ṽ = view(Q.factors, (l+2):n, l) + aᵀ = transpose(view(A, i, (l+2):n)) + aᵀv = A[i, l+1] if length(ṽ) > 0 - aᵀv += aᵀ*ṽ + aᵀv += aᵀ * ṽ end - A[i, l + 1] -= aᵀv*τl + A[i, l+1] -= aᵀv * τl aᵀ .-= aᵀv .* τl .* ṽ' end end @@ -408,16 +470,16 @@ function rmul!(A::AbstractMatrix, adjQ::AdjointQtype{<:Any,<:LinearAlgebra.Hesse throw(DimensionMismatch("")) end - for i in 1:m - for l in (n - 1):-1:1 + for i = 1:m + for l = (n-1):-1:1 τl = Q.τ[l] - ṽ = view(Q.factors, (l + 2):n, l) - aᵀ = transpose(view(A, i, (l + 2):n)) - aᵀv = A[i, l + 1] + ṽ = view(Q.factors, (l+2):n, l) + aᵀ = transpose(view(A, i, (l+2):n)) + aᵀv = A[i, l+1] if length(ṽ) > 0 - aᵀv += aᵀ*ṽ + aᵀv += aᵀ * ṽ end - A[i, l + 1] -= aᵀv*τl' + A[i, l+1] -= aᵀv * τl' aᵀ .-= aᵀv .* τl' .* ṽ' end end @@ -433,16 +495,16 @@ function rmul!(A::AbstractMatrix, Q::LinearAlgebra.LQPackedQ) throw(DimensionMismatch("")) end - for i in 1:m - for l in length(Q.τ):-1:1 + for i = 1:m + for l = length(Q.τ):-1:1 τl = Q.τ[l] - ṽ = view(Q.factors, l, (l + 1):n) - aᵀ = transpose(view(A, i, (l + 1):n)) + ṽ = view(Q.factors, l, (l+1):n) + aᵀ = transpose(view(A, i, (l+1):n)) aᵀv = A[i, l] if length(ṽ) > 0 - aᵀv += aᵀ*ṽ + aᵀv += aᵀ * ṽ end - A[i, l] -= aᵀv*τl + A[i, l] -= aᵀv * τl aᵀ .-= aᵀv .* τl .* ṽ' end end @@ -452,7 +514,8 @@ end # Overload LinearAlgebra methods -LinearAlgebra.svdvals!(B::Bidiagonal{T}; tol = eps(T), debug = false) where T<:Real = _svdvals!(B, tol = tol, debug = debug) +LinearAlgebra.svdvals!(B::Bidiagonal{T}; tol = eps(T), debug = false) where {T<:Real} = + _svdvals!(B, tol = tol, debug = debug) """ svdvals!(A [, tol, debug]) @@ -483,13 +546,15 @@ function LinearAlgebra.svdvals!(A::StridedMatrix; tol = eps(real(eltype(A))), de end # FixMe! The full keyword is redundant for Bidiagonal and should be removed from Base -LinearAlgebra.svd!(B::Bidiagonal{T}; +LinearAlgebra.svd!( + B::Bidiagonal{T}; tol = eps(T), full = false, # To avoid breaking on (real(t), imag(t)) - -@testset "General eigen problem with n=$n and element type=$T" for - n in (10, 23, 100), - T in (Float64, Complex{Float64}) - - A = randn(T, n, n) - vGLA = GenericLinearAlgebra._eigvals!(copy(A)) - vLAPACK = eigvals(A) - vBig = eigvals(big.(A)) # not defined in LinearAlgebra so will dispatch to the version in GenericLinearAlgebra - @test sort(vGLA, by = cplxord) ≈ sort(vLAPACK, by = cplxord) - @test sort(vGLA, by = cplxord) ≈ sort(complex(eltype(A)).(vBig), by = cplxord) - @test issorted(vBig, by = cplxord) -end - -@testset "make sure that solver doesn't hang" begin - for i in 1:1000 - A = randn(8, 8) - sort(abs.(GenericLinearAlgebra._eigvals!(copy(A)))) ≈ sort(abs.(eigvals(A))) + cplxord = t -> (real(t), imag(t)) + + @testset "General eigen problem with n=$n and element type=$T" for n in (10, 23, 100), + T in (Float64, Complex{Float64}) + + A = randn(T, n, n) + vGLA = GenericLinearAlgebra._eigvals!(copy(A)) + vLAPACK = eigvals(A) + vBig = eigvals(big.(A)) # not defined in LinearAlgebra so will dispatch to the version in GenericLinearAlgebra + @test sort(vGLA, by = cplxord) ≈ sort(vLAPACK, by = cplxord) + @test sort(vGLA, by = cplxord) ≈ sort(complex(eltype(A)).(vBig), by = cplxord) + @test issorted(vBig, by = cplxord) end -end -@testset "Convergence in corner cases. Issue 29." begin - function H(n::Int) - H = zeros(2n, 2n) - for i = 1 : 2 : 2n - H[i, i+1] = 1 - H[i+1, i] = 1 + @testset "make sure that solver doesn't hang" begin + for i = 1:1000 + A = randn(8, 8) + sort(abs.(GenericLinearAlgebra._eigvals!(copy(A)))) ≈ sort(abs.(eigvals(A))) end - H end - function E(n::Int) - E = zeros(2n, 2n) - for i = 1 : (n - 1) - E[2i + 1, 2i] = 1 + @testset "Convergence in corner cases. Issue 29." begin + function H(n::Int) + H = zeros(2n, 2n) + for i = 1:2:2n + H[i, i+1] = 1 + H[i+1, i] = 1 + end + H end - E[1, 2n] = 1 - E - end - my_matrix(n::Int, η::Float64 = 1e-9) = H(n) .+ η .* E(n) - - A = my_matrix(4, 1e-3); - @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(A), by = t -> (real(t), imag(t))) -end + function E(n::Int) + E = zeros(2n, 2n) + for i = 1:(n-1) + E[2i+1, 2i] = 1 + end + E[1, 2n] = 1 + E + end -@testset "Convergence in with 0s Issue 58." begin - A = [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 0.0] - @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(A), by = t -> (real(t), imag(t))) - B = [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 -1.0 0.0] - @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(B))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(B), by = t -> (real(t), imag(t))) -end + my_matrix(n::Int, η::Float64 = 1e-9) = H(n) .+ η .* E(n) -@testset "Extract Schur factor" begin - A = randn(5, 5) - @test sum(eigvals(schur(A).T)) ≈ sum(eigvals(Float64.(schur(big.(A)).T))) - @test sum(eigvals(schur(A).Schur)) ≈ sum(eigvals(Float64.(schur(big.(A)).Schur))) -end + A = my_matrix(4, 1e-3) + @test sort( + GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), + by = t -> (real(t), imag(t)), + ) ≈ sort(eigvals(A), by = t -> (real(t), imag(t))) + end -@testset "Issue 63" begin - A = [1 2 1 1 1 1 - 0 1 0 1 0 1 - 1 1 2 0 1 0 - 0 1 0 2 1 1 - 1 1 1 0 2 0 - 0 1 1 1 0 2] - - z1 = [complex(1, sqrt(big(3))), complex(1, -sqrt(big(3)))] - z2 = [ - 4*2^(big(2)/3)/complex(5, 3*sqrt(big(111)))^(big(1)/3), - complex(5, 3*sqrt(big(111)))^(big(1)/3)/2^(big(2)/3) - ] + @testset "Convergence in with 0s Issue 58." begin + A = [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 0.0] + @test sort( + GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), + by = t -> (real(t), imag(t)), + ) ≈ sort(eigvals(A), by = t -> (real(t), imag(t))) + B = [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 -1.0 0.0] + @test sort( + GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(B))), + by = t -> (real(t), imag(t)), + ) ≈ sort(eigvals(B), by = t -> (real(t), imag(t))) + end - truevals = real.([ - (7 - transpose(z1)*z2), - 3, - 3, - 3, - (7 - adjoint(z1)*z2), - (7 + [2, 2]'*z2)])/3 + @testset "Extract Schur factor" begin + A = randn(5, 5) + @test sum(eigvals(schur(A).T)) ≈ sum(eigvals(Float64.(schur(big.(A)).T))) + @test sum(eigvals(schur(A).Schur)) ≈ sum(eigvals(Float64.(schur(big.(A)).Schur))) + end - @test eigvals(big.(A)) ≈ truevals -end + @testset "Issue 63" begin + A = [ + 1 2 1 1 1 1 + 0 1 0 1 0 1 + 1 1 2 0 1 0 + 0 1 0 2 1 1 + 1 1 1 0 2 0 + 0 1 1 1 0 2 + ] -Demmel(η) = [0 1 0 0 - 1 0 η 0 - 0 -η 0 1 - 0 0 1 0] + z1 = [complex(1, sqrt(big(3))), complex(1, -sqrt(big(3)))] + z2 = [ + 4 * 2^(big(2) / 3) / complex(5, 3 * sqrt(big(111)))^(big(1) / 3), + complex(5, 3 * sqrt(big(111)))^(big(1) / 3) / 2^(big(2) / 3), + ] -@testset "Demmel matrix" for t in (1e-10, 1e-9, 1e-8) - # See "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it" - A = Demmel(t) - vals = GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(A, maxiter=35)) - @test abs.(vals) ≈ ones(4) -end + truevals = + real.([ + (7 - transpose(z1) * z2), + 3, + 3, + 3, + (7 - adjoint(z1) * z2), + (7 + [2, 2]' * z2), + ]) / 3 + + @test eigvals(big.(A)) ≈ truevals + end -function Hevil2(θ, κ, α, γ) - # Eq (13) and (14) - β = ω = 0.0 - ν = cos(θ)*cos(2γ) + cos(α + β + ω)*sin(2γ)*κ/2 - σ = 1 + κ*sin(2γ)*cos(α + β + ω - θ) + κ^2*sin(γ)^2 - μ = -sin(θ)*cos(2γ) - sin(α + β + ω)*sin(2γ)*κ/2 - ρ = sqrt(σ - ν^2) - - return [ν (cos(2θ) - ν^2)/ρ μ/ρ*(cos(2θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2) (-2*μ*ν - sin(2*θ))/sqrt(ρ^2 - μ^2) - ρ -ν - sin(2*θ)*μ/ρ^2 -μ/ρ^2*(μ*sin(2*θ) + 2*ν*ρ^2)/sqrt(ρ^2 - μ^2) -μ/ρ*(cos(2*θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2) - 0 sin(2*θ)*sqrt(ρ^2 - μ^2)/ρ^2 ν + sin(2*θ)*μ/ρ^2 (cos(2*θ) - ν^2)/ρ - 0 0 ρ -ν] -end + Demmel(η) = [ + 0 1 0 0 + 1 0 η 0 + 0 -η 0 1 + 0 0 1 0 + ] + + @testset "Demmel matrix" for t in (1e-10, 1e-9, 1e-8) + # See "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it" + A = Demmel(t) + vals = GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(A, maxiter = 35)) + @test abs.(vals) ≈ ones(4) + end -@testset "Complicate matrix from Sandia technical report" begin - H = Hevil2(0.111866322512629152, 1.08867072154101741, 0.338146383137297168, -0.313987810419091240) + function Hevil2(θ, κ, α, γ) + # Eq (13) and (14) + β = ω = 0.0 + ν = cos(θ) * cos(2γ) + cos(α + β + ω) * sin(2γ) * κ / 2 + σ = 1 + κ * sin(2γ) * cos(α + β + ω - θ) + κ^2 * sin(γ)^2 + μ = -sin(θ) * cos(2γ) - sin(α + β + ω) * sin(2γ) * κ / 2 + ρ = sqrt(σ - ν^2) + + return [ + ν (cos(2θ)-ν^2)/ρ μ/ρ*(cos(2θ)-ν^2+ρ^2)/sqrt(ρ^2 - μ^2) (-2*μ*ν-sin(2 * θ))/sqrt(ρ^2 - μ^2) + ρ -ν-sin(2 * θ)*μ/ρ^2 -μ/ρ^2*(μ*sin(2 * θ)+2*ν*ρ^2)/sqrt(ρ^2 - μ^2) -μ/ρ*(cos(2 * θ)-ν^2+ρ^2)/sqrt(ρ^2 - μ^2) + 0 sin(2 * θ)*sqrt(ρ^2 - μ^2)/ρ^2 ν+sin(2 * θ)*μ/ρ^2 (cos(2 * θ)-ν^2)/ρ + 0 0 ρ -ν + ] + end - @test Float64.(abs.(eigvals(big.(H)))) ≈ ones(4) -end + @testset "Complicate matrix from Sandia technical report" begin + H = Hevil2( + 0.111866322512629152, + 1.08867072154101741, + 0.338146383137297168, + -0.313987810419091240, + ) -@testset "Issue 67" for (A, λs) in ( - ([1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0] , - [0.0, 0.0, (3 - √big(3)*im)/2, (3 + √big(3)*im)/2, (3 - √big(3)im)/2, (3 + √big(3)im)/2]), - ([1 0 -1 0 0 0; 0 1 1 1 -1 0; 1 0 -1 0 0 0; 1 -1 0 -1 1 1; 0 0 1 0 0 0; -1 0 1 0 0 0] , - zeros(6)), - ([1 -2 -1 1 -1 0; 0 1 0 -1 0 1; -1 1 2 0 1 0; 0 -1 -2 2 -1 -1; 1 0 -1 0 0 0; 0 -1 -1 1 0 0], - [(1 - √big(3)*im)/2, (1 + √big(3)*im)/2, (1 - √big(3)im)/2, (1 + √big(3)im)/2, 2, 2]), - ([2 0 -1 0 0 1; 0 2 0 -1 -1 0; -1 0 1 0 0 -1; 0 -1 0 1 1 0; 0 1 0 -1 0 0; -1 0 1 0 0 0] , - ones(6)), - ([1 0 1 0 -1 0; -2 1 2 -1 1 1; -1 0 -1 0 1 0; 2 1 2 -1 -2 1; 1 0 1 0 -1 0; 1 -1 -2 1 0 -1] , - [-1, -1, 0, 0, 0, 0])) - - @testset "shouldn't need excessive iterations (30*n) in the Float64 case" begin - GenericLinearAlgebra._schur!(float(A)) + @test Float64.(abs.(eigvals(big.(H)))) ≈ ones(4) end - # For BigFloats, many iterations are required for convergence - # Improving this is probably a hard problem - vals = eigvals(big.(A), maxiter=1500) + @testset "Issue 67" for (A, λs) in ( + ( + [ + 1 -2 1 -1 -1 0 + 0 1 0 1 0 1 + 1 -1 2 0 -1 0 + 0 1 0 2 1 1 + 1 0 1 0 0 0 + 0 -1 1 -1 -2 0 + ], + [ + 0.0, + 0.0, + (3 - √big(3) * im) / 2, + (3 + √big(3) * im) / 2, + (3 - √big(3)im) / 2, + (3 + √big(3)im) / 2, + ], + ), + ( + [ + 1 0 -1 0 0 0 + 0 1 1 1 -1 0 + 1 0 -1 0 0 0 + 1 -1 0 -1 1 1 + 0 0 1 0 0 0 + -1 0 1 0 0 0 + ], + zeros(6), + ), + ( + [ + 1 -2 -1 1 -1 0 + 0 1 0 -1 0 1 + -1 1 2 0 1 0 + 0 -1 -2 2 -1 -1 + 1 0 -1 0 0 0 + 0 -1 -1 1 0 0 + ], + [ + (1 - √big(3) * im) / 2, + (1 + √big(3) * im) / 2, + (1 - √big(3)im) / 2, + (1 + √big(3)im) / 2, + 2, + 2, + ], + ), + ( + [ + 2 0 -1 0 0 1 + 0 2 0 -1 -1 0 + -1 0 1 0 0 -1 + 0 -1 0 1 1 0 + 0 1 0 -1 0 0 + -1 0 1 0 0 0 + ], + ones(6), + ), + ( + [ + 1 0 1 0 -1 0 + -2 1 2 -1 1 1 + -1 0 -1 0 1 0 + 2 1 2 -1 -2 1 + 1 0 1 0 -1 0 + 1 -1 -2 1 0 -1 + ], + [-1, -1, 0, 0, 0, 0], + ), + ) + + @testset "shouldn't need excessive iterations (30*n) in the Float64 case" begin + GenericLinearAlgebra._schur!(float(A)) + end - # It's hard to compare complex conjugate pairs so we compare the real and imaginary parts separately - @test sort(real(vals)) ≈ sort(real(λs)) atol=1e-25 - @test sort(imag(vals)) ≈ sort(imag(λs)) atol=1e-25 -end + # For BigFloats, many iterations are required for convergence + # Improving this is probably a hard problem + vals = eigvals(big.(A), maxiter = 1500) -@testset "_hessenberg! and Hessenberg" begin - n = 10 - A = randn(n, n) - HF = GenericLinearAlgebra._hessenberg!(copy(A)) - for i in 1:length(HF.τ) - HM = convert(Matrix, HF.τ[i]) - A[(i + 1):end, :] = HM * A[(i + 1):end, :] - A[:, (i + 1):end] = A[:, (i + 1):end] * HM' + # It's hard to compare complex conjugate pairs so we compare the real and imaginary parts separately + @test sort(real(vals)) ≈ sort(real(λs)) atol = 1e-25 + @test sort(imag(vals)) ≈ sort(imag(λs)) atol = 1e-25 end - @test tril(A, -2) ≈ zeros(n, n) atol = 1e-14 -end -end \ No newline at end of file + @testset "_hessenberg! and Hessenberg" begin + n = 10 + A = randn(n, n) + HF = GenericLinearAlgebra._hessenberg!(copy(A)) + for i = 1:length(HF.τ) + HM = convert(Matrix, HF.τ[i]) + A[(i+1):end, :] = HM * A[(i+1):end, :] + A[:, (i+1):end] = A[:, (i+1):end] * HM' + end + @test tril(A, -2) ≈ zeros(n, n) atol = 1e-14 + end + +end diff --git a/test/eigenselfadjoint.jl b/test/eigenselfadjoint.jl index e44cfd3..4599db2 100644 --- a/test/eigenselfadjoint.jl +++ b/test/eigenselfadjoint.jl @@ -8,10 +8,10 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0 # package since a BigFloat methods isn't defined in # LinearAlgebra T = SymTridiagonal(big.(randn(n)), big.(randn(n - 1))) - vals, vecs = eigen(T) + vals, vecs = eigen(T) @test issorted(vals) @testset "default" begin - @test (vecs'*T)*vecs ≈ Diagonal(vals) + @test (vecs' * T) * vecs ≈ Diagonal(vals) @test eigvals(T) ≈ vals @test vecs'vecs ≈ Matrix(I, n, n) end @@ -20,14 +20,15 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0 vals2, vecs2 = GenericLinearAlgebra.eigen2(T) @test issorted(vals2) @test vals ≈ vals2 - @test vecs[[1,n],:] == vecs2 - @test vecs2*vecs2' ≈ Matrix(I, 2, 2) + @test vecs[[1, n], :] == vecs2 + @test vecs2 * vecs2' ≈ Matrix(I, 2, 2) end @testset "QR version (QL is default)" begin - vals, vecs = GenericLinearAlgebra.eigQR!(copy(T), vectors = Matrix{eltype(T)}(I, n, n)) + vals, vecs = + GenericLinearAlgebra.eigQR!(copy(T), vectors = Matrix{eltype(T)}(I, n, n)) @test issorted(vals) - @test (vecs'*T)*vecs ≈ Diagonal(vals) + @test (vecs' * T) * vecs ≈ Diagonal(vals) @test eigvals(T) ≈ vals @test vecs'vecs ≈ Matrix(I, n, n) end @@ -37,7 +38,7 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0 A = Hermitian(big.(randn(n, n)), uplo) vals, vecs = eigen(A) @testset "default" begin - @test vecs'*A*vecs ≈ diagm(0 => vals) + @test vecs' * A * vecs ≈ diagm(0 => vals) @test eigvals(A) ≈ vals @test vecs'vecs ≈ Matrix(I, n, n) @test issorted(vals) @@ -46,63 +47,64 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0 @testset "eigen2" begin vals2, vecs2 = GenericLinearAlgebra.eigen2(A) @test vals ≈ vals2 - @test vecs[[1,n],:] ≈ vecs2 - @test vecs2*vecs2' ≈ Matrix(I, 2, 2) + @test vecs[[1, n], :] ≈ vecs2 + @test vecs2 * vecs2' ≈ Matrix(I, 2, 2) @test issorted(vals2) end end @testset "(full) Quaternion Hermitian using :$uplo" for uplo in (:L, :U) - V = qr([Quaternion(randn(4)...) for i in 1:n, j in 1:n]).Q - λ = 10 .^ range(-8, stop=0, length=n) - A = Hermitian(V*Diagonal(λ)*V' |> t -> (t + t')/2, uplo) + V = qr([Quaternion(randn(4)...) for i = 1:n, j = 1:n]).Q + λ = 10 .^ range(-8, stop = 0, length = n) + A = Hermitian(V * Diagonal(λ) * V' |> t -> (t + t') / 2, uplo) vals, vecs = eigen(A) @test issorted(vals) @testset "default" begin if uplo == :L # FixMe! Probably an conjugation is off somewhere. Don't have time to check now. - @test_broken vecs'*A*vecs ≈ diagm(0=> vals) + @test_broken vecs' * A * vecs ≈ diagm(0 => vals) else - @test vecs'*A*vecs ≈ diagm(0 => vals) + @test vecs' * A * vecs ≈ diagm(0 => vals) end - @test eigvals(A) ≈ vals - @test vals ≈ λ rtol=1e-13*n - @test vecs'vecs ≈ Matrix(I, n, n) + @test eigvals(A) ≈ vals + @test vals ≈ λ rtol = 1e-13 * n + @test vecs'vecs ≈ Matrix(I, n, n) end @testset "eigen2" begin vals2, vecs2 = GenericLinearAlgebra.eigen2(A) @test issorted(vals2) @test vals ≈ vals2 - @test vecs[[1,n],:] ≈ vecs2 - @test vecs2*vecs2' ≈ Matrix(I, 2, 2) + @test vecs[[1, n], :] ≈ vecs2 + @test vecs2 * vecs2' ≈ Matrix(I, 2, 2) end end @testset "big Hermitian{<:Complex}" begin # This one used to cause an ambiguity error. See #35 - A = complex.(randn(4,4), randn(4,4)) + A = complex.(randn(4, 4), randn(4, 4)) @test Float64.(eigen(Hermitian(big.(A))).values) ≈ eigen(Hermitian(copy(A))).values - @test Float64.(eigvals(Hermitian(big.(A)))) ≈ eigvals(Hermitian(copy(A))) + @test Float64.(eigvals(Hermitian(big.(A)))) ≈ eigvals(Hermitian(copy(A))) end @testset "generic Givens" begin x, y = randn(2) c, s, r = invoke(LinearAlgebra.givensAlgorithm, Tuple{Real,Real}, x, y) - @test c*x + s*y ≈ r - @test c*y ≈ s*x + @test c * x + s * y ≈ r + @test c * y ≈ s * x end @testset "out-of-bounds issue in 1x1 case" begin @test GenericLinearAlgebra._eigvals!(SymTridiagonal([1.0], Float64[])) == [1.0] @test GenericLinearAlgebra._eigen!(SymTridiagonal([1.0], Float64[])).values == [1.0] - @test GenericLinearAlgebra._eigen!(SymTridiagonal([1.0], Float64[])).vectors == fill(1.0, 1, 1) + @test GenericLinearAlgebra._eigen!(SymTridiagonal([1.0], Float64[])).vectors == + fill(1.0, 1, 1) end # Issue #52 @testset "convergence criterion when diagonal has zeros" begin M1 = Hermitian(zeros(3, 3)) - M1[1,1] = 0.01 + M1[1, 1] = 0.01 M2 = Hermitian(zeros(3, 3)) M2[3, 3] = 0.01 @test eigvals(M1) == GenericLinearAlgebra._eigvals!(copy(M1)) diff --git a/test/juliaBLAS.jl b/test/juliaBLAS.jl index 8ee4157..1664304 100644 --- a/test/juliaBLAS.jl +++ b/test/juliaBLAS.jl @@ -7,14 +7,14 @@ using Test, GenericLinearAlgebra, LinearAlgebra complex.(randn(5, 2), randn(5, 2)), complex.(randn(5), randn(5)), ) - @test rankUpdate!(copy(A), x) ≈ A .+ x.*x' - @test rankUpdate!(copy(Ac), xc) ≈ Ac .+ xc.*xc' + @test rankUpdate!(copy(A), x) ≈ A .+ x .* x' + @test rankUpdate!(copy(Ac), xc) ≈ Ac .+ xc .* xc' - @test rankUpdate!(copy(A), B, 0.5, 0.5) ≈ 0.5*A + 0.5*B*B' - @test rankUpdate!(copy(Ac), Bc, 0.5, 0.5) ≈ 0.5*Ac + 0.5*Bc*Bc' + @test rankUpdate!(copy(A), B, 0.5, 0.5) ≈ 0.5 * A + 0.5 * B * B' + @test rankUpdate!(copy(Ac), Bc, 0.5, 0.5) ≈ 0.5 * Ac + 0.5 * Bc * Bc' @test invoke(rankUpdate!, Tuple{Hermitian,StridedVecOrMat,Real}, copy(Ac), Bc, 1.0) ≈ - rankUpdate!(copy(Ac), Bc, 1.0, 1.0) + rankUpdate!(copy(Ac), Bc, 1.0, 1.0) end @testset "triangular multiplication: $(typeof(T))" for T ∈ ( @@ -24,6 +24,6 @@ end UnitLowerTriangular(complex.(randn(5, 5), randn(5, 5))), ) B = complex.(randn(5, 5), randn(5, 5)) - @test lmul!(T, copy(B), complex(0.5, 0.5)) ≈ T*B*complex(0.5, 0.5) - @test lmul!(T', copy(B), complex(0.5, 0.5)) ≈ T'*B*complex(0.5, 0.5) -end \ No newline at end of file + @test lmul!(T, copy(B), complex(0.5, 0.5)) ≈ T * B * complex(0.5, 0.5) + @test lmul!(T', copy(B), complex(0.5, 0.5)) ≈ T' * B * complex(0.5, 0.5) +end diff --git a/test/lapack.jl b/test/lapack.jl index 764cdd7..bc6707e 100644 --- a/test/lapack.jl +++ b/test/lapack.jl @@ -33,7 +33,7 @@ using GenericLinearAlgebra.LAPACK2 @test vals ≈ _vals @test abs.(_vecs'vecs) ≈ Matrix(I, n, n) - _vals, _vecs = LAPACK2.stedc!('V', copy(T.dv), copy(T.ev), Matrix{Float64}(I, n,n)) + _vals, _vecs = LAPACK2.stedc!('V', copy(T.dv), copy(T.ev), Matrix{Float64}(I, n, n)) @test vals ≈ _vals @test abs.(_vecs'vecs) ≈ Matrix(I, n, n) end @@ -68,4 +68,4 @@ using GenericLinearAlgebra.LAPACK2 # LAPACK's multishift algorithm (the default) seems to be broken @test !(_vals ≈ sort(eigvals(T))) end -end \ No newline at end of file +end diff --git a/test/qr.jl b/test/qr.jl index c66d1fc..736ab79 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -4,24 +4,33 @@ using GenericLinearAlgebra using GenericLinearAlgebra: qrBlocked! @testset "The QR decomposition" begin - @testset "Problem dimension ($m,$n) with block size $bz" for - (m, n) in (( 10, 5), ( 10, 10), ( 5, 10), - (100, 50), (100, 100), (50, 100)), - bz in (1, 2, 3, 4, 7, 8, 9, 15, 16, 17, 31, 32, 33) + @testset "Problem dimension ($m,$n) with block size $bz" for (m, n) in ( + (10, 5), + (10, 10), + (5, 10), + (100, 50), + (100, 100), + (50, 100), + ), + bz in (1, 2, 3, 4, 7, 8, 9, 15, 16, 17, 31, 32, 33) A = randn(m, n) Aqr = qrBlocked!(copy(A), bz) AqrQ = Aqr[Tuple{:QBlocked}] if m >= n - @test (AqrQ'A)[1:min(m,n),:] ≈ Aqr[Tuple{:R}] + @test (AqrQ'A)[1:min(m, n), :] ≈ Aqr[Tuple{:R}] else # For type stbility getindex(,Tuple{:R}) throw when the output is trapezoid @test (AqrQ'A) ≈ triu(Aqr.factors) end - @test AqrQ'*(AqrQ*A) ≈ A + @test AqrQ' * (AqrQ * A) ≈ A end @testset "Error paths" begin - @test_throws DimensionMismatch LinearAlgebra.reflectorApply!(zeros(5, 5), zeros(4), 1.0) - @test_throws ArgumentError qrBlocked!(randn(5, 10))[Tuple{:R}] + @test_throws DimensionMismatch LinearAlgebra.reflectorApply!( + zeros(5, 5), + zeros(4), + 1.0, + ) + @test_throws ArgumentError qrBlocked!(randn(5, 10))[Tuple{:R}] end end diff --git a/test/rectfullpacked.jl b/test/rectfullpacked.jl index 37e0b9a..7546d57 100644 --- a/test/rectfullpacked.jl +++ b/test/rectfullpacked.jl @@ -4,8 +4,7 @@ import GenericLinearAlgebra: Ac_mul_A_RFP, TriangularRFP @testset "Rectuangular Full Pack Format" begin - @testset "Core generic functionality" for n in (6, 7), - uplo in (:U, :L) + @testset "Core generic functionality" for n in (6, 7), uplo in (:U, :L) A = rand(10, n) @@ -17,11 +16,11 @@ import GenericLinearAlgebra: Ac_mul_A_RFP, TriangularRFP @test size(AcA_RFP, 3) == 1 @test_throws BoundsError AcA_RFP[0, 1] @test_throws BoundsError AcA_RFP[1, 0] - @test_throws BoundsError AcA_RFP[n + 1, 1] - @test_throws BoundsError AcA_RFP[1, n + 1] + @test_throws BoundsError AcA_RFP[n+1, 1] + @test_throws BoundsError AcA_RFP[1, n+1] - @test AcA_RFP[2, 1] == AcA_RFP[1, 2] - @test AcA_RFP[end - 2, end - 1] == AcA_RFP[end - 1, end - 2] + @test AcA_RFP[2, 1] == AcA_RFP[1, 2] + @test AcA_RFP[end-2, end-1] == AcA_RFP[end-1, end-2] end @testset "Triangular" begin @@ -36,40 +35,48 @@ import GenericLinearAlgebra: Ac_mul_A_RFP, TriangularRFP # @test_throws BoundsError Atr_RFP[n + 1, 1] # @test_throws BoundsError Atr_RFP[1, n + 1] - @test_broken Atr_RFP[2, 1] == Atr_RFP[1, 2] - @test_broken Atr_RFP[end - 2, end - 1] == Atr_RFP[end - 1, end - 2] + @test_broken Atr_RFP[2, 1] == Atr_RFP[1, 2] + @test_broken Atr_RFP[end-2, end-1] == Atr_RFP[end-1, end-2] end end - @testset "Hermitian with element type: $elty. Problem size: $n" for - elty in (Float32, Float64, Complex{Float32}, Complex{Float64}), - n in (6, 7), - uplo in (:L, :U) + @testset "Hermitian with element type: $elty. Problem size: $n" for elty in ( + Float32, + Float64, + Complex{Float32}, + Complex{Float64}, + ), + n in (6, 7), + uplo in (:L, :U) - A = rand(elty, 10, n) - AcA = A'A + A = rand(elty, 10, n) + AcA = A'A AcA_RFP = Ac_mul_A_RFP(A, uplo) - o = ones(elty, n) + o = ones(elty, n) - @test AcA ≈ A'A - @test AcA\o ≈ AcA_RFP\o + @test AcA ≈ A'A + @test AcA \ o ≈ AcA_RFP \ o @test inv(AcA) ≈ inv(AcA_RFP) @test inv(cholesky(AcA)) ≈ inv(factorize(AcA_RFP)) end - @testset "Hermitian with element type: $elty. Problem size: $n" for - elty in (Float32, Float64, Complex{Float32}, Complex{Float64}), - n in (6, 7), - uplo in (:L, :U) + @testset "Hermitian with element type: $elty. Problem size: $n" for elty in ( + Float32, + Float64, + Complex{Float32}, + Complex{Float64}, + ), + n in (6, 7), + uplo in (:L, :U) - A = lu(rand(elty, n, n)).U - A = uplo == :U ? A : copy(A') + A = lu(rand(elty, n, n)).U + A = uplo == :U ? A : copy(A') A_RFP = TriangularRFP(A, uplo) - o = ones(elty, n) + o = ones(elty, n) @test_broken A ≈ A_RFP - @test A ≈ Array(A_RFP) - @test A\o ≈ A_RFP\o - @test inv(A) ≈ Array(inv(A_RFP)) + @test A ≈ Array(A_RFP) + @test A \ o ≈ A_RFP \ o + @test inv(A) ≈ Array(inv(A_RFP)) end end diff --git a/test/runtests.jl b/test/runtests.jl index a81e88a..ecc35c7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,13 +1,13 @@ using Test # @testset "The LinearAlgebra Test Suite" begin - include("juliaBLAS.jl") - include("cholesky.jl") - include("qr.jl") - include("eigenselfadjoint.jl") - include("eigengeneral.jl") - include("tridiag.jl") - include("svd.jl") - include("rectfullpacked.jl") - include("lapack.jl") +include("juliaBLAS.jl") +include("cholesky.jl") +include("qr.jl") +include("eigenselfadjoint.jl") +include("eigengeneral.jl") +include("tridiag.jl") +include("svd.jl") +include("rectfullpacked.jl") +include("lapack.jl") # end diff --git a/test/svd.jl b/test/svd.jl index 6652ca0..aff1e9a 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -1,53 +1,61 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats @testset "Singular value decomposition" begin - @testset "Problem dimension ($m,$n)" for - (m,n) in ((6,5) , (6,6) , (5,6), - (60, 50) , (60, 60) , (50, 60), - (200, 150), (200, 200), (150, 200)) - - vals = reverse(collect(1:min(m,n))) - U = qr(Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:m, j = 1:min(m,n)]).Q - V = qr(Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:n, j = 1:min(m,n)]).Q + @testset "Problem dimension ($m,$n)" for (m, n) in ( + (6, 5), + (6, 6), + (5, 6), + (60, 50), + (60, 60), + (50, 60), + (200, 150), + (200, 200), + (150, 200), + ) + + vals = reverse(collect(1:min(m, n))) + U = qr(Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:m, j = 1:min(m, n)]).Q + V = qr(Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:n, j = 1:min(m, n)]).Q # FixMe! Using Array here shouldn't be necessary. Can be removed once # the bug in LinearAlgebra is fixed - A = U*Array(Diagonal(vals))*V' + A = U * Array(Diagonal(vals)) * V' @test size(A) == (m, n) @test vals ≈ svdvals(A) F = svd(A) @test vals ≈ F.S - @show norm(F.U'*A*F.V - Diagonal(F.S), Inf) - @test F.U'*A*F.V ≈ Diagonal(F.S) + @show norm(F.U' * A * F.V - Diagonal(F.S), Inf) + @test F.U' * A * F.V ≈ Diagonal(F.S) end @testset "The Ivan Slapničar Challenge" begin # This matrix used to hang (for n = 70). Thanks to Ivan Slapničar for reporting. n = 70 - J = Bidiagonal(0.5 * ones(n), ones(n-1), :U) + J = Bidiagonal(0.5 * ones(n), ones(n - 1), :U) @test GenericLinearAlgebra._svdvals!(copy(J)) ≈ svdvals(J) - @test GenericLinearAlgebra._svdvals!(copy(J))[end] / svdvals(J)[end] - 1 < n*eps() + @test GenericLinearAlgebra._svdvals!(copy(J))[end] / svdvals(J)[end] - 1 < n * eps() end - @testset "Compare to Base methods. Problem dimension ($m,$n)" for - (m, n) in ((10, 9), # tall - (10, 10), # square - (9 , 10)) # wide + @testset "Compare to Base methods. Problem dimension ($m,$n)" for (m, n) in ( + (10, 9), # tall + (10, 10), # square + (9, 10), + ) # wide - A = randn(m,n) + A = randn(m, n) Abig = big.(A) @test svdvals(A) ≈ Vector{Float64}(svdvals(Abig)) - @test cond(A) ≈ Float64(cond(Abig)) + @test cond(A) ≈ Float64(cond(Abig)) - F = svd(A) + F = svd(A) Fbig = svd(Abig) @test abs.(F.U'Float64.(Fbig.U)) ≈ I @test abs.(F.V'Float64.(Fbig.V)) ≈ I - F = svd(A, full=true) - Fbig = svd(Abig, full=true) + F = svd(A, full = true) + Fbig = svd(Abig, full = true) @test abs.(F.U'Float64.(Fbig.U)) ≈ I @test abs.(F.V'Float64.(Fbig.V)) ≈ I end @@ -57,7 +65,7 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats A = U0[:, 1:3] * V0[:, 1:3]' U, S, V = svd(A) - @test A ≈ U*Diagonal(S)*V' + @test A ≈ U * Diagonal(S) * V' end @testset "Very small matrices. Issue 79" begin @@ -66,10 +74,10 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats FAb = svd(big.(A)) FAtb = svd(big.(A')) @test FA.S ≈ Float64.(FAb.S) ≈ Float64.(FAtb.S) - @test abs.(FA.U'*Float64.(FAb.U)) ≈ I - @test abs.(FA.U'*Float64.(FAtb.V)) ≈ I - @test abs.(FA.V'*Float64.(FAb.V)) ≈ I - @test abs.(FA.V'*Float64.(FAtb.U)) ≈ I + @test abs.(FA.U' * Float64.(FAb.U)) ≈ I + @test abs.(FA.U' * Float64.(FAtb.V)) ≈ I + @test abs.(FA.V' * Float64.(FAb.V)) ≈ I + @test abs.(FA.V' * Float64.(FAtb.U)) ≈ I end @testset "Issue 81" begin @@ -77,33 +85,33 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats @test Float64.(svdvals(big.(A))) ≈ svdvals(A) A = [ - 0.3 0.0 0.0 0.0 0.0 0.2 0.3 0.0; - 0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0; - 0.0 -0.2 0.0 0.0 0.0 0.0 0.0 -0.2; - 0.3 0.0 0.0 0.0 0.0 0.2 0.4 0.0; - 0.0 0.4 -0.2 0.0 0.0 0.0 0.0 0.3; - 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.0; - 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0; - 0.0 0.3 -0.2 0.0 0.0 0.0 0.0 0.3 + 0.3 0.0 0.0 0.0 0.0 0.2 0.3 0.0 + 0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0 + 0.0 -0.2 0.0 0.0 0.0 0.0 0.0 -0.2 + 0.3 0.0 0.0 0.0 0.0 0.2 0.4 0.0 + 0.0 0.4 -0.2 0.0 0.0 0.0 0.0 0.3 + 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.0 + 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0 + 0.0 0.3 -0.2 0.0 0.0 0.0 0.0 0.3 ] @test GenericLinearAlgebra._svdvals!( - GenericLinearAlgebra.bidiagonalize!(copy(A)).bidiagonal + GenericLinearAlgebra.bidiagonalize!(copy(A)).bidiagonal, ) ≈ svdvals(A) n = 17 A = zeros(Double64, n, n) - for j in 1:n, i in 1:n + for j = 1:n, i = 1:n A[i, j] = 1 / Double64(i + j - 1) end @test svdvals(A) ≈ svdvals(Float64.(A)) # From https://github.com/JuliaMath/DoubleFloats.jl/issues/149 n = 64 - c = Complex{BigFloat}(3//1 + 1im//1) + c = Complex{BigFloat}(3 // 1 + 1im // 1) A = diagm( - 1 => c*ones(BigFloat, n - 1), - -1 => c*ones(BigFloat,n - 1), - -2 => ones(BigFloat, n - 2) + 1 => c * ones(BigFloat, n - 1), + -1 => c * ones(BigFloat, n - 1), + -2 => ones(BigFloat, n - 2), ) @test svdvals(A) ≈ svdvals(Complex{Double64}.(A)) end diff --git a/test/tridiag.jl b/test/tridiag.jl index 572b4a8..63c3c84 100644 --- a/test/tridiag.jl +++ b/test/tridiag.jl @@ -2,7 +2,6 @@ using GenericLinearAlgebra @testset "Test sign of eigenvalues" begin n = 20 - T = SymTridiagonal(randn(n), randn(n-1)) - @test numnegevals(T) == count(x->x<0, eigvals(T)) + T = SymTridiagonal(randn(n), randn(n - 1)) + @test numnegevals(T) == count(x -> x < 0, eigvals(T)) end - From a89dada0b3335013b353b4ecf4a31da7654cfce9 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Fri, 24 Feb 2023 13:26:00 +0100 Subject: [PATCH 4/4] Set version to 0.3.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a2b58e0..d87dcb7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "GenericLinearAlgebra" uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a" -version = "0.3.6-DEV" +version = "0.3.6" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"