From 66ba2efea2d420508e65f86408a59008fa1ba8e4 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Fri, 18 Apr 2025 17:42:33 +0200
Subject: [PATCH 1/2] Stop pirating eigen and svd functions from LinearAlgebra
Instead, define separate unexported functions in the GenericLinearAlgebra
module. Hence, users would have to explicitly opt in to using
GenericLinearAlgebra. A consequence is that the implementations
in this module will be used even for the element types support by
LAPACK.
Also, delete unused functionality.
---
Project.toml | 4 +-
src/GenericLinearAlgebra.jl | 10 +-
src/cholesky.jl | 2 -
src/eigenGeneral.jl | 43 ++--
src/eigenSelfAdjoint.jl | 398 ++++++++++++++++++++++--------------
src/householder.jl | 19 +-
src/juliaBLAS.jl | 7 -
src/lapack.jl | 12 +-
src/qr.jl | 13 +-
src/svd.jl | 58 ++++--
src/tridiag.jl | 36 ----
test/cholesky.jl | 24 +--
test/eigengeneral.jl | 45 ++--
test/eigenselfadjoint.jl | 91 +++++----
test/juliaBLAS.jl | 31 +--
test/lapack.jl | 21 +-
test/qr.jl | 11 +-
test/runtests.jl | 1 -
test/svd.jl | 115 +++++------
test/tridiag.jl | 7 -
20 files changed, 499 insertions(+), 449 deletions(-)
delete mode 100644 src/tridiag.jl
delete mode 100644 test/tridiag.jl
diff --git a/Project.toml b/Project.toml
index e65e4ae..ce0cbfc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
name = "GenericLinearAlgebra"
uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
-version = "0.3.15"
+version = "1.0.0"
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -10,7 +10,7 @@ libblastrampoline_jll = "8e850b90-86db-534c-a0d3-1478176c7d93"
[compat]
Quaternions = "0.7.0"
-julia = "1.6"
+julia = "1.10"
[extras]
DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
diff --git a/src/GenericLinearAlgebra.jl b/src/GenericLinearAlgebra.jl
index c39aa32..c27ac3c 100644
--- a/src/GenericLinearAlgebra.jl
+++ b/src/GenericLinearAlgebra.jl
@@ -1,6 +1,12 @@
module GenericLinearAlgebra
-import LinearAlgebra: mul!, ldiv!
+using LinearAlgebra: LinearAlgebra,
+ Adjoint, Bidiagonal, Diagonal, Factorization, Givens, HermOrSym, Hermitian, I, LowerTriangular,
+ Rotation, SVD, SymTridiagonal, Symmetric, UnitLowerTriangular, UnitUpperTriangular,
+ UpperTriangular,
+ BLAS,
+ abs2, axpy!, diag, dot, eigencopy_oftype, givens, ishermitian, mul!, rdiv!, tril, triu
+using LinearAlgebra.BLAS: BlasFloat, BlasReal
include("juliaBLAS.jl")
include("lapack.jl")
@@ -9,6 +15,6 @@ include("householder.jl")
include("qr.jl")
include("eigenSelfAdjoint.jl")
include("eigenGeneral.jl")
-include("tridiag.jl")
include("svd.jl")
+
end
diff --git a/src/cholesky.jl b/src/cholesky.jl
index 9029fcf..3073146 100644
--- a/src/cholesky.jl
+++ b/src/cholesky.jl
@@ -1,5 +1,3 @@
-using LinearAlgebra: rdiv!
-
function cholUnblocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}) where {T<:Number}
n = LinearAlgebra.checksquare(A)
A[1, 1] = sqrt(A[1, 1])
diff --git a/src/eigenGeneral.jl b/src/eigenGeneral.jl
index 893800b..58ce203 100644
--- a/src/eigenGeneral.jl
+++ b/src/eigenGeneral.jl
@@ -1,9 +1,3 @@
-using Printf
-using LinearAlgebra
-using LinearAlgebra: Givens, Rotation, givens
-
-import Base: \
-
# Hessenberg Matrix
struct HessenbergMatrix{T,S<:StridedMatrix} <: AbstractMatrix{T}
data::S
@@ -29,14 +23,6 @@ function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat)
end
(\)(H::HessenbergMatrix, B::AbstractVecOrMat) = ldiv!(copy(H), copy(B))
-if VERSION < v"1.10"
- # ensure tests pass on Julia v1.6
- copy_similar(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T, size(A)), A)
- eigtype(T) = promote_type(Float32, typeof(zero(T)/sqrt(abs2(one(T)))))
- eigencopy_oftype(A, S) = copy_similar(A, S)
- LinearAlgebra.eigvals(A::HessenbergMatrix{T}; kws...) where T = LinearAlgebra.eigvals!(eigencopy_oftype(A, eigtype(T)); kws...)
-end
-
# Hessenberg factorization
struct HessenbergFactorization{T,S<:StridedMatrix,U} <: Factorization{T}
data::S
@@ -59,7 +45,7 @@ function _hessenberg!(A::StridedMatrix{T}) where {T}
end
return HessenbergFactorization{T,typeof(A),eltype(τ)}(A, τ)
end
-LinearAlgebra.hessenberg!(A::StridedMatrix) = _hessenberg!(A)
+hessenberg!(A::StridedMatrix) = _hessenberg!(A)
Base.size(H::HessenbergFactorization, args...) = size(H.data, args...)
@@ -179,10 +165,7 @@ function _schur!(
return Schur{T,typeof(HH)}(HH, τ)
end
-_schur!(A::StridedMatrix; kwargs...) = _schur!(_hessenberg!(A); kwargs...)
-
-# FIXME! Move this method to piracy extension
-LinearAlgebra.schur!(A::StridedMatrix; kwargs...) = _schur!(A; kwargs...)
+schur!(A::StridedMatrix; kwargs...) = _schur!(_hessenberg!(A); kwargs...)
function singleShiftQR!(
HH::StridedMatrix,
@@ -278,11 +261,11 @@ function doubleShiftQR!(
return HH
end
-_eigvals!(A::StridedMatrix; kwargs...) = _eigvals!(_schur!(A; 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!(
+function eigvals!(
A::StridedMatrix;
sortby::Union{Function,Nothing} = LinearAlgebra.eigsortby,
kwargs...,
@@ -294,13 +277,13 @@ function LinearAlgebra.eigvals!(
LinearAlgebra.sorteig!(_eigvals!(A; kwargs...), sortby)
end
-LinearAlgebra.eigvals!(
+eigvals!(
H::HessenbergMatrix;
sortby::Union{Function,Nothing} = LinearAlgebra.eigsortby,
kwargs...,
) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby)
-LinearAlgebra.eigvals!(
+eigvals!(
H::HessenbergFactorization;
sortby::Union{Function,Nothing} = LinearAlgebra.eigsortby,
kwargs...,
@@ -338,3 +321,17 @@ function _eigvals!(S::Schur{T}; tol = eps(real(T))) where {T}
end
return vals
end
+
+## eigen!
+function eigen!(
+ A::StridedMatrix;
+ sortby::Union{Function,Nothing} = LinearAlgebra.eigsortby,
+ kwargs...,
+)
+
+ if ishermitian(A)
+ return eigen!(Hermitian(A); sortby)
+ end
+
+ throw(ArgumentError("eigen! for general matrices not yet supported. Consider using schur!"))
+end
\ No newline at end of file
diff --git a/src/eigenSelfAdjoint.jl b/src/eigenSelfAdjoint.jl
index 140866c..ecfccfa 100644
--- a/src/eigenSelfAdjoint.jl
+++ b/src/eigenSelfAdjoint.jl
@@ -1,6 +1,4 @@
-using Printf
-using LinearAlgebra
-using LinearAlgebra: givensAlgorithm
+eigeltype(T::Type) = typeof(sqrt(zero(T)))
## EigenQ
struct EigenQ{T} <: AbstractMatrix{T}
@@ -162,18 +160,28 @@ function eig2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matrix
return c, s
end
-function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) where {T<:Real}
+function eigvalsPWK!(
+ S::SymTridiagonal{T};
+ tol = eps(T),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(S, 1)
+) where {T<:Real}
d = S.dv
e = S.ev
n = length(d)
blockstart = 1
blockend = n
iter = 0
+ i = 0
@inbounds begin
for i = 1:n-1
e[i] = abs2(e[i])
end
while true
+ i += 1
+ if i > maxiter
+ throw(ArgumentError("iteration limit $maxiter reached"))
+ end
# 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]))
@@ -225,7 +233,8 @@ end
function eigQL!(
S::SymTridiagonal{T};
vectors::Matrix = zeros(T, 0, size(S, 1)),
- tol = eps(T),
+ tol::AbstractFloat = eps(T),
+ maxiter::Int = 30 * size(S, 1)
) where {T<:Real}
d = S.dv
e = S.ev
@@ -248,43 +257,46 @@ function eigQL!(
blockstart = 1
blockend = n
- @inbounds begin
- 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]))
- blockend = be - 1
- break
- end
- blockend = n
+ i = 0
+ @inbounds while true
+ i += 1
+ if i > maxiter
+ throw(ArgumentError("iteration limit $maxiter reached"))
+ end
+ # 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]))
+ blockend = be - 1
+ break
end
+ blockend = n
+ end
- @debug "submatrix" blockstart blockend
-
- # Deflate?
- if blockstart == blockend
- @debug "Deflate? Yes!"
- blockstart += 1
- elseif blockstart + 1 == blockend
- @debug "Deflate? Yes, but after solving 2x2 block"
- eig2x2!(d, e, blockstart, vectors)
- blockstart += 2
- else
- @debug "Deflate? No!"
- # Calculate shift
- μ = (d[blockstart + 1] - d[blockstart]) / (2 * e[blockstart])
- r = hypot(μ, one(T))
- μ = d[blockstart] - (e[blockstart] / (μ + copysign(r, μ)))
-
- # QL bulk chase
- @debug "Values before bulge chase" e[blockstart] e[blockend - 1] d[blockstart] μ
- singleShiftQL!(S, μ, blockstart, blockend, vectors)
- @debug "Values after bulge chase" e[blockstart] e[blockend - 1] d[blockstart]
- end
+ @debug "submatrix" blockstart blockend
+
+ # Deflate?
+ if blockstart == blockend
+ @debug "Deflate? Yes!"
+ blockstart += 1
+ elseif blockstart + 1 == blockend
+ @debug "Deflate? Yes, but after solving 2x2 block"
+ eig2x2!(d, e, blockstart, vectors)
+ blockstart += 2
+ else
+ @debug "Deflate? No!"
+ # Calculate shift
+ μ = (d[blockstart + 1] - d[blockstart]) / (2 * e[blockstart])
+ r = hypot(μ, one(T))
+ μ = d[blockstart] - (e[blockstart] / (μ + copysign(r, μ)))
+
+ # QL bulk chase
+ @debug "Values before bulge chase" e[blockstart] e[blockend - 1] d[blockstart] μ
+ singleShiftQL!(S, μ, blockstart, blockend, vectors)
+ @debug "Values after bulge chase" e[blockstart] e[blockend - 1] d[blockstart]
+ end
- if blockstart >= n
- break
- end
+ if blockstart >= n
+ break
end
end
p = sortperm(d)
@@ -294,49 +306,53 @@ end
function eigQR!(
S::SymTridiagonal{T};
vectors::Matrix = zeros(T, 0, size(S, 1)),
- tol = eps(T),
+ tol::AbstractFloat = eps(T),
+ maxiter::Int = 30 * size(S, 1)
) where {T<:Real}
d = S.dv
e = S.ev
n = length(d)
blockstart = 1
blockend = n
- @inbounds begin
- 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]))
- blockend = be - 1
- break
- end
- blockend = n
- end
-
- @debug "submatrix" blockstart blockend
-
- # Deflate?
- if blockstart == blockend
- @debug "Deflate? Yes!"
- blockstart += 1
- elseif blockstart + 1 == blockend
- @debug "Deflate? Yes, but after solving 2x2 block!"
- eig2x2!(d, e, blockstart, vectors)
- blockstart += 2
- else
- @debug "Deflate? No!"
- # Calculate shift
- μ = (d[blockend - 1] - d[blockend]) / (2 * e[blockend - 1])
- r = hypot(μ, one(T))
- μ = d[blockend] - (e[blockend - 1] / (μ + copysign(r, μ)))
-
- # QR bulk chase
- @debug "Values before QR bulge chase" e[blockstart] e[blockend - 1] d[blockend] μ
- singleShiftQR!(S, μ, blockstart, blockend, vectors)
- @debug "Values after QR bulge chase" e[blockstart] e[blockend - 1] d[blockend]
- end
- if blockstart >= n
+ i = 0
+ @inbounds while true
+ i += 1
+ if i > maxiter
+ throw(ArgumentError("iteration limit $maxiter reached"))
+ end
+ # 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]))
+ blockend = be - 1
break
end
+ blockend = n
+ end
+
+ @debug "submatrix" blockstart blockend
+
+ # Deflate?
+ if blockstart == blockend
+ @debug "Deflate? Yes!"
+ blockstart += 1
+ elseif blockstart + 1 == blockend
+ @debug "Deflate? Yes, but after solving 2x2 block!"
+ eig2x2!(d, e, blockstart, vectors)
+ blockstart += 2
+ else
+ @debug "Deflate? No!"
+ # Calculate shift
+ μ = (d[blockend - 1] - d[blockend]) / (2 * e[blockend - 1])
+ r = hypot(μ, one(T))
+ μ = d[blockend] - (e[blockend - 1] / (μ + copysign(r, μ)))
+
+ # QR bulk chase
+ @debug "Values before QR bulge chase" e[blockstart] e[blockend - 1] d[blockend] μ
+ singleShiftQR!(S, μ, blockstart, blockend, vectors)
+ @debug "Values after QR bulge chase" e[blockstart] e[blockend - 1] d[blockend]
+ end
+ if blockstart >= n
+ break
end
end
p = sortperm(d)
@@ -576,114 +592,188 @@ function symtriUpper!(
end
-_eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
- eigvalsPWK!(A.diagonals; tol, sortby)
-
-_eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvalsPWK!(A; tol, sortby)
-
-_eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvals!(symtri!(A); tol, sortby)
-
-_eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvals!(symtri!(A); tol, sortby)
-
-LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
- _eigvals!(A; tol, sortby)
-
-LinearAlgebra.eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
- _eigvals!(A; tol, sortby)
-
-LinearAlgebra.eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby)
-
-LinearAlgebra.eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby)
-
-_eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
- LinearAlgebra.Eigen(LinearAlgebra.sorteig!(eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol)..., sortby)...)
-
-_eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = LinearAlgebra.Eigen(
- eigQL!(A, vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), tol = tol)...,
+eigvals!(
+ A::SymmetricTridiagonalFactorization;
+ tol = eps(real(eltype(A))),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigvalsPWK!(
+ A.diagonals;
+ tol,
+ sortby,
+ maxiter
)
-_eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(symtri!(A), tol = tol)
-
-_eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(symtri!(A), tol = tol)
+eigvals!(
+ A::SymTridiagonal;
+ tol::AbstractFloat = eps(real(eltype(A))),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigvalsPWK!(A; tol, sortby, maxiter)
+
+eigvals!(
+ A::Hermitian;
+ tol::AbstractFloat = eps(real(eltype(A))),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigvals!(symtri!(A); tol, sortby, maxiter)
+
+eigvals!(
+ A::Symmetric{<:Real};
+ tol::AbstractFloat = eps(eltype(A)),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigvals!(
+ symtri!(A);
+ tol,
+ sortby,
+ maxiter
+)
-LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
- _eigen!(A; tol, sortby)
+eigvals(
+ A::AbstractMatrix;
+ tol::AbstractFloat = eps(real(eigeltype(eltype(A)))),
+ sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigvals!(
+ _eigencopy_oftype(
+ A,
+ eigeltype(eltype(A)),
+ );
+ tol,
+ sortby,
+ maxiter
+)
-LinearAlgebra.eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
+eigen!(
+ A::SymmetricTridiagonalFactorization;
+ tol::AbstractFloat = eps(real(eltype(A))),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = LinearAlgebra.Eigen(
+ LinearAlgebra.sorteig!(
+ eigQL!(
+ A.diagonals;
+ vectors = Array(A.Q),
+ tol,
+ maxiter
+ )...,
+ sortby
+ )...
+)
-LinearAlgebra.eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
+function eigen!(
+ A::SymTridiagonal;
+ tol = eps(real(eltype(A))),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+)
+ return LinearAlgebra.Eigen(
+ LinearAlgebra.sorteig!(
+ eigQL!(
+ A;
+ vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)),
+ tol,
+ maxiter
+ )...,
+ sortby
+ )...,
+ )
+end
-LinearAlgebra.eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
+eigen!(
+ A::Hermitian;
+ tol = eps(real(eltype(A))),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigen!(symtri!(A); tol, sortby, maxiter)
+
+eigen!(
+ A::Symmetric{<:Real};
+ tol = eps(eltype(A)),
+ sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigen!(symtri!(A); tol, sortby, maxiter)
+
+eigen(
+ A::AbstractMatrix;
+ tol::AbstractFloat = eps(real(eigeltype(eltype(A)))),
+ sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+ maxiter::Int = 30 * size(A, 1)
+) = eigen!(
+ _eigencopy_oftype(
+ A,
+ eigeltype(eltype(A)),
+ );
+ tol,
+ sortby,
+ maxiter
+)
+# We need this method because eigencopy_oftype converts SymTridiagonal to Matrix.
+# That is probably a bad idea so maybe this should be filed upstream
+eigen(
+ A::SymTridiagonal;
+ tol::AbstractFloat = eps(real(eigeltype(eltype(A)))),
+ sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby
+) = eigen!(
+ LinearAlgebra.copy_oftype(
+ A,
+ eigeltype(eltype(A)),
+ );
+ tol,
+ sortby
+)
function eigen2!(
A::SymmetricTridiagonalFactorization;
- tol = eps(real(float(one(eltype(A))))),
+ tol::AbstractFloat = eps(real(float(one(eltype(A))))),
+ sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby
)
V = zeros(eltype(A), 2, size(A, 1))
V[1] = 1
V[end] = 1
- eigQL!(A.diagonals, vectors = rmul!(V, A.Q), tol = tol)
+ return LinearAlgebra.sorteig!(
+ eigQL!(A.diagonals; vectors = rmul!(V, A.Q), tol)...,
+ sortby
+ )
end
-function eigen2!(A::SymTridiagonal; tol = eps(real(float(one(eltype(A))))))
+function eigen2!(
+ A::SymTridiagonal;
+ tol = eps(real(float(one(eltype(A))))),
+ sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby
+)
V = zeros(eltype(A), 2, size(A, 1))
V[1] = 1
V[end] = 1
- eigQL!(A, vectors = V, tol = tol)
+ return LinearAlgebra.sorteig!(
+ eigQL!(A; vectors = V, tol)...,
+ sortby
+ )
end
-eigen2!(A::Hermitian; tol = eps(float(real(one(eltype(A)))))) =
- eigen2!(symtri!(A), tol = tol)
-
-eigen2!(A::Symmetric{<:Real}; tol = eps(float(one(eltype(A))))) =
- eigen2!(symtri!(A), tol = tol)
-
+eigen2!(A::Hermitian; tol = eps(float(real(one(eltype(A))))), sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+ eigen2!(symtri!(A); tol, sortby)
-eigen2(A::SymTridiagonal; tol = eps(float(real(one(eltype(A)))))) =
- eigen2!(copy(A), tol = tol)
+eigen2!(A::Symmetric{<:Real}; tol = eps(float(one(eltype(A)))), sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+ eigen2!(symtri!(A); tol, sortby)
-eigen2(A::Hermitian, tol = eps(float(real(one(eltype(A)))))) = eigen2!(copy(A), tol = tol)
-eigen2(A::Symmetric{<:Real}, tol = eps(float(one(eltype(A))))) = eigen2!(copy(A), tol = tol)
+eigen2(A::SymTridiagonal; tol = eps(float(real(one(eltype(A))))), sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+ eigen2!(copy(A); tol, sortby)
-# First method of each type here is identical to the method defined in
-# LinearAlgebra but is needed for disambiguation
-const _eigencopy_oftype = if VERSION >= v"1.9"
- LinearAlgebra.eigencopy_oftype
-else
- LinearAlgebra.copy_oftype
-end
+eigen2(A::Hermitian, tol = eps(float(real(one(eltype(A))))), sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+ eigen2!(copy(A); tol, sortby)
-if VERSION < v"1.7"
- function LinearAlgebra.eigvals(A::Hermitian{<:Real})
- T = typeof(sqrt(zero(eltype(A))))
- return eigvals!(_eigencopy_oftype(A, T))
- end
- function LinearAlgebra.eigvals(A::Hermitian{<:Complex})
- T = typeof(sqrt(zero(eltype(A))))
- return eigvals!(_eigencopy_oftype(A, T))
- end
- function LinearAlgebra.eigvals(A::Hermitian)
- T = typeof(sqrt(zero(eltype(A))))
- return eigvals!(_eigencopy_oftype(A, T))
- end
- function LinearAlgebra.eigen(A::Hermitian{<:Real})
- T = typeof(sqrt(zero(eltype(A))))
- return eigen!(_eigencopy_oftype(A, T))
- end
- function LinearAlgebra.eigen(A::Hermitian{<:Complex})
- T = typeof(sqrt(zero(eltype(A))))
- return eigen!(_eigencopy_oftype(A, T))
- end
- function LinearAlgebra.eigen(A::Hermitian)
- T = typeof(sqrt(zero(eltype(A))))
- return eigen!(_eigencopy_oftype(A, T))
- end
-end
+eigen2(A::Symmetric{<:Real}; tol = eps(float(one(eltype(A)))), sortby::Union{Function,Nothing} = sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+ eigen2!(copy(A); tol, sortby)
# Aux (should go somewhere else at some point)
-function LinearAlgebra.givensAlgorithm(f::Real, g::Real)
+function givensAlgorithm(f::Real, g::Real)
h = hypot(f, g)
return f / h, g / h, h
end
+
+_eigencopy_oftype(A::SymTridiagonal, ::Type{T}) where T = LinearAlgebra.copy_oftype(A, T)
+_eigencopy_oftype(A::AbstractMatrix, ::Type{T}) where T = LinearAlgebra.eigencopy_oftype(A, T)
diff --git a/src/householder.jl b/src/householder.jl
index b34e1f1..bf68760 100644
--- a/src/householder.jl
+++ b/src/householder.jl
@@ -1,6 +1,3 @@
-import Base: *, eltype, size
-import LinearAlgebra: adjoint, mul!, rmul!
-
struct Householder{T,S<:StridedVector}
v::S
τ::T
@@ -10,14 +7,14 @@ struct HouseholderBlock{T,S<:StridedMatrix,U<:StridedMatrix}
T::UpperTriangular{T,U}
end
-size(H::Householder) = (length(H.v) + 1, length(H.v) + 1)
-size(H::Householder, i::Integer) = i <= 2 ? length(H.v) + 1 : 1
+Base.size(H::Householder) = (length(H.v) + 1, length(H.v) + 1)
+Base.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
+Base.eltype(H::Householder{T}) where {T} = T
+Base.eltype(H::HouseholderBlock{T}) where {T} = T
-adjoint(H::Householder{T}) where {T} = Adjoint{T,typeof(H)}(H)
-adjoint(H::HouseholderBlock{T}) where {T} = Adjoint{T,typeof(H)}(H)
+LinearAlgebra.adjoint(H::Householder{T}) where {T} = Adjoint{T,typeof(H)}(H)
+LinearAlgebra.adjoint(H::HouseholderBlock{T}) where {T} = Adjoint{T,typeof(H)}(H)
function lmul!(H::Householder, A::StridedMatrix)
m, n = size(A)
@@ -113,7 +110,7 @@ function lmul!(H::HouseholderBlock{T}, A::StridedMatrix{T}, M::StridedMatrix{T})
return A
end
-(*)(H::HouseholderBlock{T}, A::StridedMatrix{T}) where {T} =
+Base.:(*)(H::HouseholderBlock{T}, A::StridedMatrix{T}) where {T} =
lmul!(H, copy(A), similar(A, (min(size(H.V)...), size(A, 2))))
function lmul!(
@@ -155,7 +152,7 @@ function lmul!(
return A
end
-(*)(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}) where {T} =
+Base.:(*)(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}) where {T} =
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)
diff --git a/src/juliaBLAS.jl b/src/juliaBLAS.jl
index cb40e8b..6c79afe 100644
--- a/src/juliaBLAS.jl
+++ b/src/juliaBLAS.jl
@@ -1,10 +1,3 @@
-using LinearAlgebra
-using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, HermOrSym
-
-import LinearAlgebra: lmul!, mul!
-
-export rankUpdate!
-
# Rank one update
## General
diff --git a/src/lapack.jl b/src/lapack.jl
index 39a10f7..4cfe0c3 100644
--- a/src/lapack.jl
+++ b/src/lapack.jl
@@ -1,12 +1,10 @@
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
+using LinearAlgebra: LinearAlgebra
+using LinearAlgebra.BLAS: BlasInt,
+ chkstride1
+using LinearAlgebra.LAPACK: @blasfunc
+liblapack_name = LinearAlgebra.libblastrampoline_jll.libblastrampoline
## Standard QR/QL
function steqr!(
diff --git a/src/qr.jl b/src/qr.jl
index 323a5d0..b8567e9 100644
--- a/src/qr.jl
+++ b/src/qr.jl
@@ -1,8 +1,3 @@
-using LinearAlgebra
-
-import Base: getindex, size
-import LinearAlgebra: reflectorApply!
-
struct QR2{T,S<:AbstractMatrix{T},V<:AbstractVector{T}} <: Factorization{T}
factors::S
τ::V
@@ -13,7 +8,7 @@ struct QR2{T,S<:AbstractMatrix{T},V<:AbstractVector{T}} <: Factorization{T}
end
QR2(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = QR2{T,typeof(factors)}(factors, τ)
-size(F::QR2, i::Integer...) = size(F.factors, i...)
+Base.size(F::QR2, i::Integer...) = size(F.factors, i...)
# Similar to the definition in base but applies the reflector from the right
@inline function reflectorApply!(A::StridedMatrix, x::AbstractVector, τ::Number) # apply conjugate transpose reflector from right.
@@ -41,6 +36,8 @@ size(F::QR2, i::Integer...) = size(F.factors, i...)
return A
end
+reflectorApply!(x::AbstractVector, τ::Number, A::StridedVecOrMat) = LinearAlgebra.reflectorApply!(x, τ, A)
+
# FixMe! Consider how to represent Q
# immutable Q{T,S<:QR2} <: AbstractMatrix{T}
@@ -52,7 +49,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 Base.getindex(A::QR2{T}, ::Type{Tuple{:R}}) where {T}
m, n = size(A)
if m >= n
UpperTriangular(view(A.factors, 1:n, 1:n))
@@ -61,7 +58,7 @@ function getindex(A::QR2{T}, ::Type{Tuple{:R}}) where {T}
end
end
-function getindex(A::QR2{T}, ::Type{Tuple{:QBlocked}}) where {T}
+function Base.getindex(A::QR2{T}, ::Type{Tuple{:QBlocked}}) where {T}
m, n = size(A)
mmn = min(m, n)
F = A.factors
diff --git a/src/svd.jl b/src/svd.jl
index 656ebb0..8254c7f 100644
--- a/src/svd.jl
+++ b/src/svd.jl
@@ -1,11 +1,11 @@
-using LinearAlgebra
-
-import LinearAlgebra: mul!, rmul!
-
AdjointQtype = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint
+svdeltype(T::Type) = typeof(sqrt(zero(T)))
+
lmul!(G::LinearAlgebra.Givens, ::Nothing) = nothing
+lmul!(A::Any, B::Any) = LinearAlgebra.lmul!(A, B)
rmul!(::Nothing, G::LinearAlgebra.Givens) = nothing
+rmul!(A::Any, B::Any) = LinearAlgebra.rmul!(A, B)
function svdvals2x2(d1, d2, e)
d1sq = d1 * d1
@@ -266,7 +266,7 @@ function _sort_and_adjust!(U, s, Vᴴ)
n = length(s)
# Adjust sign of singular values if necessary
- if any(!isposdef, s)
+ if any(!LinearAlgebra.isposdef, s)
rmul!(U, Diagonal(_sign.(s)))
map!(abs, s, s)
end
@@ -322,7 +322,7 @@ function bidiagonalize!(A::AbstractMatrix)
conj!(x)
τi = LinearAlgebra.reflector!(x)
push!(τr, τi)
- LinearAlgebra.reflectorApply!(view(A, (i+1):m, (i+1):n), x, τi)
+ reflectorApply!(view(A, (i+1):m, (i+1):n), x, τi)
end
end
@@ -341,7 +341,7 @@ function bidiagonalize!(A::AbstractMatrix)
conj!(x)
τi = LinearAlgebra.reflector!(x)
push!(τr, τi)
- LinearAlgebra.reflectorApply!(view(A, (i+1):m, i:n), x, τi)
+ reflectorApply!(view(A, (i+1):m, i:n), x, τi)
if i < m
x = view(A, i+1:m, i)
τi = LinearAlgebra.reflector!(x)
@@ -410,8 +410,9 @@ function Base.getproperty(F::BidiagonalFactorization, s::Symbol)
end
end
+# PIRACY ALLERT!
# For now, we define a generic lmul! and rmul! for HessenbergQ here
-function lmul!(Q::LinearAlgebra.HessenbergQ, B::AbstractVecOrMat)
+function LinearAlgebra.lmul!(Q::LinearAlgebra.HessenbergQ, B::AbstractVecOrMat)
m, n = size(B, 1), size(B, 2)
@@ -436,7 +437,7 @@ function lmul!(Q::LinearAlgebra.HessenbergQ, B::AbstractVecOrMat)
return B
end
-function lmul!(adjQ::AdjointQtype{<:Any,<:LinearAlgebra.HessenbergQ}, B::AbstractVecOrMat)
+function LinearAlgebra.lmul!(adjQ::AdjointQtype{<:Any,<:LinearAlgebra.HessenbergQ}, B::AbstractVecOrMat)
Q = parent(adjQ)
m, n = size(B, 1), size(B, 2)
@@ -462,7 +463,7 @@ function lmul!(adjQ::AdjointQtype{<:Any,<:LinearAlgebra.HessenbergQ}, B::Abstrac
return B
end
-function rmul!(A::AbstractMatrix, Q::LinearAlgebra.HessenbergQ)
+function LinearAlgebra.rmul!(A::AbstractMatrix, Q::LinearAlgebra.HessenbergQ)
m, n = size(A)
@@ -487,7 +488,7 @@ function rmul!(A::AbstractMatrix, Q::LinearAlgebra.HessenbergQ)
return A
end
-function rmul!(A::AbstractMatrix, adjQ::AdjointQtype{<:Any,<:LinearAlgebra.HessenbergQ})
+function LinearAlgebra.rmul!(A::AbstractMatrix, adjQ::AdjointQtype{<:Any,<:LinearAlgebra.HessenbergQ})
m, n = size(A)
Q = parent(adjQ)
@@ -541,7 +542,7 @@ end
# Overload LinearAlgebra methods
-LinearAlgebra.svdvals!(B::Bidiagonal{T}; tol = eps(T)) where {T<:Real} =
+svdvals!(B::Bidiagonal{T}; tol = eps(T)) where {T<:Real} =
_svdvals!(B, tol = tol)
"""
@@ -566,17 +567,21 @@ julia> svdvals(A) ≈ [3, 2, 1]
true
```
"""
-function LinearAlgebra.svdvals!(A::StridedMatrix; tol = eps(real(eltype(A))))
+function svdvals!(A::StridedMatrix; tol = eps(real(eltype(A))))
B = bidiagonalize!(A).bidiagonal
# It doesn't matter that we transpose the bidiagonal matrix when we are only interested in the values
return svdvals!(Bidiagonal(B.dv, B.ev, :U), tol = tol)
end
+svdvals(A::AbstractMatrix; tol = eps(real(svdeltype(eltype(A))))) =
+ svdvals!(LinearAlgebra.copy_oftype(A, svdeltype(eltype(A))))
+
+
# FixMe! The full keyword is redundant for Bidiagonal and should be removed from Base
-LinearAlgebra.svd!(
+svd!(
B::Bidiagonal{T};
- tol = eps(T),
- full = false,
+ tol::AbstractFloat = eps(T),
+ full::Bool = false,
# To avoid breaking on <Julia 1.3, the `alg` keyword doesn't do anything. Once we drop support for Julia 1.2
# and below, we can make the keyword argument work correctly
alg = nothing,
@@ -613,13 +618,10 @@ Vt factor:
-0.817416 0.576048
```
"""
-function LinearAlgebra.svd!(
+function svd!(
A::StridedMatrix{T};
- tol = eps(real(eltype(A))),
- full = false,
- # To avoid breaking on <Julia 1.3, the `alg` keyword doesn't do anything. Once we drop support for Julia 1.2
- # and below, we can make the keyword argument work correctly
- alg = nothing,
+ tol::AbstractFloat = eps(real(eltype(A))),
+ full::Bool = false,
) where {T}
m, n = size(A)
@@ -664,3 +666,15 @@ function LinearAlgebra.svd!(
return SVD(U, s, Vᴴ)
end
+
+svd(
+ A::AbstractMatrix{T};
+ tol::AbstractFloat = eps(real(svdeltype(T))),
+ full::Bool = false,
+) where {T} = svd!(LinearAlgebra.copy_oftype(A, svdeltype(T)); tol, full)
+
+## Definitions that rely on the SVD
+function cond(A::AbstractMatrix)
+ vals = svdvals(A)
+ return first(vals) / last(vals)
+end
diff --git a/src/tridiag.jl b/src/tridiag.jl
deleted file mode 100644
index a863383..0000000
--- a/src/tridiag.jl
+++ /dev/null
@@ -1,36 +0,0 @@
-################################################
-## Specialized routines for tridiagonal matrices
-################################################
-
-export numnegevals
-
-"""
-Computes the number of negative eigenvalues of T - σI, a.k.a. spectrum slicing
-
-Inputs:
- T: A SymTridiagonal{<:Real} matrix
- σ: The shift parameter
-
-Outputs:
- ν: The number of negative eigenvalues
-
-Reference:
- B. N. Parlett, "The symmetric eigenvalue problem", Section 3.3.1, p. 52.
-"""
-function numnegevals(T::SymTridiagonal{S}, σ::S = zero(S)) where {S}
- α = T.dv
- β = T.ev
- n = length(α)
- ϵ = eps(S)
- δ = α[1] - σ
- ν = δ < 0 ? 1 : 0
- for k = 1:n-1
- if δ == 0
- info("zero in iteration $k")
- δ = ϵ * (β[k] + ϵ) #Parlett prefers adjusting σ and starting again
- end
- δ = (α[k+1] - σ) - β[k] * (β[k] / δ)
- ν += (δ < 0)
- end
- ν
-end
diff --git a/test/cholesky.jl b/test/cholesky.jl
index 0f92f45..81444c7 100644
--- a/test/cholesky.jl
+++ b/test/cholesky.jl
@@ -1,5 +1,5 @@
-using Test, LinearAlgebra, Random
-using GenericLinearAlgebra: cholUnblocked!, cholBlocked!, cholRecursive!
+using Test, Random, GenericLinearAlgebra
+using LinearAlgebra: LinearAlgebra
@testset "Cholesky" begin
@@ -15,16 +15,16 @@ using GenericLinearAlgebra: cholUnblocked!, cholBlocked!, cholRecursive!
A = T.(rand(n, n))
end
AcA = A'A
- @test LowerTriangular(cholUnblocked!(copy(AcA), Val{:L})) ≈
- cholesky(Hermitian(AcA)).L
- @test LowerTriangular(cholBlocked!(copy(AcA), Val{:L}, 5)) ≈
- cholesky(Hermitian(AcA)).L
- @test LowerTriangular(cholBlocked!(copy(AcA), Val{:L}, 10)) ≈
- cholesky(Hermitian(AcA)).L
- @test LowerTriangular(cholRecursive!(copy(AcA), Val{:L}, 1)) ≈
- cholesky(Hermitian(AcA)).L
- @test LowerTriangular(cholRecursive!(copy(AcA), Val{:L}, 4)) ≈
- cholesky(Hermitian(AcA)).L
+ @test LinearAlgebra.LowerTriangular(GenericLinearAlgebra.cholUnblocked!(copy(AcA), Val{:L})) ≈
+ LinearAlgebra.cholesky(LinearAlgebra.Hermitian(AcA)).L
+ @test LinearAlgebra.LowerTriangular(GenericLinearAlgebra.cholBlocked!(copy(AcA), Val{:L}, 5)) ≈
+ LinearAlgebra.cholesky(LinearAlgebra.Hermitian(AcA)).L
+ @test LinearAlgebra.LowerTriangular(GenericLinearAlgebra.cholBlocked!(copy(AcA), Val{:L}, 10)) ≈
+ LinearAlgebra.cholesky(LinearAlgebra.Hermitian(AcA)).L
+ @test LinearAlgebra.LowerTriangular(GenericLinearAlgebra.cholRecursive!(copy(AcA), Val{:L}, 1)) ≈
+ LinearAlgebra.cholesky(LinearAlgebra.Hermitian(AcA)).L
+ @test LinearAlgebra.LowerTriangular(GenericLinearAlgebra.cholRecursive!(copy(AcA), Val{:L}, 4)) ≈
+ LinearAlgebra.cholesky(LinearAlgebra.Hermitian(AcA)).L
end
end
end
diff --git a/test/eigengeneral.jl b/test/eigengeneral.jl
index cfc7acd..c13e94a 100644
--- a/test/eigengeneral.jl
+++ b/test/eigengeneral.jl
@@ -1,4 +1,5 @@
-using Test, GenericLinearAlgebra, LinearAlgebra
+using Test, GenericLinearAlgebra
+using LinearAlgebra: LinearAlgebra
@testset "The General eigenvalue problem" begin
@@ -9,8 +10,8 @@ using Test, GenericLinearAlgebra, LinearAlgebra
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
+ vLAPACK = LinearAlgebra.eigvals(A)
+ vBig = GenericLinearAlgebra.eigvals(big.(A))
@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)
@@ -19,10 +20,10 @@ using Test, GenericLinearAlgebra, LinearAlgebra
@testset "Rayleigh shifts" begin
@test sort(
GenericLinearAlgebra._eigvals!(
- GenericLinearAlgebra._schur!(copy(A), shiftmethod = :Rayleigh),
+ GenericLinearAlgebra.schur!(copy(A), shiftmethod = :Rayleigh),
),
by = t -> (real(t), imag(t)),
- ) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
+ ) ≈ sort(LinearAlgebra.eigvals(A), by = t -> (real(t), imag(t)))
end
end
end
@@ -30,7 +31,7 @@ using Test, GenericLinearAlgebra, LinearAlgebra
@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)))
+ sort(abs.(GenericLinearAlgebra._eigvals!(copy(A)))) ≈ sort(abs.(LinearAlgebra.eigvals(A)))
end
end
@@ -57,28 +58,28 @@ using Test, GenericLinearAlgebra, LinearAlgebra
A = my_matrix(4, 1e-3)
@test sort(
- GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))),
+ GenericLinearAlgebra._eigvals!(GenericLinearAlgebra.schur!(copy(A))),
by = t -> (real(t), imag(t)),
- ) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
+ ) ≈ sort(LinearAlgebra.eigvals(A), by = t -> (real(t), imag(t)))
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))),
+ GenericLinearAlgebra._eigvals!(GenericLinearAlgebra.schur!(copy(A))),
by = t -> (real(t), imag(t)),
- ) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
+ ) ≈ sort(LinearAlgebra.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))),
+ GenericLinearAlgebra._eigvals!(GenericLinearAlgebra.schur!(copy(B))),
by = t -> (real(t), imag(t)),
- ) ≈ sort(eigvals(B), by = t -> (real(t), imag(t)))
+ ) ≈ sort(LinearAlgebra.eigvals(B), by = t -> (real(t), imag(t)))
end
@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)))
+ @test sum(LinearAlgebra.eigvals(LinearAlgebra.schur(A).T)) ≈ sum(LinearAlgebra.eigvals(GenericLinearAlgebra.schur!(copy(A)).T))
+ @test sum(LinearAlgebra.eigvals(LinearAlgebra.schur(A).Schur)) ≈ sum(LinearAlgebra.eigvals(GenericLinearAlgebra.schur!(copy(A)).Schur))
end
@testset "Issue 63" begin
@@ -107,7 +108,7 @@ using Test, GenericLinearAlgebra, LinearAlgebra
(7 + [2, 2]' * z2),
]) / 3
- @test eigvals(big.(A)) ≈ truevals
+ @test GenericLinearAlgebra.eigvals(big.(A)) ≈ truevals
end
Demmel(η) = [
@@ -120,7 +121,7 @@ using Test, GenericLinearAlgebra, LinearAlgebra
@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))
+ vals = GenericLinearAlgebra._eigvals!(GenericLinearAlgebra.schur!(A, maxiter = 35))
@test abs.(vals) ≈ ones(4)
end
@@ -148,7 +149,7 @@ using Test, GenericLinearAlgebra, LinearAlgebra
-0.313987810419091240,
)
- @test Float64.(abs.(eigvals(big.(H)))) ≈ ones(4)
+ @test Float64.(abs.(GenericLinearAlgebra.eigvals(big.(H)))) ≈ ones(4)
end
@testset "Issue 67" for (A, λs) in (
@@ -224,12 +225,12 @@ using Test, GenericLinearAlgebra, LinearAlgebra
)
@testset "shouldn't need excessive iterations (30*n) in the Float64 case" begin
- GenericLinearAlgebra._schur!(float(A))
+ GenericLinearAlgebra.schur!(float(A))
end
# For BigFloats, many iterations are required for convergence
# Improving this is probably a hard problem
- vals = eigvals(big.(A), maxiter = 1500)
+ vals = GenericLinearAlgebra.eigvals(big.(A), maxiter = 1500)
# 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
@@ -245,10 +246,10 @@ using Test, GenericLinearAlgebra, LinearAlgebra
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
+ @test LinearAlgebra.tril(A, -2) ≈ zeros(n, n) atol = 1e-14
- @test eigvals(HF.H) ≈ eigvals(A)
- @test eigvals(HF.H) ≈ eigvals!(copy(HF))
+ @test GenericLinearAlgebra.eigvals(HF.H) ≈ LinearAlgebra.eigvals(A)
+ @test GenericLinearAlgebra.eigvals(HF.H) ≈ GenericLinearAlgebra.eigvals!(copy(HF))
@test HF.H \ ones(n) ≈ Matrix(HF.H) \ ones(n)
end
end
diff --git a/test/eigenselfadjoint.jl b/test/eigenselfadjoint.jl
index 3fe8eb7..6e4e047 100644
--- a/test/eigenselfadjoint.jl
+++ b/test/eigenselfadjoint.jl
@@ -1,4 +1,6 @@
-using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
+using Test, GenericLinearAlgebra, Quaternions
+using LinearAlgebra: LinearAlgebra,
+ Diagonal, Hermitian, I, SymTridiagonal, Symmetric, Tridiagonal
@testset "The selfadjoint eigen problem" begin
n = 50
@@ -7,11 +9,11 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
# 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 = GenericLinearAlgebra.eigen(T)
@test issorted(vals)
@testset "default" begin
@test (vecs' * T) * vecs ≈ Diagonal(vals)
- @test eigvals(T) ≈ vals
+ @test GenericLinearAlgebra.eigvals(T) ≈ vals
@test vecs'vecs ≈ Matrix(I, n, n)
end
@@ -28,17 +30,17 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
GenericLinearAlgebra.eigQR!(copy(T), vectors = Matrix{eltype(T)}(I, n, n))
@test issorted(vals)
@test (vecs' * T) * vecs ≈ Diagonal(vals)
- @test eigvals(T) ≈ vals
+ @test GenericLinearAlgebra.eigvals(T) ≈ vals
@test vecs'vecs ≈ Matrix(I, n, n)
end
end
@testset "(full) Symmetric" for uplo in (:L, :U)
A = Symmetric(big.(randn(n, n)), uplo)
- vals, vecs = @inferred(eigen(A))
+ vals, vecs = GenericLinearAlgebra.eigen(A)
@testset "default" begin
- @test vecs' * A * vecs ≈ diagm(0 => vals)
- @test @inferred(eigvals(A)) ≈ vals
+ @test vecs' * A * vecs ≈ LinearAlgebra.diagm(0 => vals)
+ @test @inferred(GenericLinearAlgebra.eigvals(A)) ≈ vals
@test vecs'vecs ≈ Matrix(I, n, n)
@test issorted(vals)
end
@@ -53,19 +55,19 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
end
@testset "(full) Quaternion Hermitian using :$uplo" for uplo in (:L, :U)
- V = qr([Quaternion(randn(4)...) for i = 1:n, j = 1:n]).Q
+ V = LinearAlgebra.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)
+ vals, vecs = GenericLinearAlgebra.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 ≈ LinearAlgebra.diagm(0 => vals)
else
- @test vecs' * A * vecs ≈ diagm(0 => vals)
+ @test vecs' * A * vecs ≈ LinearAlgebra.diagm(0 => vals)
end
- @test eigvals(A) ≈ vals
+ @test GenericLinearAlgebra.eigvals(A) ≈ vals
@test vals ≈ λ rtol = 1e-13 * n
@test vecs'vecs ≈ Matrix(I, n, n)
end
@@ -82,21 +84,21 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
@testset "big Hermitian{<:Complex}" begin
# This one used to cause an ambiguity error. See #35
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.(GenericLinearAlgebra.eigen(Hermitian(big.(A))).values) ≈ LinearAlgebra.eigen(Hermitian(copy(A))).values
+ @test Float64.(GenericLinearAlgebra.eigvals(Hermitian(big.(A)))) ≈ LinearAlgebra.eigvals(Hermitian(copy(A)))
end
@testset "generic Givens" begin
x, y = randn(2)
- c, s, r = invoke(LinearAlgebra.givensAlgorithm, Tuple{Real,Real}, x, y)
+ c, s, r = GenericLinearAlgebra.givensAlgorithm(x, y)
@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 ==
+ @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)
end
@@ -106,16 +108,16 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
M1[1, 1] = 0.01
M2 = Hermitian(zeros(3, 3))
M2[3, 3] = 0.01
- @test eigvals(M1) == GenericLinearAlgebra._eigvals!(copy(M1))
- @test eigvals(M2) == GenericLinearAlgebra._eigvals!(copy(M2))
- @test eigen(M1).values == GenericLinearAlgebra._eigen!(copy(M1)).values
- @test eigen(M2).values == GenericLinearAlgebra._eigen!(copy(M2)).values
+ @test LinearAlgebra.eigvals(M1) == GenericLinearAlgebra.eigvals!(copy(M1))
+ @test LinearAlgebra.eigvals(M2) == GenericLinearAlgebra.eigvals!(copy(M2))
+ @test LinearAlgebra.eigen(M1).values == GenericLinearAlgebra.eigen!(copy(M1)).values
+ @test LinearAlgebra.eigen(M2).values == GenericLinearAlgebra.eigen!(copy(M2)).values
end
@testset "Sorting of `ishermitian(T)==true` matrices on pre-1.2" begin
T = big.(randn(5, 5))
T = T + T'
- @test issorted(eigvals(T))
+ @test issorted(GenericLinearAlgebra.eigvals(T))
end
@testset "issue 123" begin
@@ -125,7 +127,7 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
big"-1.0795157507124452910839896877334667387210301781514938067860918240771876343947" big"0.1882735697944729855991976669864503854920622386133987141371224931350749728226066" big"0.9688026817149176598146701814747478080649943014810992426739997593840858865727305" big"-1.672789745967021000172452940954243617442140494364475046869527486458478435262502";
big"0.9172086645212876240254394768180975107502376572771647296150618931226550446699544" big"-0.1599663084136352437739757607131301560774255778371317602542426234968564801904052" big"-1.672789745967021000172452940954243617442140494364475046869527486458478435262502" big"0.4212828742060771422472975116067336073573584644697624467523583310058490760719874"
])
- F = eigen(M)
+ F = GenericLinearAlgebra.eigen(M)
@test M * F.vectors ≈ F.vectors * Diagonal(F.values)
@testset "2x2 cases, d: $d, e: $e" for (d, e) in (
@@ -147,28 +149,27 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
@testset "#133" begin
A = SymTridiagonal{BigFloat}(randn(5), randn(4))
- T = Tridiagonal(A)
- if VERSION ≥ v"1.10"
- # The following depends on more recent behaviour of eigvals/eigen in stdlib
- @test eigvals(A) == eigvals(T) == eigvals(A; sortby=LinearAlgebra.eigsortby) == eigvals(T; sortby=LinearAlgebra.eigsortby) == eigvals!(deepcopy(A); sortby=LinearAlgebra.eigsortby)
- @test eigen(A).values == eigen(T).values == eigen(A; sortby=LinearAlgebra.eigsortby).values == eigen(T; sortby=LinearAlgebra.eigsortby).values
- # compare abs to avoid sign issues
- @test abs.(eigen(A).vectors) == abs.(eigen(T).vectors) == abs.(eigen(A; sortby=LinearAlgebra.eigsortby).vectors) == abs.(eigen(T; sortby=LinearAlgebra.eigsortby).vectors)
- end
+ T = Matrix(Tridiagonal(A))
+ # The following depends on more recent behaviour of eigvals/eigen in stdlib
+ @test GenericLinearAlgebra.eigvals(A) ==
+ GenericLinearAlgebra.eigvals(T) ==
+ GenericLinearAlgebra.eigvals(A; sortby=LinearAlgebra.eigsortby) ==
+ GenericLinearAlgebra.eigvals(T; sortby=LinearAlgebra.eigsortby) ==
+ GenericLinearAlgebra.eigvals!(deepcopy(A); sortby=LinearAlgebra.eigsortby)
+ @test GenericLinearAlgebra.eigen(A).values ==
+ GenericLinearAlgebra.eigen(T).values ==
+ GenericLinearAlgebra.eigen(A; sortby=LinearAlgebra.eigsortby).values ==
+ GenericLinearAlgebra.eigen(T; sortby=LinearAlgebra.eigsortby).values
+ # compare abs to avoid sign issues
+ @test abs.(GenericLinearAlgebra.eigen(A).vectors) ==
+ abs.(GenericLinearAlgebra.eigen(T).vectors) ==
+ abs.(GenericLinearAlgebra.eigen(A; sortby=LinearAlgebra.eigsortby).vectors) ==
+ abs.(GenericLinearAlgebra.eigen(T; sortby=LinearAlgebra.eigsortby).vectors)
end
- if VERSION >= v"1.9"
- @testset "#139" begin
- A = Hermitian(Diagonal([1.0, 2.0]))
- @test eigvals(A) == diag(A)
- @test eigen(A).values == diag(A)
- end
- end
-
- if VERSION >= v"1.11"
- @testset "Method ambiguity in eigen with Julia 1.11 #141" begin
- M = Hermitian(Tridiagonal(ones(ComplexF64, 2), ones(ComplexF64, 3), ones(ComplexF64, 2)))
- @test eigen(M).values ≈ Float64.(eigen(big.(M)).values)
- end
+ @testset "#139" begin
+ A = Hermitian(Diagonal([1.0, 2.0]))
+ @test GenericLinearAlgebra.eigvals(A) == LinearAlgebra.diag(A)
+ @test GenericLinearAlgebra.eigen(A).values == LinearAlgebra.diag(A)
end
end
diff --git a/test/juliaBLAS.jl b/test/juliaBLAS.jl
index 1664304..dba84dd 100644
--- a/test/juliaBLAS.jl
+++ b/test/juliaBLAS.jl
@@ -1,29 +1,30 @@
-using Test, GenericLinearAlgebra, LinearAlgebra
+using Test, GenericLinearAlgebra
+using LinearAlgebra: LinearAlgebra
@testset "rankUpdate!" begin
- A, B, x = (Hermitian(randn(5, 5)), randn(5, 2), randn(5))
+ A, B, x = (LinearAlgebra.Hermitian(randn(5, 5)), randn(5, 2), randn(5))
Ac, Bc, xc = (
- Hermitian(complex.(randn(5, 5), randn(5, 5))),
+ LinearAlgebra.Hermitian(complex.(randn(5, 5), randn(5, 5))),
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 GenericLinearAlgebra.rankUpdate!(copy(A), x) ≈ A .+ x .* x'
+ @test GenericLinearAlgebra.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 GenericLinearAlgebra.rankUpdate!(copy(A), B, 0.5, 0.5) ≈ 0.5 * A + 0.5 * B * B'
+ @test GenericLinearAlgebra.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)
+ @test GenericLinearAlgebra.rankUpdate!(copy(Ac), Bc, 1.0) ≈
+ GenericLinearAlgebra.rankUpdate!(copy(Ac), Bc, 1.0, 1.0)
end
@testset "triangular multiplication: $(typeof(T))" for T ∈ (
- UpperTriangular(complex.(randn(5, 5), randn(5, 5))),
- UnitUpperTriangular(complex.(randn(5, 5), randn(5, 5))),
- LowerTriangular(complex.(randn(5, 5), randn(5, 5))),
- UnitLowerTriangular(complex.(randn(5, 5), randn(5, 5))),
+ LinearAlgebra.UpperTriangular(complex.(randn(5, 5), randn(5, 5))),
+ LinearAlgebra.UnitUpperTriangular(complex.(randn(5, 5), randn(5, 5))),
+ LinearAlgebra.LowerTriangular(complex.(randn(5, 5), randn(5, 5))),
+ LinearAlgebra.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)
+ @test GenericLinearAlgebra.lmul!(T, copy(B), complex(0.5, 0.5)) ≈ T * B * complex(0.5, 0.5)
+ @test GenericLinearAlgebra.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 9d7ec1e..18e9150 100644
--- a/test/lapack.jl
+++ b/test/lapack.jl
@@ -1,12 +1,13 @@
-using Test, GenericLinearAlgebra, LinearAlgebra
+using Test, GenericLinearAlgebra
using GenericLinearAlgebra.LAPACK2
+using LinearAlgebra: LinearAlgebra
@testset "LAPACK wrappers" begin
n = 100
T = SymTridiagonal(randn(n), randn(n - 1))
- vals, vecs = eigen(T)
+ vals, vecs = LinearAlgebra.eigen(T)
@testset "steqr" begin
_vals, _vecs = LAPACK2.steqr!('N', copy(T.dv), copy(T.ev))
@test vals ≈ _vals
@@ -66,7 +67,7 @@ using GenericLinearAlgebra.LAPACK2
_vals = sort(LAPACK2.lahqr!(copy(T))[1])
@test _vals ≈ sort(real.(GenericLinearAlgebra._eigvals!(copy(T))))
# LAPACK's multishift algorithm (the default) seems to be broken
- @test !(_vals ≈ sort(eigvals(T)))
+ @test !(_vals ≈ sort(LinearAlgebra.eigvals(T)))
end
@testset "syevd: eltype=$eltype, uplo=$uplo" for eltype in (
@@ -84,7 +85,7 @@ using GenericLinearAlgebra.LAPACK2
else
vals, vecs = LAPACK2.heevd!('V', uplo, copy(A))
end
- @test diag(vecs' * A * vecs) ≈ eigvals(A)
+ @test LinearAlgebra.diag(vecs' * A * vecs) ≈ LinearAlgebra.eigvals(A)
end
@testset "tgevc: eltype=$eltype, side=$side, howmny=$howmny" for eltype in
@@ -93,16 +94,16 @@ using GenericLinearAlgebra.LAPACK2
howmny in ('A', 'S')
select = ones(Int, n)
- S, P = triu(randn(eltype, n, n)), triu(randn(eltype, n, n))
+ S, P = LinearAlgebra.triu(randn(eltype, n, n)), LinearAlgebra.triu(randn(eltype, n, n))
VL, VR, m = LAPACK2.tgevc!(side, howmny, select, copy(S), copy(P))
if side ∈ ('R', 'B')
- w = diag(S * VR) ./ diag(P * VR)
- @test S * VR ≈ P * VR * Diagonal(w) rtol = sqrt(eps(eltype)) atol =
+ w = LinearAlgebra.diag(S * VR) ./ LinearAlgebra.diag(P * VR)
+ @test S * VR ≈ P * VR * LinearAlgebra.Diagonal(w) rtol = sqrt(eps(eltype)) atol =
sqrt(eps(eltype))
end
if side ∈ ('L', 'B')
- w = w = diag(VL' * S) ./ diag(VL' * P)
- @test VL' * S ≈ Diagonal(w) * VL' * P rtol = sqrt(eps(eltype)) atol =
+ w = w = LinearAlgebra.diag(VL' * S) ./ LinearAlgebra.diag(VL' * P)
+ @test VL' * S ≈ LinearAlgebra.Diagonal(w) * VL' * P rtol = sqrt(eps(eltype)) atol =
sqrt(eps(eltype))
end
end
@@ -111,7 +112,7 @@ using GenericLinearAlgebra.LAPACK2
d = fill(10.0, n)
e = fill(1.0, n - 1)
vals, vecs = LAPACK2.pteqr!('I', copy(d), copy(e))
- @test SymTridiagonal(d, e) ≈ vecs * Diagonal(vals) * vecs'
+ @test LinearAlgebra.SymTridiagonal(d, e) ≈ vecs * LinearAlgebra.Diagonal(vals) * vecs'
vals2, _ = LAPACK2.pteqr!('N', copy(d), copy(e))
@test vals ≈ vals2
diff --git a/test/qr.jl b/test/qr.jl
index 1b7bbb3..4bec48c 100644
--- a/test/qr.jl
+++ b/test/qr.jl
@@ -1,7 +1,6 @@
using Test
-using LinearAlgebra
+using LinearAlgebra: LinearAlgebra
using GenericLinearAlgebra
-using GenericLinearAlgebra: qrBlocked!
@testset "The QR decomposition" begin
@testset "Problem dimension ($m,$n) with block size $bz" for (m, n) in (
@@ -15,22 +14,22 @@ using GenericLinearAlgebra: qrBlocked!
bz in (1, 2, 3, 4, 7, 8, 9, 15, 16, 17, 31, 32, 33)
A = randn(m, n)
- Aqr = qrBlocked!(copy(A), bz)
+ Aqr = GenericLinearAlgebra.qrBlocked!(copy(A), bz)
AqrQ = Aqr[Tuple{:QBlocked}]
if m >= n
@test (AqrQ'A)[1:min(m, n), :] ≈ Aqr[Tuple{:R}]
else # For type stability getindex(,Tuple{:R}) throw when the output is trapezoid
- @test (AqrQ'A) ≈ triu(Aqr.factors)
+ @test (AqrQ'A) ≈ LinearAlgebra.triu(Aqr.factors)
end
@test AqrQ' * (AqrQ * A) ≈ A
end
@testset "Error paths" begin
- @test_throws DimensionMismatch LinearAlgebra.reflectorApply!(
+ @test_throws DimensionMismatch GenericLinearAlgebra.reflectorApply!(
zeros(5, 5),
zeros(4),
1.0,
)
- @test_throws ArgumentError qrBlocked!(randn(5, 10))[Tuple{:R}]
+ @test_throws ArgumentError GenericLinearAlgebra.qrBlocked!(randn(5, 10))[Tuple{:R}]
end
end
diff --git a/test/runtests.jl b/test/runtests.jl
index 4de5615..dc2552e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,6 @@ include("cholesky.jl")
include("qr.jl")
include("eigenselfadjoint.jl")
include("eigengeneral.jl")
-include("tridiag.jl")
include("svd.jl")
include("lapack.jl")
# end
diff --git a/test/svd.jl b/test/svd.jl
index bf211b4..7567fb7 100644
--- a/test/svd.jl
+++ b/test/svd.jl
@@ -1,4 +1,5 @@
-using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
+using Test, GenericLinearAlgebra, Quaternions, DoubleFloats
+using LinearAlgebra: LinearAlgebra
@testset "Singular value decomposition" begin
@testset "Problem dimension ($m,$n)" for (m, n) in (
@@ -14,28 +15,27 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
)
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
+ U = LinearAlgebra.qr(Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:m, j = 1:min(m, n)]).Q
+ V = LinearAlgebra.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'
@test size(A) == (m, n)
- @test vals ≈ svdvals(A)
+ @test vals ≈ GenericLinearAlgebra.svdvals(A)
- F = svd(A)
+ F = GenericLinearAlgebra.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)
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)
- @test GenericLinearAlgebra._svdvals!(copy(J)) ≈ svdvals(J)
- @test GenericLinearAlgebra._svdvals!(copy(J))[end] / svdvals(J)[end] - 1 < n * eps()
+ J = LinearAlgebra.Bidiagonal(0.5 * ones(n), ones(n - 1), :U)
+ @test GenericLinearAlgebra.svdvals(J) ≈ LinearAlgebra.svdvals(J)
+ @test GenericLinearAlgebra.svdvals(J)[end] / LinearAlgebra.svdvals(J)[end] - 1 < n * eps()
end
@testset "Compare to Base methods. Problem dimension ($m,$n)" for (m, n) in (
@@ -45,57 +45,60 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
) # wide
A = randn(m, n)
- Abig = big.(A)
- @test svdvals(A) ≈ Vector{Float64}(svdvals(Abig))
- @test cond(A) ≈ Float64(cond(Abig))
+ @test GenericLinearAlgebra.svdvals(A) ≈ GenericLinearAlgebra.svdvals(A)
+ @test LinearAlgebra.cond(A) ≈ GenericLinearAlgebra.cond(A)
- F = svd(A)
- Fbig = svd(Abig)
- @test abs.(F.U'Float64.(Fbig.U)) ≈ I
- @test abs.(F.V'Float64.(Fbig.V)) ≈ I
+ F_LAPACK = LinearAlgebra.svd(A)
+ F_GLA = GenericLinearAlgebra.svd(A)
+ @test abs.(F_LAPACK.U'F_GLA.U) ≈ I
+ @test abs.(F_LAPACK.V'F_GLA.V) ≈ I
- 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
+ F_LAPACK = LinearAlgebra.svd(A, full = true)
+ F_GLA = GenericLinearAlgebra.svd(A, full = true)
+ @test abs.(F_LAPACK.U'F_GLA.U) ≈ I
+ @test abs.(F_LAPACK.V'F_GLA.V) ≈ I
end
@testset "Issue 54" begin
- U0, _, V0 = svd(big.(reshape(0:15, 4, 4)))
+ U0, _, V0 = GenericLinearAlgebra.svd(big.(reshape(0:15, 4, 4)))
A = U0[:, 1:3] * V0[:, 1:3]'
- U, S, V = svd(A)
+ U, S, V = GenericLinearAlgebra.svd(A)
@test A ≈ U * Diagonal(S) * V'
end
- @testset "Empty matrices. Issue 70" begin
- @test svdvals!(Float16.(ones(10, 0))) == Float16[]
- @test svdvals!(Float16.(ones(0, 10))) == Float16[]
- U, s, Vt = svd!(Float16.(ones(10, 0)))
- @test U == Matrix{Float16}(undef, 10, 0)
- @test s == Float16[]
- @test Vt == Matrix{Float16}(undef, 0, 0)
- U, s, Vt = svd!(Float16.(ones(0, 10)))
- @test U == Matrix{Float16}(undef, 0, 0)
- @test s == Float16[]
- @test Vt == Matrix{Float16}(undef, 10, 0)
+ @testset "Empty matrices. Issue 70. Eltype: $T" for T in (Float16, Float64)
+ @test eltype(GenericLinearAlgebra.svdvals(ones(T, 10, 0))) == T
+ @test eltype(GenericLinearAlgebra.svdvals(ones(T, 0, 10))) == T
+ if T == Float16
+ U, s, Vt = GenericLinearAlgebra.svd(ones(T, 10, 0))
+ @test U == Matrix{T}(undef, 10, 0)
+ @test eltype(s) == T
+ @test Vt == Matrix{T}(undef, 0, 0)
+ U, s, Vt = GenericLinearAlgebra.svd(ones(T, 0, 10))
+ @test U == Matrix{T}(undef, 0, 0)
+ @test eltype(s) == T
+ @test Vt == Matrix{T}(undef, 10, 0)
+ else
+ @test_broken false
+ end
end
@testset "Very small matrices. Issue 79" begin
A = randn(1, 2)
- FA = svd(A)
- 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
+ F_LAPACK = LinearAlgebra.svd(A)
+ F_GLA = GenericLinearAlgebra.svd(A)
+ Ft_GLA = GenericLinearAlgebra.svd(A')
+ @test F_LAPACK.S ≈ F_GLA.S ≈ Ft_GLA.S
+ @test abs.(F_LAPACK.U' * F_GLA.U) ≈ I
+ @test abs.(F_LAPACK.U' * Ft_GLA.V) ≈ I
+ @test abs.(F_LAPACK.V' * F_GLA.V) ≈ I
+ @test abs.(F_LAPACK.V' * Ft_GLA.U) ≈ I
end
@testset "Issue 81" begin
A = [1 0 0 0; 0 2 1 0; 0 1 2 0; 0 0 0 -1]
- @test Float64.(svdvals(big.(A))) ≈ svdvals(A)
+ @test Float64.(GenericLinearAlgebra.svdvals(big.(A))) ≈ LinearAlgebra.svdvals(A)
A = [
0.3 0.0 0.0 0.0 0.0 0.2 0.3 0.0
@@ -107,26 +110,24 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
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,
- ) ≈ svdvals(A)
+ @test GenericLinearAlgebra.svdvals(A) ≈ LinearAlgebra.svdvals(A)
n = 17
A = zeros(Double64, n, n)
for j = 1:n, i = 1:n
A[i, j] = 1 / Double64(i + j - 1)
end
- @test svdvals(A) ≈ svdvals(Float64.(A))
+ @test GenericLinearAlgebra.svdvals(A) ≈ LinearAlgebra.svdvals(Float64.(A))
# From https://github.com/JuliaMath/DoubleFloats.jl/issues/149
n = 64
c = Complex{BigFloat}(3 // 1 + 1im // 1)
- A = diagm(
+ A = LinearAlgebra.diagm(
1 => c * ones(BigFloat, n - 1),
-1 => c * ones(BigFloat, n - 1),
-2 => ones(BigFloat, n - 2),
)
- @test svdvals(A) ≈ svdvals(Complex{Double64}.(A))
+ @test GenericLinearAlgebra.svdvals(A) ≈ GenericLinearAlgebra.svdvals(Complex{Double64}.(A))
end
@testset "Issue 104. Trailing zero in bidiagonal." begin
@@ -165,9 +166,9 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
-2.169544740239464e-155,
2.3352910098001318e-232,
]
- B = Bidiagonal(dv, ev, :U)
- F = GenericLinearAlgebra._svd!(copy(B))
- @test diag(F.U' * B * F.Vt') ≈ F.S rtol = 5e-15
+ B = LinearAlgebra.Bidiagonal(dv, ev, :U)
+ F = GenericLinearAlgebra.svd(B)
+ @test LinearAlgebra.diag(F.U' * B * F.Vt') ≈ F.S rtol = 5e-15
dv = [
-2.130128478463753,
@@ -244,9 +245,9 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
-1.8291280657969477e-46,
1.9436291711597153e-50,
]
- B = Bidiagonal(dv, ev, :U)
- F = GenericLinearAlgebra._svd!(copy(B))
- @test diag(F.U' * B * F.Vt') ≈ F.S rtol = 5e-15
+ B = LinearAlgebra.Bidiagonal(dv, ev, :U)
+ F = GenericLinearAlgebra.svd(B)
+ @test LinearAlgebra.diag(F.U' * B * F.Vt') ≈ F.S rtol = 5e-15
end
@testset "Generic HessenbergQ multiplication" begin
@@ -256,14 +257,14 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
end
@testset "Issue 119" begin
- F = svd(zeros(BigFloat, 2, 2))
+ F = GenericLinearAlgebra.svd(zeros(BigFloat, 2, 2))
@test F.S == zeros(2)
@test F.U == I
@test F.Vt == I
end
@testset "Issue 121" begin
- @test svdvals(BigFloat[0 0; 1 -1]) ≈ [sqrt(2), 0]
- @test svdvals(BigFloat[1 0 0; 0 0 0; 0 1 -1]) ≈ [sqrt(2), 1, 0]
+ @test GenericLinearAlgebra.svdvals(BigFloat[0 0; 1 -1]) ≈ [sqrt(2), 0]
+ @test GenericLinearAlgebra.svdvals(BigFloat[1 0 0; 0 0 0; 0 1 -1]) ≈ [sqrt(2), 1, 0]
end
end
diff --git a/test/tridiag.jl b/test/tridiag.jl
deleted file mode 100644
index 63c3c84..0000000
--- a/test/tridiag.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-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))
-end
From e037d77a3d22d7937164468db6c6f8e3d4d5f017 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Fri, 18 Apr 2025 22:04:14 +0200
Subject: [PATCH 2/2] Fixup method extension for Hessenberg
---
src/eigenGeneral.jl | 12 +++++++-----
test/eigengeneral.jl | 4 ++--
2 files changed, 9 insertions(+), 7 deletions(-)
diff --git a/src/eigenGeneral.jl b/src/eigenGeneral.jl
index 58ce203..e5057bc 100644
--- a/src/eigenGeneral.jl
+++ b/src/eigenGeneral.jl
@@ -19,9 +19,10 @@ function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat)
lmul!(G, view(Hd, 1:n, i:n))
lmul!(G, B)
end
- ldiv!(UpperTriangular(Hd), B)
+ LinearAlgebra.ldiv!(UpperTriangular(Hd), B)
end
-(\)(H::HessenbergMatrix, B::AbstractVecOrMat) = ldiv!(copy(H), copy(B))
+LinearAlgebra.:\(H::HessenbergMatrix, B::AbstractVecOrMat) =
+ LinearAlgebra.ldiv!(copy(H), copy(B))
# Hessenberg factorization
struct HessenbergFactorization{T,S<:StridedMatrix,U} <: Factorization{T}
@@ -32,7 +33,7 @@ end
Base.copy(HF::HessenbergFactorization{T,S,U}) where {T,S,U} =
HessenbergFactorization{T,S,U}(copy(HF.data), copy(HF.τ))
-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
@@ -45,7 +46,6 @@ function _hessenberg!(A::StridedMatrix{T}) where {T}
end
return HessenbergFactorization{T,typeof(A),eltype(τ)}(A, τ)
end
-hessenberg!(A::StridedMatrix) = _hessenberg!(A)
Base.size(H::HessenbergFactorization, args...) = size(H.data, args...)
@@ -57,6 +57,8 @@ function Base.getproperty(F::HessenbergFactorization, s::Symbol)
end
end
+Base.propertynames(F::HessenbergFactorization) = (fieldnames(typeof(F))..., :H)
+
# Schur
struct Schur{T,S<:StridedMatrix} <: Factorization{T}
data::S
@@ -165,7 +167,7 @@ function _schur!(
return Schur{T,typeof(HH)}(HH, τ)
end
-schur!(A::StridedMatrix; kwargs...) = _schur!(_hessenberg!(A); kwargs...)
+schur!(A::StridedMatrix; kwargs...) = _schur!(hessenberg!(A); kwargs...)
function singleShiftQR!(
HH::StridedMatrix,
diff --git a/test/eigengeneral.jl b/test/eigengeneral.jl
index c13e94a..443fd8c 100644
--- a/test/eigengeneral.jl
+++ b/test/eigengeneral.jl
@@ -237,10 +237,10 @@ using LinearAlgebra: LinearAlgebra
@test sort(imag(vals)) ≈ sort(imag(λs)) atol = 1e-25
end
- @testset "_hessenberg! and Hessenberg" begin
+ @testset "hessenberg! and Hessenberg" begin
n = 10
A = randn(n, n)
- HF = GenericLinearAlgebra._hessenberg!(copy(A))
+ HF = GenericLinearAlgebra.hessenberg!(copy(A))
for i = 1:length(HF.τ)
HM = convert(Matrix, HF.τ[i])
A[(i+1):end, :] = HM * A[(i+1):end, :]