From c70a35fdc974854139b8e85ab43758a19f2e170f Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Thu, 25 Apr 2019 15:01:17 -0700 Subject: [PATCH 01/18] Add matrix symmetrizing functions This adds functions `symmetrize` and `symmetrize!`, which produce symmetric matrices based on the input `X` via efficient computation of `(X + X') / 2`. --- NEWS.md | 1 + stdlib/LinearAlgebra/docs/src/index.md | 2 + stdlib/LinearAlgebra/src/LinearAlgebra.jl | 2 + stdlib/LinearAlgebra/src/symmetric.jl | 46 ++++++++++++++++ stdlib/LinearAlgebra/test/symmetric.jl | 67 +++++++++++++++++++++++ 5 files changed, 118 insertions(+) diff --git a/NEWS.md b/NEWS.md index fc0fd55d24f8f..0feb55d082d0c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -90,6 +90,7 @@ Standard library changes #### Package Manager #### LinearAlgebra +* New functions `symmetrize` and `symmetrize!` for constructing symmetric matrices ([#31836]). #### Markdown diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 38c48bfe6d8d2..8ffb5bd01b459 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -460,6 +460,8 @@ Base.copy(::Union{Transpose,Adjoint}) LinearAlgebra.stride1 LinearAlgebra.checksquare LinearAlgebra.peakflops +LinearAlgebra.symmetrize +LinearAlgebra.symmetrize! ``` ## Low-level matrix operations diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index ae8b9d461ffc2..53e0f298cdb54 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -138,6 +138,8 @@ export svdvals!, svdvals, sylvester, + symmetrize, + symmetrize!, tr, transpose, transpose!, diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 7347dd6f78639..143a416af2185 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -854,3 +854,49 @@ for func in (:log, :sqrt) end end end + +""" + symmetrize!(X::AbstractMatrix; conjugate::Bool=true) + +Make the matrix `X` symmetric in-place using the formula `(X + X') / 2`. +If `conjugate` is `true`, `X` is made Hermitian rather than symmetric. + +See also: [`symmetrize`](@ref) + +!!! compat "Julia 1.3" + This function requires Julia 1.3 or later. +""" +symmetrize!(X::AbstractMatrix; conjugate::Bool=true) = + conjugate ? _symmetrize!(conj, real, X) : _symmetrize!(identity, identity, X) + +function _symmetrize!(f::Function, g::Function, X::AbstractMatrix) + inds = axes(X, 1) + inds == axes(X, 2) || throw(DimensionMismatch("matrix is not square")) + r = first(inds) + s = step(inds) + @inbounds for j in inds + X[j,j] = g(X[j,j]) + for i in r:s:j-s + X[i,j] = (X[i,j] + f(X[j,i])) / 2 + X[j,i] = f(X[i,j]) + end + end + X +end + +""" + symmetrize(X::AbstractMatrix; conjugate::Bool=true) + +Construct a symmetric matrix based on `X` using the formula `(X + X') / 2`. +If `conjugate` is `true`, the result is Hermitian rather than symmetric. + +See also: [`symmetrize!`](@ref) + +!!! compat "Julia 1.3" + This function requires Julia 1.3 or later. +""" +symmetrize(X::AbstractMatrix{T}; conjugate::Bool=true) where {T<:Number} = + symmetrize!(copyto!(similar(X, typeof(zero(T) / 2)), X), conjugate=conjugate) +symmetrize(X::Symmetric{<:Real}; conjugate::Bool=true) = X +symmetrize(X::Symmetric{<:Complex}; conjugate::Bool=true) = conjugate ? Hermitian(X) : X +symmetrize(X::Hermitian; conjugate::Bool=true) = X diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 169dfb0071718..234703a745b00 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -4,6 +4,12 @@ module TestSymmetric using Test, LinearAlgebra, SparseArrays, Random +const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") +if !isdefined(Main, :OffsetArrays) + @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl")) +end +using .Main.OffsetArrays + Random.seed!(1010) @testset "Pauli σ-matrices: $σ" for σ in map(Hermitian, @@ -776,4 +782,65 @@ end end end +@testset "symmetrize" begin + @testset "in-place" begin + for T in [Float32, Complex{Float32}], conj in [true, false], offset in [true, false] + X = rand(T, 4, 4) + offset && (X = OffsetArray(X, -2, -2)) + Y = copy(X) + symmetrize!(X, conjugate=conj) + if conj + @test X ≈ (Y + Y') / 2 + @test ishermitian(X) + else + @test X ≈ (Y + transpose(Y)) / 2 + @test issymmetric(X) + end + end + end + @testset "out-of-place" begin + for conj in [true, false] + for T in [Float32, Complex{Float32}], offset in [true, false] + X = randn(T, 4, 4) + offset && (X = OffsetArray(X, -2, -2)) + Y = symmetrize(X, conjugate=conj) + @test eltype(Y) === T + if conj + @test Y ≈ (X + X') / 2 + @test ishermitian(Y) + else + @test Y ≈ (X + transpose(X)) / 2 + @test issymmetric(Y) + end + end + for T in [Int32, Rational{Int32}, Complex{Int32}, Complex{Rational{Int32}}] + X = T[1 2 3; 4 5 6; 7 8 9] + Y = symmetrize(X, conjugate=conj) + @test Y ≈ T[1 3 5; 3 5 7; 5 7 9] + end + end + end + @testset "Symmetric and Hermitian" begin + S = Symmetric([1 2; 2 1]) + @test symmetrize(S) === S + @test symmetrize(S, conjugate=false) === S + SC = Symmetric([1+1im 2+2im; 2+2im 1+1im]) + @test symmetrize(SC) == Hermitian(SC) + @test symmetrize(SC, conjugate=false) === SC + H = Hermitian([1 2; 2 1]) + @test symmetrize(H) === H + @test symmetrize(H, conjugate=false) === H + HC = Hermitian([1+1im 2+2im; 2+2im 1+1im]) + @test symmetrize(HC) === HC + @test symmetrize(HC, conjugate=false) === HC + end + @testset "error conditions" begin + # Not square + @test_throws DimensionMismatch symmetrize!([1.0 2.0 3.0; 4.0 5.0 6.0]) + # Mismatched indexing + bad = OffsetArray([1.0 2.0; 2.0 1.0], 0, 1) + @test_throws DimensionMismatch symmetrize!(bad) + end +end + end # module TestSymmetric From 099b93a37189fa08ea6cfd4a265431660f6e1a4c Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Thu, 11 Nov 2021 20:23:08 -0800 Subject: [PATCH 02/18] Rename new functions to `(symmetric|hermitian)part(!)` --- NEWS.md | 2 +- stdlib/LinearAlgebra/docs/src/index.md | 6 +- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 6 +- stdlib/LinearAlgebra/src/symmetric.jl | 83 +++++++++++++---------- stdlib/LinearAlgebra/test/symmetric.jl | 83 ++++++----------------- 5 files changed, 79 insertions(+), 101 deletions(-) diff --git a/NEWS.md b/NEWS.md index 0feb55d082d0c..2c77f76d4e0d9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -90,7 +90,7 @@ Standard library changes #### Package Manager #### LinearAlgebra -* New functions `symmetrize` and `symmetrize!` for constructing symmetric matrices ([#31836]). +* New functions `symmetricpart`, `symmetricpart!`, `hermitianpart`, and `hermitianpart!` for extracting the symmetric or Hermitian part of a matrix ([#31836]). #### Markdown diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 8ffb5bd01b459..67af579153663 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -460,8 +460,10 @@ Base.copy(::Union{Transpose,Adjoint}) LinearAlgebra.stride1 LinearAlgebra.checksquare LinearAlgebra.peakflops -LinearAlgebra.symmetrize -LinearAlgebra.symmetrize! +LinearAlgebra.symmetricpart +LinearAlgebra.symmetricpart! +LinearAlgebra.hermitianpart +LinearAlgebra.hermitianpart! ``` ## Low-level matrix operations diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 53e0f298cdb54..c663df9f10840 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -89,6 +89,8 @@ export eigvecs, factorize, givens, + hermitianpart, + hermitianpart!, hessenberg, hessenberg!, isdiag, @@ -138,8 +140,8 @@ export svdvals!, svdvals, sylvester, - symmetrize, - symmetrize!, + symmetricpart, + symmetricpart!, tr, transpose, transpose!, diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 143a416af2185..0b0c24abdc56c 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -856,47 +856,62 @@ for func in (:log, :sqrt) end """ - symmetrize!(X::AbstractMatrix; conjugate::Bool=true) + symmetricpart(A, uplo=:U) -Make the matrix `X` symmetric in-place using the formula `(X + X') / 2`. -If `conjugate` is `true`, `X` is made Hermitian rather than symmetric. +Return the symmetric part of the square matrix `A`, defined as `(A + transpose(A)) / 2`, +as a [`Symmetric`](@ref) matrix. -See also: [`symmetrize`](@ref) +!!! compat "Julia 1.8" + This function requires Julia 1.8 or later. +""" +symmetricpart(A::AbstractMatrix{T}, uplo=:U) where {T} = + symmetricpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) -!!! compat "Julia 1.3" - This function requires Julia 1.3 or later. """ -symmetrize!(X::AbstractMatrix; conjugate::Bool=true) = - conjugate ? _symmetrize!(conj, real, X) : _symmetrize!(identity, identity, X) - -function _symmetrize!(f::Function, g::Function, X::AbstractMatrix) - inds = axes(X, 1) - inds == axes(X, 2) || throw(DimensionMismatch("matrix is not square")) - r = first(inds) - s = step(inds) - @inbounds for j in inds - X[j,j] = g(X[j,j]) - for i in r:s:j-s - X[i,j] = (X[i,j] + f(X[j,i])) / 2 - X[j,i] = f(X[i,j]) - end - end - X -end + symmetricpart!(A, uplo=:U) + +Overwrite the square matrix `A` with its symmetric part, `(A + transpose(A)) / 2`, +and return `Symmetric(A, uplo)`. +!!! compat "Julia 1.8" + This function requires Julia 1.8 or later. """ - symmetrize(X::AbstractMatrix; conjugate::Bool=true) +symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) = + Symmetric(_hermorsympart!(identity, identity, A, char_uplo(uplo)), uplo) -Construct a symmetric matrix based on `X` using the formula `(X + X') / 2`. -If `conjugate` is `true`, the result is Hermitian rather than symmetric. +""" + hermitianpart(A, uplo=:U) + +Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, +as a [`Hermitian`](@ref) matrix. -See also: [`symmetrize!`](@ref) +!!! compat "Julia 1.8" + This function requires Julia 1.8 or later. +""" +hermitianpart(A::AbstractMatrix{T}, uplo=:U) where {T} = + hermitianpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) -!!! compat "Julia 1.3" - This function requires Julia 1.3 or later. """ -symmetrize(X::AbstractMatrix{T}; conjugate::Bool=true) where {T<:Number} = - symmetrize!(copyto!(similar(X, typeof(zero(T) / 2)), X), conjugate=conjugate) -symmetrize(X::Symmetric{<:Real}; conjugate::Bool=true) = X -symmetrize(X::Symmetric{<:Complex}; conjugate::Bool=true) = conjugate ? Hermitian(X) : X -symmetrize(X::Hermitian; conjugate::Bool=true) = X + hermitianpart!(A, uplo=:U) + +Overwrite the square matrix `A` with its Hermitian part, `(A + A') / 2`, +and return `Hermitian(A, uplo)`. + +!!! compat "Julia 1.8" + This function requires Julia 1.8 or later. +""" +hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = + Hermitian(_hermorsympart!(real, conj, A, char_uplo(uplo)), uplo) + +function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix, uplo::Char) + require_one_based_indexing(A) + n = checksquare(A) + triangle = uplo === 'U' ? (j -> 1:(j - 1)) : (j -> (j + 1):n) + @inbounds for j in 1:n + A[j, j] = real(A[j, j]) + for i in triangle(j) + A[i, j] = (A[i, j] + conj(A[j, i])) / 2 + end + end + return A +end diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 234703a745b00..91049adf65a08 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -4,12 +4,6 @@ module TestSymmetric using Test, LinearAlgebra, SparseArrays, Random -const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") -if !isdefined(Main, :OffsetArrays) - @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl")) -end -using .Main.OffsetArrays - Random.seed!(1010) @testset "Pauli σ-matrices: $σ" for σ in map(Hermitian, @@ -782,65 +776,30 @@ end end end -@testset "symmetrize" begin - @testset "in-place" begin - for T in [Float32, Complex{Float32}], conj in [true, false], offset in [true, false] - X = rand(T, 4, 4) - offset && (X = OffsetArray(X, -2, -2)) - Y = copy(X) - symmetrize!(X, conjugate=conj) - if conj - @test X ≈ (Y + Y') / 2 - @test ishermitian(X) - else - @test X ≈ (Y + transpose(Y)) / 2 - @test issymmetric(X) +@testset "symmetric and hermitian parts" begin + for T in [Float32, Complex{Float32}, Int32, Rational{Int32}, + Complex{Int32}, Complex{Rational{Int32}}] + for (f, f!, t) in [(hermitianpart, hermitianpart!, adjoint), + (symmetricpart, symmetricpart!, transpose)] + X = T[1 2 3; 4 5 6; 7 8 9] + T <: Complex && (X .+= im .* X) + Xc = copy(X) + Y = (X + t(X)) / 2 + U = f(X) + L = f(X, :L) + @test U == L == Y + if T <: AbstractFloat || real(T) <: AbstractFloat + HU = f!(X) + @test HU == Y + @test triu(X) == triu(Y) + HL = f!(Xc, :L) + @test HL == Y + @test tril(Xc) == tril(Y) end end end - @testset "out-of-place" begin - for conj in [true, false] - for T in [Float32, Complex{Float32}], offset in [true, false] - X = randn(T, 4, 4) - offset && (X = OffsetArray(X, -2, -2)) - Y = symmetrize(X, conjugate=conj) - @test eltype(Y) === T - if conj - @test Y ≈ (X + X') / 2 - @test ishermitian(Y) - else - @test Y ≈ (X + transpose(X)) / 2 - @test issymmetric(Y) - end - end - for T in [Int32, Rational{Int32}, Complex{Int32}, Complex{Rational{Int32}}] - X = T[1 2 3; 4 5 6; 7 8 9] - Y = symmetrize(X, conjugate=conj) - @test Y ≈ T[1 3 5; 3 5 7; 5 7 9] - end - end - end - @testset "Symmetric and Hermitian" begin - S = Symmetric([1 2; 2 1]) - @test symmetrize(S) === S - @test symmetrize(S, conjugate=false) === S - SC = Symmetric([1+1im 2+2im; 2+2im 1+1im]) - @test symmetrize(SC) == Hermitian(SC) - @test symmetrize(SC, conjugate=false) === SC - H = Hermitian([1 2; 2 1]) - @test symmetrize(H) === H - @test symmetrize(H, conjugate=false) === H - HC = Hermitian([1+1im 2+2im; 2+2im 1+1im]) - @test symmetrize(HC) === HC - @test symmetrize(HC, conjugate=false) === HC - end - @testset "error conditions" begin - # Not square - @test_throws DimensionMismatch symmetrize!([1.0 2.0 3.0; 4.0 5.0 6.0]) - # Mismatched indexing - bad = OffsetArray([1.0 2.0; 2.0 1.0], 0, 1) - @test_throws DimensionMismatch symmetrize!(bad) - end + @test_throws DimensionMismatch symmetricpart([1.0 2.0 3.0; 4.0 5.0 6.0]) + @test_throws DimensionMismatch hermitianpart([1.0 2.0 3.0; 4.0 5.0 6.0]) end end # module TestSymmetric From ca1ab7b487dbb62d9da4ddf4034f3e1cf4a3c45a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 18 Jan 2023 16:12:17 +0100 Subject: [PATCH 03/18] Apply suggestions from code review --- stdlib/LinearAlgebra/src/symmetric.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 0b0c24abdc56c..13aabb2976a7e 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -861,8 +861,8 @@ end Return the symmetric part of the square matrix `A`, defined as `(A + transpose(A)) / 2`, as a [`Symmetric`](@ref) matrix. -!!! compat "Julia 1.8" - This function requires Julia 1.8 or later. +!!! compat "Julia 1.10" + This function requires Julia 1.10 or later. """ symmetricpart(A::AbstractMatrix{T}, uplo=:U) where {T} = symmetricpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) @@ -873,8 +873,8 @@ symmetricpart(A::AbstractMatrix{T}, uplo=:U) where {T} = Overwrite the square matrix `A` with its symmetric part, `(A + transpose(A)) / 2`, and return `Symmetric(A, uplo)`. -!!! compat "Julia 1.8" - This function requires Julia 1.8 or later. +!!! compat "Julia 1.10" + This function requires Julia 1.10 or later. """ symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) = Symmetric(_hermorsympart!(identity, identity, A, char_uplo(uplo)), uplo) @@ -885,8 +885,8 @@ symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) = Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, as a [`Hermitian`](@ref) matrix. -!!! compat "Julia 1.8" - This function requires Julia 1.8 or later. +!!! compat "Julia 1.10" + This function requires Julia 1.10 or later. """ hermitianpart(A::AbstractMatrix{T}, uplo=:U) where {T} = hermitianpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) @@ -897,8 +897,8 @@ hermitianpart(A::AbstractMatrix{T}, uplo=:U) where {T} = Overwrite the square matrix `A` with its Hermitian part, `(A + A') / 2`, and return `Hermitian(A, uplo)`. -!!! compat "Julia 1.8" - This function requires Julia 1.8 or later. +!!! compat "Julia 1.10" + This function requires Julia 1.10 or later. """ hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermorsympart!(real, conj, A, char_uplo(uplo)), uplo) From a44597f8514dea6e603f8321d76bb4309859fd72 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 18 Jan 2023 16:35:59 +0100 Subject: [PATCH 04/18] generalize to matrix of matrices case --- stdlib/LinearAlgebra/src/symmetric.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 809fcaa208464..e9bcd1f9760c5 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -868,6 +868,7 @@ as a [`Symmetric`](@ref) matrix. """ symmetricpart(A::AbstractMatrix{T}, uplo=:U) where {T} = symmetricpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) +symmetricpart(x::Number) = x """ symmetricpart!(A, uplo=:U) @@ -879,7 +880,7 @@ and return `Symmetric(A, uplo)`. This function requires Julia 1.10 or later. """ symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) = - Symmetric(_hermorsympart!(identity, identity, A, char_uplo(uplo)), uplo) + Symmetric(_hermorsympart!(symmetricpart, transpose, A, char_uplo(uplo)), uplo) """ hermitianpart(A, uplo=:U) @@ -892,6 +893,7 @@ as a [`Hermitian`](@ref) matrix. """ hermitianpart(A::AbstractMatrix{T}, uplo=:U) where {T} = hermitianpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) +hermitianpart(x::Number) = real(x) """ hermitianpart!(A, uplo=:U) @@ -903,7 +905,7 @@ and return `Hermitian(A, uplo)`. This function requires Julia 1.10 or later. """ hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = - Hermitian(_hermorsympart!(real, conj, A, char_uplo(uplo)), uplo) + Hermitian(_hermorsympart!(hermitianpart, adjoint, A, char_uplo(uplo)), uplo) function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix, uplo::Char) require_one_based_indexing(A) From 764913cc517f8e1a94b5a22044a921d8a55ffbe3 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 20 Jan 2023 10:35:42 +0100 Subject: [PATCH 05/18] Update stdlib/LinearAlgebra/src/symmetric.jl --- stdlib/LinearAlgebra/src/symmetric.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index e9bcd1f9760c5..5c720425f40c7 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -910,10 +910,9 @@ hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix, uplo::Char) require_one_based_indexing(A) n = checksquare(A) - triangle = uplo === 'U' ? (j -> 1:(j - 1)) : (j -> (j + 1):n) @inbounds for j in 1:n A[j, j] = real(A[j, j]) - for i in triangle(j) + for i in (uplo === 'U' ? (1:(j - 1)) : ((j + 1):n)) A[i, j] = (A[i, j] + conj(A[j, i])) / 2 end end From c32571aca3040e22539ab8b86ea584096b63fc53 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 20 Jan 2023 11:08:07 +0100 Subject: [PATCH 06/18] expand docstrings --- stdlib/LinearAlgebra/src/symmetric.jl | 94 ++++++++++++++++----------- 1 file changed, 55 insertions(+), 39 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 5c720425f40c7..afcbe3b6fce6c 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -17,34 +17,37 @@ end Construct a `Symmetric` view of the upper (if `uplo = :U`) or lower (if `uplo = :L`) triangle of the matrix `A`. +To compute the Hermitian part `(A + A') / 2` of `A`, use [`hermitianpart`](@ref). + # Examples ```jldoctest -julia> A = [1 0 2 0 3; 0 4 0 5 0; 6 0 7 0 8; 0 9 0 1 0; 2 0 3 0 4] -5×5 Matrix{Int64}: - 1 0 2 0 3 - 0 4 0 5 0 - 6 0 7 0 8 - 0 9 0 1 0 - 2 0 3 0 4 +julia> A = [1 2 3; 4 5 6; 7 8 9] +3×3 Matrix{Int64}: + 1 2 3 + 4 5 6 + 7 8 9 julia> Supper = Symmetric(A) -5×5 Symmetric{Int64, Matrix{Int64}}: - 1 0 2 0 3 - 0 4 0 5 0 - 2 0 7 0 8 - 0 5 0 1 0 - 3 0 8 0 4 +3×3 Symmetric{Int64, Matrix{Int64}}: + 1 2 3 + 2 5 6 + 3 6 9 julia> Slower = Symmetric(A, :L) -5×5 Symmetric{Int64, Matrix{Int64}}: - 1 0 6 0 2 - 0 4 0 9 0 - 6 0 7 0 3 - 0 9 0 1 0 - 2 0 3 0 4 +3×3 Symmetric{Int64, Matrix{Int64}}: + 1 4 7 + 4 5 8 + 7 8 9 + +julia> hermitianpart(A) +3×3 Hermitian{Float64, Matrix{Float64}}: + 1.0 3.0 5.0 + 3.0 5.0 7.0 + 5.0 7.0 9.0 ``` -Note that `Supper` will not be equal to `Slower` unless `A` is itself symmetric (e.g. if `A == transpose(A)`). +Note that `Supper` will not be equal to `Slower` unless `A` is itself symmetric (e.g. if +`A == transpose(A)`). """ function Symmetric(A::AbstractMatrix, uplo::Symbol=:U) checksquare(A) @@ -99,25 +102,33 @@ end Construct a `Hermitian` view of the upper (if `uplo = :U`) or lower (if `uplo = :L`) triangle of the matrix `A`. +To compute the Hermitian part of `A`, use [`hermitianpart`](@ref). + # Examples ```jldoctest -julia> A = [1 0 2+2im 0 3-3im; 0 4 0 5 0; 6-6im 0 7 0 8+8im; 0 9 0 1 0; 2+2im 0 3-3im 0 4]; +julia> A = [1 2+2im 3-3im; 4 5 6-6im; 7 8+8im 9] +3×3 Matrix{Complex{Int64}}: + 1+0im 2+2im 3-3im + 4+0im 5+0im 6-6im + 7+0im 8+8im 9+0im julia> Hupper = Hermitian(A) -5×5 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}: - 1+0im 0+0im 2+2im 0+0im 3-3im - 0+0im 4+0im 0+0im 5+0im 0+0im - 2-2im 0+0im 7+0im 0+0im 8+8im - 0+0im 5+0im 0+0im 1+0im 0+0im - 3+3im 0+0im 8-8im 0+0im 4+0im +3×3 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}: + 1+0im 2+2im 3-3im + 2-2im 5+0im 6-6im + 3+3im 6+6im 9+0im julia> Hlower = Hermitian(A, :L) -5×5 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}: - 1+0im 0+0im 6+6im 0+0im 2-2im - 0+0im 4+0im 0+0im 9+0im 0+0im - 6-6im 0+0im 7+0im 0+0im 3+3im - 0+0im 9+0im 0+0im 1+0im 0+0im - 2+2im 0+0im 3-3im 0+0im 4+0im +3×3 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}: + 1+0im 4+0im 7+0im + 4+0im 5+0im 8-8im + 7+0im 8+8im 9+0im + +julia> hermitianpart(A) +3×3 Hermitian{ComplexF64, Matrix{ComplexF64}}: + 1.0+0.0im 3.0+1.0im 5.0-1.5im + 3.0-1.0im 5.0+0.0im 7.0-7.0im + 5.0+1.5im 7.0+7.0im 9.0+0.0im ``` Note that `Hupper` will not be equal to `Hlower` unless `A` is itself Hermitian (e.g. if `A == adjoint(A)`). @@ -883,23 +894,28 @@ symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) = Symmetric(_hermorsympart!(symmetricpart, transpose, A, char_uplo(uplo)), uplo) """ - hermitianpart(A, uplo=:U) + hermitianpart(A, uplo=:U) -> Hermitian Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, -as a [`Hermitian`](@ref) matrix. +as a [`Hermitian`](@ref) matrix. For real matrices `A`, this is also known as the symmetric +part of `A`. The optional argument `uplo` controls whether the upper +(`uplo = :U`) or lower (`uplo = :L`) triangle of the underlying matrix is explicitly +filled. The opposite triangle is then only implicitly filled via the [`Hermitian`](@ref) +view, which, for real matrices, is equivalent to a [`Symmetric`](@ref) view. !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ hermitianpart(A::AbstractMatrix{T}, uplo=:U) where {T} = - hermitianpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) + hermitianpart!(copy_similar(A, typeof(one(T) / 2)), uplo) hermitianpart(x::Number) = real(x) """ - hermitianpart!(A, uplo=:U) + hermitianpart!(A, uplo=:U) -> Hermitian -Overwrite the square matrix `A` with its Hermitian part, `(A + A') / 2`, -and return `Hermitian(A, uplo)`. +Overwrite the upper (`uplo = :U`) or lower (`uplo = :L`) triangle of the square matrix `A` +with its Hermitian part `(A + A') / 2`, and return `Hermitian(A, uplo)`. For real matrices +`A`, this is also known as the symmetric part of `A`. !!! compat "Julia 1.10" This function requires Julia 1.10 or later. From 49618a4a733dae6a8587450928d2a00efd6f1431 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 20 Jan 2023 17:31:56 +0100 Subject: [PATCH 07/18] Update stdlib/LinearAlgebra/src/symmetric.jl Co-authored-by: Steven G. Johnson --- stdlib/LinearAlgebra/src/symmetric.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index afcbe3b6fce6c..8bc54426be1db 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -17,7 +17,15 @@ end Construct a `Symmetric` view of the upper (if `uplo = :U`) or lower (if `uplo = :L`) triangle of the matrix `A`. -To compute the Hermitian part `(A + A') / 2` of `A`, use [`hermitianpart`](@ref). +`Symmetric` views are mainly useful for real-symmetric matrices, for which +specialized algorithms (e.g. for eigenproblems) are enabled for `Symmetric` types. +More generally, see also `Hermitian(A)` for Hermitian matrices `A == A'`, which +is effectively equivalent to `Symmetric` for real matrices but is also useful for +complex matrices. (Whereas complex `Symmetric` matrices are supported but have few +if any specialized algorithms.) + +To compute the symmetric part of a real matrix, or more generally the Hermitian part `(A + A') / 2` of +a real or complex matrix `A`, use [`hermitianpart`](@ref). # Examples ```jldoctest From 46e3b943b4b7e19278eb1f7a5c48f29f299f1f22 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 23 Jan 2023 17:50:58 +0100 Subject: [PATCH 08/18] Update stdlib/LinearAlgebra/src/symmetric.jl Co-authored-by: Steven G. Johnson --- stdlib/LinearAlgebra/src/symmetric.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 8bc54426be1db..26e6efd187c14 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -19,7 +19,7 @@ triangle of the matrix `A`. `Symmetric` views are mainly useful for real-symmetric matrices, for which specialized algorithms (e.g. for eigenproblems) are enabled for `Symmetric` types. -More generally, see also `Hermitian(A)` for Hermitian matrices `A == A'`, which +More generally, see also [`Hermitian(A)`](@ref) for Hermitian matrices `A == A'`, which is effectively equivalent to `Symmetric` for real matrices but is also useful for complex matrices. (Whereas complex `Symmetric` matrices are supported but have few if any specialized algorithms.) From de9f4a86bf36ac23209e2617554a8d70e5ccf266 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 23 Jan 2023 17:55:19 +0100 Subject: [PATCH 09/18] more xrefs --- stdlib/LinearAlgebra/src/symmetric.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 26e6efd187c14..ee8af51bba5a3 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -911,6 +911,8 @@ part of `A`. The optional argument `uplo` controls whether the upper filled. The opposite triangle is then only implicitly filled via the [`Hermitian`](@ref) view, which, for real matrices, is equivalent to a [`Symmetric`](@ref) view. +See also [`hermitianpart!`](@ref) for the corresponding in-place operation. + !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ @@ -925,6 +927,8 @@ Overwrite the upper (`uplo = :U`) or lower (`uplo = :L`) triangle of the square with its Hermitian part `(A + A') / 2`, and return `Hermitian(A, uplo)`. For real matrices `A`, this is also known as the symmetric part of `A`. +See also [`hermitianpart`](@ref). + !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ From 64124b96ef038097493002a102fd92a5eca51db4 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 27 Jan 2023 09:37:16 +0100 Subject: [PATCH 10/18] allow full hermitianization --- stdlib/LinearAlgebra/src/symmetric.jl | 29 ++++++++++--------- stdlib/LinearAlgebra/test/symmetric.jl | 39 +++++++++++++++----------- 2 files changed, 38 insertions(+), 30 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index ee8af51bba5a3..f8eee53f7f7ff 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -902,26 +902,28 @@ symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) = Symmetric(_hermorsympart!(symmetricpart, transpose, A, char_uplo(uplo)), uplo) """ - hermitianpart(A, uplo=:U) -> Hermitian + hermitianpart(A, uplo=:U, full=false) -> Hermitian Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, as a [`Hermitian`](@ref) matrix. For real matrices `A`, this is also known as the symmetric -part of `A`. The optional argument `uplo` controls whether the upper -(`uplo = :U`) or lower (`uplo = :L`) triangle of the underlying matrix is explicitly -filled. The opposite triangle is then only implicitly filled via the [`Hermitian`](@ref) -view, which, for real matrices, is equivalent to a [`Symmetric`](@ref) view. +part of `A`. The first optional argument `uplo` controls whether the upper +(`uplo = :U`) or lower (`uplo = :L`) triangle of the matrix underlying the [`Hermitian`](@ref) +view is explicitly filled. Unless the second optional argument `full=true` (in which case the +full underlying matrix is filled with the Hermitian part), the opposite triangle is then +only implicitly filled via the [`Hermitian`](@ref) view. For real matrices, the latter is +equivalent to a [`Symmetric`](@ref) view. See also [`hermitianpart!`](@ref) for the corresponding in-place operation. !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ -hermitianpart(A::AbstractMatrix{T}, uplo=:U) where {T} = - hermitianpart!(copy_similar(A, typeof(one(T) / 2)), uplo) -hermitianpart(x::Number) = real(x) +hermitianpart(A::AbstractMatrix{T}, uplo=:U, full::Bool = false) where {T} = + hermitianpart!(copy_similar(A, typeof(one(T) / 2)), uplo, full) +hermitianpart(x::Number, uplo=:U, full::Bool=false) = real(x) """ - hermitianpart!(A, uplo=:U) -> Hermitian + hermitianpart!(A, uplo=:U, full=false) -> Hermitian Overwrite the upper (`uplo = :U`) or lower (`uplo = :L`) triangle of the square matrix `A` with its Hermitian part `(A + A') / 2`, and return `Hermitian(A, uplo)`. For real matrices @@ -932,16 +934,17 @@ See also [`hermitianpart`](@ref). !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ -hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = - Hermitian(_hermorsympart!(hermitianpart, adjoint, A, char_uplo(uplo)), uplo) +hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U, full::Bool = false) = + Hermitian(_hermorsympart!(a -> hermitianpart(a, uplo, full), adjoint, A, char_uplo(uplo), full), uplo) -function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix, uplo::Char) +@inline function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix, uplo::Char, full::Bool) require_one_based_indexing(A) n = checksquare(A) @inbounds for j in 1:n A[j, j] = real(A[j, j]) for i in (uplo === 'U' ? (1:(j - 1)) : ((j + 1):n)) - A[i, j] = (A[i, j] + conj(A[j, i])) / 2 + A[i, j] = val = (A[i, j] + conj(A[j, i])) / 2 + full && (A[j, i] = conj(val)) end end return A diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 26727e48a6299..dbf95252435b8 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -779,26 +779,31 @@ end @testset "symmetric and hermitian parts" begin for T in [Float32, Complex{Float32}, Int32, Rational{Int32}, Complex{Int32}, Complex{Rational{Int32}}] - for (f, f!, t) in [(hermitianpart, hermitianpart!, adjoint), - (symmetricpart, symmetricpart!, transpose)] - X = T[1 2 3; 4 5 6; 7 8 9] - T <: Complex && (X .+= im .* X) + f, f!, t = hermitianpart, hermitianpart!, T <: Real ? transpose : adjoint + X = T[1 2 3; 4 5 6; 7 8 9] + T <: Complex && (X .+= im .* X) + Xc = copy(X) + Y = (X + t(X)) / 2 + U = f(X) + L = f(X, :L) + @test U == L == Y + U = f(X, :U, true) + L = f(X, :L, true) + @test parent(U) == U == parent(L) == L == Y + if T <: AbstractFloat || real(T) <: AbstractFloat + HU = f!(copy(X)) + @test HU == Y + HU = f!(Xc, :U, true) + @test parent(HU) == HU == Y + @test triu(Xc) == triu(Y) + HL = f!(copy(X), :L) + @test HL == Y Xc = copy(X) - Y = (X + t(X)) / 2 - U = f(X) - L = f(X, :L) - @test U == L == Y - if T <: AbstractFloat || real(T) <: AbstractFloat - HU = f!(X) - @test HU == Y - @test triu(X) == triu(Y) - HL = f!(Xc, :L) - @test HL == Y - @test tril(Xc) == tril(Y) - end + HL = f!(Xc, :L, true) + @test parent(HL) == HL == Y + @test tril(Xc) == tril(Y) end end - @test_throws DimensionMismatch symmetricpart([1.0 2.0 3.0; 4.0 5.0 6.0]) @test_throws DimensionMismatch hermitianpart([1.0 2.0 3.0; 4.0 5.0 6.0]) end From 62bc5638884f7abcce2d1edb44039b1809e09867 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 30 Jan 2023 12:16:22 +0100 Subject: [PATCH 11/18] rm symmetricpart, write always full, add matofmat tests --- NEWS.md | 4 +- stdlib/LinearAlgebra/docs/src/index.md | 2 - stdlib/LinearAlgebra/src/symmetric.jl | 64 ++++++++------------------ stdlib/LinearAlgebra/test/symmetric.jl | 20 ++++---- 4 files changed, 29 insertions(+), 61 deletions(-) diff --git a/NEWS.md b/NEWS.md index 2f47267397070..6e692814227f1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -43,9 +43,11 @@ Standard library changes packages being loaded in the Julia session. This has similar applications as the Requires.jl package but also supports precompilation and setting compatibility. + #### LinearAlgebra -* New functions `symmetricpart`, `symmetricpart!`, `hermitianpart`, and `hermitianpart!` for extracting the symmetric or Hermitian part of a matrix ([#31836]). +- New functions `hermitianpart`, and `hermitianpart!` for extracting the Hermitian + (real symmetric) part of a matrix ([#31836]). #### Printf diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index b7dff5072783d..8eb5af5605b91 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -467,8 +467,6 @@ Base.copy(::Union{Transpose,Adjoint}) LinearAlgebra.stride1 LinearAlgebra.checksquare LinearAlgebra.peakflops -LinearAlgebra.symmetricpart -LinearAlgebra.symmetricpart! LinearAlgebra.hermitianpart LinearAlgebra.hermitianpart! ``` diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index f8eee53f7f7ff..58e93b77b4a34 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -877,74 +877,46 @@ for func in (:log, :sqrt) end """ - symmetricpart(A, uplo=:U) + hermitianpart(A, uplo=:U) -> Hermitian -Return the symmetric part of the square matrix `A`, defined as `(A + transpose(A)) / 2`, -as a [`Symmetric`](@ref) matrix. - -!!! compat "Julia 1.10" - This function requires Julia 1.10 or later. -""" -symmetricpart(A::AbstractMatrix{T}, uplo=:U) where {T} = - symmetricpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo) -symmetricpart(x::Number) = x - -""" - symmetricpart!(A, uplo=:U) - -Overwrite the square matrix `A` with its symmetric part, `(A + transpose(A)) / 2`, -and return `Symmetric(A, uplo)`. - -!!! compat "Julia 1.10" - This function requires Julia 1.10 or later. -""" -symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) = - Symmetric(_hermorsympart!(symmetricpart, transpose, A, char_uplo(uplo)), uplo) - -""" - hermitianpart(A, uplo=:U, full=false) -> Hermitian - -Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, -as a [`Hermitian`](@ref) matrix. For real matrices `A`, this is also known as the symmetric -part of `A`. The first optional argument `uplo` controls whether the upper -(`uplo = :U`) or lower (`uplo = :L`) triangle of the matrix underlying the [`Hermitian`](@ref) -view is explicitly filled. Unless the second optional argument `full=true` (in which case the -full underlying matrix is filled with the Hermitian part), the opposite triangle is then -only implicitly filled via the [`Hermitian`](@ref) view. For real matrices, the latter is -equivalent to a [`Symmetric`](@ref) view. +Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, as a +[`Hermitian`](@ref) matrix. For real matrices `A`, this is also known as the symmetric part +of `A`. The optional argument `uplo` controls the corresponding argument of the +[`Hermitian`](@ref) view. For real matrices, the latter is equivalent to a +[`Symmetric`](@ref) view. See also [`hermitianpart!`](@ref) for the corresponding in-place operation. !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ -hermitianpart(A::AbstractMatrix{T}, uplo=:U, full::Bool = false) where {T} = - hermitianpart!(copy_similar(A, typeof(one(T) / 2)), uplo, full) -hermitianpart(x::Number, uplo=:U, full::Bool=false) = real(x) +hermitianpart(A::AbstractMatrix{T}, uplo::Symbol=:U) where {T} = + hermitianpart!(copy_similar(A, Base.promote_op(/, T, Int)), uplo) +hermitianpart(x::Number, uplo::Symbol=:U) = real(x) """ - hermitianpart!(A, uplo=:U, full=false) -> Hermitian + hermitianpart!(A, uplo=:U) -> Hermitian -Overwrite the upper (`uplo = :U`) or lower (`uplo = :L`) triangle of the square matrix `A` -with its Hermitian part `(A + A') / 2`, and return `Hermitian(A, uplo)`. For real matrices -`A`, this is also known as the symmetric part of `A`. +Overwrite the square matrix `A` with its Hermitian part `(A + A') / 2`, and return +[`Hermitian(A, uplo)`](@ref). For real matrices `A`, this is also known as the symmetric +part of `A`. See also [`hermitianpart`](@ref). !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ -hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U, full::Bool = false) = - Hermitian(_hermorsympart!(a -> hermitianpart(a, uplo, full), adjoint, A, char_uplo(uplo), full), uplo) +hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = + Hermitian(_hermorsympart!(a -> hermitianpart(a, uplo), adjoint, A), uplo) -@inline function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix, uplo::Char, full::Bool) +@inline function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix) require_one_based_indexing(A) n = checksquare(A) @inbounds for j in 1:n A[j, j] = real(A[j, j]) - for i in (uplo === 'U' ? (1:(j - 1)) : ((j + 1):n)) + for i in 1:j-1 A[i, j] = val = (A[i, j] + conj(A[j, i])) / 2 - full && (A[j, i] = conj(val)) + A[j, i] = conj(val) end end return A diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index dbf95252435b8..9d32eaad11e77 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -787,24 +787,20 @@ end U = f(X) L = f(X, :L) @test U == L == Y - U = f(X, :U, true) - L = f(X, :L, true) - @test parent(U) == U == parent(L) == L == Y if T <: AbstractFloat || real(T) <: AbstractFloat - HU = f!(copy(X)) + HU = f!(X) @test HU == Y - HU = f!(Xc, :U, true) - @test parent(HU) == HU == Y - @test triu(Xc) == triu(Y) - HL = f!(copy(X), :L) + @test triu(X) == triu(Y) + HL = f!(Xc, :L) @test HL == Y - Xc = copy(X) - HL = f!(Xc, :L, true) - @test parent(HL) == HL == Y @test tril(Xc) == tril(Y) end end - @test_throws DimensionMismatch hermitianpart([1.0 2.0 3.0; 4.0 5.0 6.0]) + @test_throws DimensionMismatch hermitianpart(ones(1,2)) + for T in (Float64, ComplexF64) + A = [randn(T, 2, 2) for _ in 1:2, _ in 1:2] + @test hermitianpart(A) == (A + A')/2 + end end end # module TestSymmetric From 62235bf5f6e2de43f600737021ae2a3f2f5076b5 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 31 Jan 2023 18:04:05 +0100 Subject: [PATCH 12/18] more tests --- stdlib/LinearAlgebra/test/symmetric.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 9d32eaad11e77..90ea9c368b958 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -776,7 +776,7 @@ end end end -@testset "symmetric and hermitian parts" begin +@testset "symmetric and hermitian part" begin for T in [Float32, Complex{Float32}, Int32, Rational{Int32}, Complex{Int32}, Complex{Rational{Int32}}] f, f!, t = hermitianpart, hermitianpart!, T <: Real ? transpose : adjoint @@ -786,6 +786,10 @@ end Y = (X + t(X)) / 2 U = f(X) L = f(X, :L) + @test U isa Hermitian + @test L isa Hermitian + @test U.uplo == 'U' + @test L.uplo == 'L' @test U == L == Y if T <: AbstractFloat || real(T) <: AbstractFloat HU = f!(X) @@ -797,9 +801,12 @@ end end end @test_throws DimensionMismatch hermitianpart(ones(1,2)) - for T in (Float64, ComplexF64) + for T in (Float64, ComplexF64), uplo in (:U, :L) A = [randn(T, 2, 2) for _ in 1:2, _ in 1:2] - @test hermitianpart(A) == (A + A')/2 + Aherm = hermitianpart(A, uplo) + @test Aherm == Aherm.data == (A + A')/2 + @test Aherm isa Hermitian + @test Aherm.uplo == LinearAlgebra.char_uplo(uplo) end end From 97838cad2e4f0333117e8b7b43371d53f9e6e3f1 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 1 Feb 2023 09:17:58 +0100 Subject: [PATCH 13/18] Update stdlib/LinearAlgebra/src/LinearAlgebra.jl Co-authored-by: Alex Arslan --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index ea8ee00be8aaf..96c2e38d187e8 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -146,8 +146,6 @@ export svdvals!, svdvals, sylvester, - symmetricpart, - symmetricpart!, tr, transpose, transpose!, From e40148b5c0d034f13535ad0868b5b17d720e9ef6 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 1 Feb 2023 09:32:58 +0100 Subject: [PATCH 14/18] rearrange code Co-authored-by: Steven G. Johnson --- stdlib/LinearAlgebra/src/symmetric.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 58e93b77b4a34..ce5de0138ae72 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -890,9 +890,8 @@ See also [`hermitianpart!`](@ref) for the corresponding in-place operation. !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ -hermitianpart(A::AbstractMatrix{T}, uplo::Symbol=:U) where {T} = - hermitianpart!(copy_similar(A, Base.promote_op(/, T, Int)), uplo) -hermitianpart(x::Number, uplo::Symbol=:U) = real(x) +hermitianpart(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart(A), uplo) +hermitianpart(a::Number) = _hermitianpart(a) """ hermitianpart!(A, uplo=:U) -> Hermitian @@ -906,17 +905,19 @@ See also [`hermitianpart`](@ref). !!! compat "Julia 1.10" This function requires Julia 1.10 or later. """ -hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = - Hermitian(_hermorsympart!(a -> hermitianpart(a, uplo), adjoint, A), uplo) +hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart!(A), uplo) -@inline function _hermorsympart!(real::Function, conj::Function, A::AbstractMatrix) +_hermitianpart(A::AbstractMatrix) = _hermitianpart!(copy_similar(A, Base.promote_op(/, eltype(A), Int))) +_hermitianpart(a::Number) = real(a) + +function _hermitianpart!(A) require_one_based_indexing(A) n = checksquare(A) @inbounds for j in 1:n - A[j, j] = real(A[j, j]) + A[j, j] = _hermitianpart(A[j, j]) for i in 1:j-1 - A[i, j] = val = (A[i, j] + conj(A[j, i])) / 2 - A[j, i] = conj(val) + A[i, j] = val = (A[i, j] + adjoint(A[j, i])) / 2 + A[j, i] = adjoint(val) end end return A From 94e5ceea9f3f2ba7b84a3473da50d15473e81db0 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 1 Feb 2023 10:12:06 +0100 Subject: [PATCH 15/18] make returned eltype consistent --- stdlib/LinearAlgebra/src/symmetric.jl | 2 +- stdlib/LinearAlgebra/test/symmetric.jl | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index ce5de0138ae72..528d7c24ef09b 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -908,7 +908,7 @@ See also [`hermitianpart`](@ref). hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart!(A), uplo) _hermitianpart(A::AbstractMatrix) = _hermitianpart!(copy_similar(A, Base.promote_op(/, eltype(A), Int))) -_hermitianpart(a::Number) = real(a) +_hermitianpart(a::T) where {T<:Number} = convert(Base.promote_op(/, T, Int), real(a)) function _hermitianpart!(A) require_one_based_indexing(A) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 90ea9c368b958..08e171edae4b2 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -808,6 +808,8 @@ end @test Aherm isa Hermitian @test Aherm.uplo == LinearAlgebra.char_uplo(uplo) end + z = 3 + im + hermitianpart(z) === only(hermitianpart([z;;])) end end # module TestSymmetric From 064f74b70a869e05531158bbfedb848b6edcb9f4 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 1 Feb 2023 18:10:21 +0100 Subject: [PATCH 16/18] Apply suggestions from code review Co-authored-by: Steven G. Johnson --- stdlib/LinearAlgebra/src/symmetric.jl | 8 ++++---- stdlib/LinearAlgebra/test/symmetric.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 528d7c24ef09b..28c1a694cd330 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -881,7 +881,7 @@ end Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, as a [`Hermitian`](@ref) matrix. For real matrices `A`, this is also known as the symmetric part -of `A`. The optional argument `uplo` controls the corresponding argument of the +of `A`; it is also sometimes called the "operator real part". The optional argument `uplo` controls the corresponding argument of the [`Hermitian`](@ref) view. For real matrices, the latter is equivalent to a [`Symmetric`](@ref) view. @@ -896,11 +896,11 @@ hermitianpart(a::Number) = _hermitianpart(a) """ hermitianpart!(A, uplo=:U) -> Hermitian -Overwrite the square matrix `A` with its Hermitian part `(A + A') / 2`, and return +Overwrite the square matrix `A` in-place with its Hermitian part `(A + A') / 2`, and return [`Hermitian(A, uplo)`](@ref). For real matrices `A`, this is also known as the symmetric part of `A`. -See also [`hermitianpart`](@ref). +See also [`hermitianpart`](@ref) for the corresponding out-of-place operation. !!! compat "Julia 1.10" This function requires Julia 1.10 or later. @@ -910,7 +910,7 @@ hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart!(A _hermitianpart(A::AbstractMatrix) = _hermitianpart!(copy_similar(A, Base.promote_op(/, eltype(A), Int))) _hermitianpart(a::T) where {T<:Number} = convert(Base.promote_op(/, T, Int), real(a)) -function _hermitianpart!(A) +function _hermitianpart!(A::AbstractMatrix) require_one_based_indexing(A) n = checksquare(A) @inbounds for j in 1:n diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 08e171edae4b2..bb5561860e380 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -776,7 +776,7 @@ end end end -@testset "symmetric and hermitian part" begin +@testset "hermitian part" begin for T in [Float32, Complex{Float32}, Int32, Rational{Int32}, Complex{Int32}, Complex{Rational{Int32}}] f, f!, t = hermitianpart, hermitianpart!, T <: Real ? transpose : adjoint From 96d9061b0845afa6d735f6084d1f0b5df0382cc7 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 2 Feb 2023 09:09:04 +0100 Subject: [PATCH 17/18] rm hermitianpart(::Number) --- stdlib/LinearAlgebra/src/symmetric.jl | 3 +-- stdlib/LinearAlgebra/test/symmetric.jl | 2 -- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 28c1a694cd330..fdbfe59302a87 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -891,7 +891,6 @@ See also [`hermitianpart!`](@ref) for the corresponding in-place operation. This function requires Julia 1.10 or later. """ hermitianpart(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart(A), uplo) -hermitianpart(a::Number) = _hermitianpart(a) """ hermitianpart!(A, uplo=:U) -> Hermitian @@ -908,7 +907,7 @@ See also [`hermitianpart`](@ref) for the corresponding out-of-place operation. hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart!(A), uplo) _hermitianpart(A::AbstractMatrix) = _hermitianpart!(copy_similar(A, Base.promote_op(/, eltype(A), Int))) -_hermitianpart(a::T) where {T<:Number} = convert(Base.promote_op(/, T, Int), real(a)) +_hermitianpart(a::Number) = real(a) function _hermitianpart!(A::AbstractMatrix) require_one_based_indexing(A) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index bb5561860e380..284264d726202 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -808,8 +808,6 @@ end @test Aherm isa Hermitian @test Aherm.uplo == LinearAlgebra.char_uplo(uplo) end - z = 3 + im - hermitianpart(z) === only(hermitianpart([z;;])) end end # module TestSymmetric From 08e50dd22db1845da9751cda1c27dfa7a90da56d Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 14 Feb 2023 08:00:39 -0500 Subject: [PATCH 18/18] Update NEWS.md --- NEWS.md | 1 - 1 file changed, 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index d875828502f36..ede34d3238ebe 100644 --- a/NEWS.md +++ b/NEWS.md @@ -46,7 +46,6 @@ Standard library changes packages being loaded in the Julia session. This has similar applications as the Requires.jl package but also supports precompilation and setting compatibility. - * `Pkg.precompile` now accepts `timing` as a keyword argument which displays per package timing information for precompilation (e.g. `Pkg.precompile(timing=true)`) #### LinearAlgebra