Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove a few redundant diagonal ldiv! methods #38484

Merged
merged 3 commits into from
Apr 6, 2021

Conversation

mcognetta
Copy link
Contributor

Three ldiv! methods for diagonal matrices were effectively the exact same but had different type specializations. This PR replaces them with a single ldiv! that subsumes their functionality.

@dkarrasch dkarrasch added the linear algebra Linear algebra label Nov 18, 2020
Copy link
Member

@dkarrasch dkarrasch left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Two comments (summarizing what @KlausC said), and then this is ready to go.

stdlib/LinearAlgebra/src/diagonal.jl Outdated Show resolved Hide resolved
stdlib/LinearAlgebra/src/diagonal.jl Outdated Show resolved Hide resolved
@dkarrasch
Copy link
Member

Hm, actually, this really requires benchmarking. The strided version is optimized towards good access pattern for B, whereas the abstract version is optimized towards less reading of the diagonal elements of D. I think before changing anything here, we need to see some benchmarks.

@mcognetta
Copy link
Contributor Author

@dkarrasch I used the changes you and @KlausC suggested. Benchmarking revealed only minor performance differences:

Julia 1.5.0

julia> D = Diagonal(rand(500)); M = rand(500, 500); @btime ldiv!(D, M)
  580.819 μs (0 allocations: 0 bytes)

julia> D = Diagonal(rand(500)); M = UpperTriangular(rand(500, 500)); @btime ldiv!(D, M)
  533.071 μs (0 allocations: 0 bytes)

This PR:

julia> D = Diagonal(rand(500)); M = rand(500, 500); @btime ldiv!(D, M)
  536.923 μs (0 allocations: 0 bytes)

julia> D = Diagonal(rand(500)); M = UpperTriangular(rand(500, 500)); @btime ldiv!(D, M)
  521.033 μs (0 allocations: 0 bytes)

Copy link
Member

@dkarrasch dkarrasch left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The benchmark suggests that it makes sense to switch to the previous strided matrix loop order also for abstract matrices. This LGTM. The only question I have is whether we should try to preserve structure when the matrix has structure.

#Linear solver
function ldiv!(D::Diagonal, B::StridedVecOrMat)
(\)(D::Diagonal, A::AbstractMatrix) =
ldiv!(D, (typeof(oneunit(eltype(D))/oneunit(eltype(A)))).(A))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we change this to something like copy_oftype or copyto!(similar(A, typeof(...), size(A)), A) or so, in order to preserve the structure whenever possible?

@vtjnash vtjnash dismissed dkarrasch’s stale review April 3, 2021 18:04

I'm assuming the LGTM meant to dismiss this

@vtjnash vtjnash added the merge me PR is reviewed. Merge when all tests are passing label Apr 3, 2021
@vtjnash vtjnash merged commit 7352320 into JuliaLang:master Apr 6, 2021
@simeonschaub simeonschaub removed the merge me PR is reviewed. Merge when all tests are passing label Apr 28, 2021
ElOceanografo pushed a commit to ElOceanografo/julia that referenced this pull request May 4, 2021
antoine-levitt pushed a commit to antoine-levitt/julia that referenced this pull request May 9, 2021
johanmon pushed a commit to johanmon/julia that referenced this pull request Jul 5, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants