Repository: JuliaLang/julia
File: stdlib/LinearAlgebra/src/symmetric.jl (Julia)

Code before:
        function ($func)(A::Hermitian{T}; $(rtolarg...)) where {T<:Complex}
            n = checksquare(A)
            F = eigen(A)
            λ₀ = $rtolval # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff
            if all(λ -> λ ≥ λ₀, F.values)
                retmat = (F.vectors * Diagonal(($func).(max.(0, F.values)))) * F.vectors'
                for i = 1:n
                    retmat[i,i] = real(retmat[i,i])
                end
                return Hermitian(retmat)
            else
                retmat = (F.vectors * Diagonal(($func).(complex(F.values)))) * F.vectors'
                return retmat
            end
        end
    end
end

"""
    symmetricpart(A, uplo=:U)

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.
"""
symmetricpart(A::AbstractMatrix{T}, uplo=:U) where {T} =
    symmetricpart!(copyto!(similar(A, typeof(one(T) / 2)), A), uplo)

"""
    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.
"""
symmetricpart!(A::AbstractMatrix, uplo::Symbol=:U) =
    Symmetric(_hermorsympart!(identity, identity, A, char_uplo(uplo)), uplo)

"""
    hermitianpart(A, uplo=: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.
"""

Review comment (suggestion):
```suggestion
!!! compat "Julia 1.10"
    This function requires Julia 1.10 or later.
```

Code after:
    for i in 1:n
        S.diag[i], C.diag[i] = sincos(F.values[i])
    end
    retmatS, retmatC = (F.vectors * S) * F.vectors', (F.vectors * C) * F.vectors'
    for i = 1:n
        retmatS[i,i] = real(retmatS[i,i])
        retmatC[i,i] = real(retmatC[i,i])
    end
    return Hermitian(retmatS), Hermitian(retmatC)
end


for func in (:log, :sqrt)
    # sqrt has rtol arg to handle matrices that are semidefinite up to roundoff errors
    rtolarg = func === :sqrt ? Any[Expr(:kw, :(rtol::Real), :(eps(real(float(one(T))))*size(A,1)))] : Any[]
    rtolval = func === :sqrt ? :(-maximum(abs, F.values) * rtol) : 0
    @eval begin
        function ($func)(A::HermOrSym{T}; $(rtolarg...)) where {T<:Real}
            F = eigen(A)
            λ₀ = $rtolval # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff
            if all(λ -> λ ≥ λ₀, F.values)
                retmat = (F.vectors * Diagonal(($func).(max.(0, F.values)))) * F.vectors'
            else
                retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
            end
            return Symmetric(retmat)
        end

        function ($func)(A::Hermitian{T}; $(rtolarg...)) where {T<:Complex}
            n = checksquare(A)
            F = eigen(A)
            λ₀ = $rtolval # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff
            if all(λ -> λ ≥ λ₀, F.values)
                retmat = (F.vectors * Diagonal(($func).(max.(0, F.values)))) * F.vectors'
                for i = 1:n
                    retmat[i,i] = real(retmat[i,i])
                end
                return Hermitian(retmat)
            else
                retmat = (F.vectors * Diagonal(($func).(complex(F.values)))) * F.vectors'
                return retmat
            end
        end
    end
end

"""
    hermitianpart(A, uplo=:U) -> 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