Skip to content

Commit 4e7c3f4

Browse files
dkarraschdkarrasch
authored
Backports release v1.12 (#1237)
2 parents f0f7a46 + e909af2 commit 4e7c3f4

11 files changed

+83
-57
lines changed

Diff for: ‎src/bidiag.jl

+4-4
Original file line numberDiff line numberDiff line change
@@ -1183,25 +1183,25 @@ function _dibimul!(C::Bidiagonal, A::Diagonal, B::Bidiagonal, _add)
11831183
C
11841184
end
11851185

1186-
function *(A::UpperOrUnitUpperTriangular, B::Bidiagonal)
1186+
function mul(A::UpperOrUnitUpperTriangular, B::Bidiagonal)
11871187
TS = promote_op(matprod, eltype(A), eltype(B))
11881188
C = mul!(similar(A, TS, size(A)), A, B)
11891189
return B.uplo == 'U' ? UpperTriangular(C) : C
11901190
end
11911191

1192-
function *(A::LowerOrUnitLowerTriangular, B::Bidiagonal)
1192+
function mul(A::LowerOrUnitLowerTriangular, B::Bidiagonal)
11931193
TS = promote_op(matprod, eltype(A), eltype(B))
11941194
C = mul!(similar(A, TS, size(A)), A, B)
11951195
return B.uplo == 'L' ? LowerTriangular(C) : C
11961196
end
11971197

1198-
function *(A::Bidiagonal, B::UpperOrUnitUpperTriangular)
1198+
function mul(A::Bidiagonal, B::UpperOrUnitUpperTriangular)
11991199
TS = promote_op(matprod, eltype(A), eltype(B))
12001200
C = mul!(similar(B, TS, size(B)), A, B)
12011201
return A.uplo == 'U' ? UpperTriangular(C) : C
12021202
end
12031203

1204-
function *(A::Bidiagonal, B::LowerOrUnitLowerTriangular)
1204+
function mul(A::Bidiagonal, B::LowerOrUnitLowerTriangular)
12051205
TS = promote_op(matprod, eltype(A), eltype(B))
12061206
C = mul!(similar(B, TS, size(B)), A, B)
12071207
return A.uplo == 'L' ? LowerTriangular(C) : C

Diff for: ‎src/diagonal.jl

+7-17
Original file line numberDiff line numberDiff line change
@@ -322,7 +322,7 @@ Base.literal_pow(::typeof(^), D::Diagonal, valp::Val) =
322322
Diagonal(Base.literal_pow.(^, D.diag, valp)) # for speed
323323
Base.literal_pow(::typeof(^), D::Diagonal, ::Val{-1}) = inv(D) # for disambiguation
324324

325-
function (*)(Da::Diagonal, Db::Diagonal)
325+
function mul(Da::Diagonal, Db::Diagonal)
326326
matmul_size_check(size(Da), size(Db))
327327
return Diagonal(Da.diag .* Db.diag)
328328
end
@@ -332,26 +332,19 @@ function (*)(D::Diagonal, V::AbstractVector)
332332
return D.diag .* V
333333
end
334334

335-
function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal)
335+
function mul(A::AdjOrTransAbsMat, D::Diagonal)
336336
adj = wrapperop(A)
337337
copy(adj(adj(D) * adj(A)))
338338
end
339-
function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number})
340-
@invoke *(A::AbstractMatrix, D::AbstractMatrix)
339+
function mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number})
340+
@invoke mul(A::AbstractMatrix, D::AbstractMatrix)
341341
end
342-
function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat)
342+
function mul(D::Diagonal, A::AdjOrTransAbsMat)
343343
adj = wrapperop(A)
344344
copy(adj(adj(A) * adj(D)))
345345
end
346-
function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix})
347-
@invoke *(D::AbstractMatrix, A::AbstractMatrix)
348-
end
349-
350-
function (*)(A::AdjOrTransAbsMat, D::Diagonal)
351-
_diag_adj_mul(A, D)
352-
end
353-
function (*)(D::Diagonal, A::AdjOrTransAbsMat)
354-
_diag_adj_mul(D, A)
346+
function mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix})
347+
@invoke mul(D::AbstractMatrix, A::AbstractMatrix)
355348
end
356349

357350
function rmul!(A::AbstractMatrix, D::Diagonal)
@@ -1037,9 +1030,6 @@ end
10371030
*(x::TransposeAbsVec, D::Diagonal, y::AbstractVector) = _mapreduce_prod(*, x, D, y)
10381031
/(u::AdjointAbsVec, D::Diagonal) = (D' \ u')'
10391032
/(u::TransposeAbsVec, D::Diagonal) = transpose(transpose(D) \ transpose(u))
1040-
# disambiguation methods: Call unoptimized version for user defined AbstractTriangular.
1041-
*(A::AbstractTriangular, D::Diagonal) = @invoke *(A::AbstractMatrix, D::Diagonal)
1042-
*(D::Diagonal, A::AbstractTriangular) = @invoke *(D::Diagonal, A::AbstractMatrix)
10431033

10441034
_opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag)
10451035
_opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag)

Diff for: ‎src/hessenberg.jl

+7-3
Original file line numberDiff line numberDiff line change
@@ -137,7 +137,7 @@ for T = (:UniformScaling, :Diagonal, :Bidiagonal, :Tridiagonal, :SymTridiagonal,
137137
end
138138
end
139139

140-
for T = (:Number, :UniformScaling, :Diagonal)
140+
for T = (:Number, :UniformScaling)
141141
@eval begin
142142
*(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data * x)
143143
*(x::$T, H::UpperHessenberg) = UpperHessenberg(x * H.data)
@@ -146,19 +146,23 @@ for T = (:Number, :UniformScaling, :Diagonal)
146146
end
147147
end
148148

149-
function *(H::UpperHessenberg, U::UpperOrUnitUpperTriangular)
149+
mul(H::UpperHessenberg, D::Diagonal) = UpperHessenberg(H.data * D)
150+
mul(D::Diagonal, H::UpperHessenberg) = UpperHessenberg(D * H.data)
151+
function mul(H::UpperHessenberg, U::UpperOrUnitUpperTriangular)
150152
HH = mul!(matprod_dest(H, U, promote_op(matprod, eltype(H), eltype(U))), H, U)
151153
UpperHessenberg(HH)
152154
end
153-
function *(U::UpperOrUnitUpperTriangular, H::UpperHessenberg)
155+
function mul(U::UpperOrUnitUpperTriangular, H::UpperHessenberg)
154156
HH = mul!(matprod_dest(U, H, promote_op(matprod, eltype(U), eltype(H))), U, H)
155157
UpperHessenberg(HH)
156158
end
157159

160+
/(H::UpperHessenberg, D::Diagonal) = UpperHessenberg(H.data / D)
158161
function /(H::UpperHessenberg, U::UpperTriangular)
159162
HH = _rdiv!(matprod_dest(H, U, promote_op(/, eltype(H), eltype(U))), H, U)
160163
UpperHessenberg(HH)
161164
end
165+
\(D::Diagonal, H::UpperHessenberg) = UpperHessenberg(D \ H.data)
162166
function /(H::UpperHessenberg, U::UnitUpperTriangular)
163167
HH = _rdiv!(matprod_dest(H, U, promote_op(/, eltype(H), eltype(U))), H, U)
164168
UpperHessenberg(HH)

Diff for: ‎src/special.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -120,12 +120,12 @@ _mul!(C::AbstractMatrix, A::AbstractTriangular, B::BandedMatrix, alpha::Number,
120120
_mul!(C::AbstractMatrix, A::BandedMatrix, B::AbstractTriangular, alpha::Number, beta::Number) =
121121
@stable_muladdmul _mul!(C, A, B, MulAddMul(alpha, beta))
122122

123-
function *(H::UpperHessenberg, B::Bidiagonal)
123+
function mul(H::UpperHessenberg, B::Bidiagonal)
124124
T = promote_op(matprod, eltype(H), eltype(B))
125125
A = mul!(similar(H, T, size(H)), H, B)
126126
return B.uplo == 'U' ? UpperHessenberg(A) : A
127127
end
128-
function *(B::Bidiagonal, H::UpperHessenberg)
128+
function mul(B::Bidiagonal, H::UpperHessenberg)
129129
T = promote_op(matprod, eltype(B), eltype(H))
130130
A = mul!(similar(H, T, size(H)), B, H)
131131
return B.uplo == 'U' ? UpperHessenberg(A) : A

Diff for: ‎src/symmetric.jl

+14-1
Original file line numberDiff line numberDiff line change
@@ -702,7 +702,20 @@ for f in (:+, :-)
702702
end
703703
end
704704

705-
*(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B)
705+
mul(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B)
706+
# catch a few potential BLAS-cases
707+
function mul(A::HermOrSym{<:BlasFloat,<:StridedMatrix}, B::AdjOrTrans{<:BlasFloat,<:StridedMatrix})
708+
T = promote_type(eltype(A), eltype(B))
709+
mul!(similar(B, T, (size(A, 1), size(B, 2))),
710+
convert(AbstractMatrix{T}, A),
711+
copy_oftype(B, T)) # make sure the AdjOrTrans wrapper is resolved
712+
end
713+
function mul(A::AdjOrTrans{<:BlasFloat,<:StridedMatrix}, B::HermOrSym{<:BlasFloat,<:StridedMatrix})
714+
T = promote_type(eltype(A), eltype(B))
715+
mul!(similar(B, T, (size(A, 1), size(B, 2))),
716+
copy_oftype(A, T), # make sure the AdjOrTrans wrapper is resolved
717+
convert(AbstractMatrix{T}, B))
718+
end
706719

707720
function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
708721
require_one_based_indexing(x, y)

Diff for: ‎src/symmetriceigen.jl

+19-13
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,16 @@ eigencopy_oftype(A::Hermitian, S) = Hermitian(copytrito!(similar(parent(A), S, s
66
eigencopy_oftype(A::Symmetric, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo))
77
eigencopy_oftype(A::Symmetric{<:Complex}, S) = copyto!(similar(parent(A), S), A)
88

9-
default_eigen_alg(A) = DivideAndConquer()
9+
"""
10+
default_eigen_alg(A)
11+
12+
Return the default algorithm used to solve the eigensystem `A v = λ v` for a symmetric matrix `A`.
13+
Defaults to `LinearAlegbra.DivideAndConquer()`, which corresponds to the LAPACK function `LAPACK.syevd!`.
14+
"""
15+
default_eigen_alg(@nospecialize(A)) = DivideAndConquer()
1016

1117
# Eigensolvers for symmetric and Hermitian matrices
12-
function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
18+
function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
1319
if alg === DivideAndConquer()
1420
Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...)
1521
elseif alg === QRIteration()
@@ -22,7 +28,7 @@ function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algo
2228
end
2329

2430
"""
25-
eigen(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A)) -> Eigen
31+
eigen(A::Union{Hermitian, Symmetric}; alg::LinearAlgebra.Algorithm = LinearAlgebra.default_eigen_alg(A)) -> Eigen
2632
2733
Compute the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F`
2834
which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
@@ -45,19 +51,19 @@ The default `alg` used may change in the future.
4551
4652
The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).
4753
"""
48-
function eigen(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
49-
_eigen(A, alg; sortby)
54+
function eigen(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
55+
_eigen(A; alg, sortby)
5056
end
5157

5258
# we dispatch on the eltype in an internal method to avoid ambiguities
53-
function _eigen(A::RealHermSymComplexHerm, alg::Algorithm; sortby)
59+
function _eigen(A::RealHermSymComplexHerm; alg::Algorithm, sortby)
5460
S = eigtype(eltype(A))
55-
eigen!(eigencopy_oftype(A, S), alg; sortby)
61+
eigen!(eigencopy_oftype(A, S); alg, sortby)
5662
end
5763

58-
function _eigen(A::RealHermSymComplexHerm{Float16}, alg::Algorithm; sortby::Union{Function,Nothing}=nothing)
64+
function _eigen(A::RealHermSymComplexHerm{Float16}; alg::Algorithm, sortby::Union{Function,Nothing}=nothing)
5965
S = eigtype(eltype(A))
60-
E = eigen!(eigencopy_oftype(A, S), alg, sortby=sortby)
66+
E = eigen!(eigencopy_oftype(A, S); alg, sortby)
6167
values = convert(AbstractVector{Float16}, E.values)
6268
vectors = convert(AbstractMatrix{isreal(E.vectors) ? Float16 : Complex{Float16}}, E.vectors)
6369
return Eigen(values, vectors)
@@ -114,7 +120,7 @@ function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
114120
end
115121

116122

117-
function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
123+
function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
118124
vals::Vector{real(eltype(A))} = if alg === DivideAndConquer()
119125
LAPACK.syevd!('N', A.uplo, A.data)
120126
elseif alg === QRIteration()
@@ -129,7 +135,7 @@ function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Al
129135
end
130136

131137
"""
132-
eigvals(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A))) -> values
138+
eigvals(A::Union{Hermitian, Symmetric}; alg::Algorithm = default_eigen_alg(A))) -> values
133139
134140
Return the eigenvalues of `A`.
135141
@@ -143,9 +149,9 @@ a comparison of the accuracy and performance of different methods.
143149
144150
The default `alg` used may change in the future.
145151
"""
146-
function eigvals(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
152+
function eigvals(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
147153
S = eigtype(eltype(A))
148-
eigvals!(eigencopy_oftype(A, S), alg; sortby)
154+
eigvals!(eigencopy_oftype(A, S); alg, sortby)
149155
end
150156

151157

Diff for: ‎src/triangular.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -1886,7 +1886,7 @@ end
18861886

18871887
## Some Triangular-Triangular cases. We might want to write tailored methods
18881888
## for these cases, but I'm not sure it is worth it.
1889-
for f in (:*, :\)
1889+
for f in (:mul, :\)
18901890
@eval begin
18911891
($f)(A::LowerTriangular, B::LowerTriangular) =
18921892
LowerTriangular(@invoke $f(A::LowerTriangular, B::AbstractMatrix))

Diff for: ‎test/matmul.jl

+8
Original file line numberDiff line numberDiff line change
@@ -325,6 +325,10 @@ end
325325
@test 0 == @allocations mul!(C, At, Bt)
326326
end
327327
# syrk/herk
328+
mul!(C, transpose(A), A)
329+
mul!(C, adjoint(A), A)
330+
mul!(C, A, transpose(A))
331+
mul!(C, A, adjoint(A))
328332
@test 0 == @allocations mul!(C, transpose(A), A)
329333
@test 0 == @allocations mul!(C, adjoint(A), A)
330334
@test 0 == @allocations mul!(C, A, transpose(A))
@@ -334,6 +338,7 @@ end
334338
Ac = complex(A)
335339
for t in (identity, adjoint, transpose)
336340
Bt = t(B)
341+
mul!(Cc, Ac, Bt)
337342
@test 0 == @allocations mul!(Cc, Ac, Bt)
338343
end
339344
end
@@ -356,6 +361,9 @@ end
356361
A = rand(-10:10, n, n)
357362
B = ones(Float64, n, n)
358363
C = zeros(Float64, n, n)
364+
mul!(C, A, B)
365+
mul!(C, A, transpose(B))
366+
mul!(C, adjoint(A), B)
359367
@test 0 == @allocations mul!(C, A, B)
360368
@test 0 == @allocations mul!(C, A, transpose(B))
361369
@test 0 == @allocations mul!(C, adjoint(A), B)

Diff for: ‎test/symmetric.jl

-15
Original file line numberDiff line numberDiff line change
@@ -720,21 +720,6 @@ end
720720
end
721721
end
722722

723-
@testset "eigendecomposition Algorithms" begin
724-
using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations
725-
for T in (Float64, ComplexF64, Float32, ComplexF32)
726-
n = 4
727-
A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n))
728-
d, v = eigen(A)
729-
for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations())
730-
@test (@inferred eigvals(A, alg)) ≈ d
731-
d2, v2 = @inferred eigen(A, alg)
732-
@test d2 ≈ d
733-
@test A * v2 ≈ v2 * Diagonal(d2)
734-
end
735-
end
736-
end
737-
738723
const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
739724
isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl"))
740725
using .Main.ImmutableArrays

Diff for: ‎test/symmetriceigen.jl

+16-1
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
module TestSymmetricEigen
44

55
using Test, LinearAlgebra
6+
using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations
67

78
@testset "chol-eigen-eigvals" begin
89
## Cholesky decomposition based
@@ -173,7 +174,7 @@ end
173174
@test D.vectors ≈ D32.vectors
174175

175176
# ensure that different algorithms dispatch correctly
176-
λ, V = eigen(C, LinearAlgebra.QRIteration())
177+
λ, V = eigen(C; alg=LinearAlgebra.QRIteration())
177178
@test λ isa Vector{Float16}
178179
@test C * V ≈ V * Diagonal(λ)
179180
end
@@ -184,4 +185,18 @@ end
184185
@test S * v ≈ v * Diagonal(λ)
185186
end
186187

188+
@testset "eigendecomposition Algorithms" begin
189+
for T in (Float64, ComplexF64, Float32, ComplexF32)
190+
n = 4
191+
A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n))
192+
d, v = eigen(A)
193+
for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations())
194+
@test (@inferred eigvals(A; alg)) ≈ d
195+
d2, v2 = @inferred eigen(A; alg)
196+
@test d2 ≈ d
197+
@test A * v2 ≈ v2 * Diagonal(d2)
198+
end
199+
end
200+
end
201+
187202
end # module TestSymmetricEigen

Diff for: ‎test/triangular.jl

+5
Original file line numberDiff line numberDiff line change
@@ -703,6 +703,7 @@ end
703703
@testset "(l/r)mul! and (l/r)div! for generic triangular" begin
704704
@testset for T in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular)
705705
M = MyTriangular(T(rand(4,4)))
706+
D = Diagonal(randn(4))
706707
A = rand(4,4)
707708
Ac = similar(A)
708709
@testset "lmul!" begin
@@ -725,6 +726,10 @@ end
725726
rdiv!(Ac, M)
726727
@test Ac ≈ A / M
727728
end
729+
@testset "diagonal mul" begin
730+
@test D * M ≈ D * M.data
731+
@test M * D ≈ M.data * D
732+
end
728733
end
729734
end
730735

0 commit comments

Comments
 (0)