Skip to content

Conversation

@jishnub
Copy link
Member

@jishnub jishnub commented Dec 4, 2023

The current implementation assumes that the adjoint matrix may be copied to the destination, but this is not necessary in general. This PR changes the implementation to use mul! instead of copy + (l/r)mul! removes the specialized methods, so that the multiplication generalizes to arbitrary element types. In particular, the following works after this:

julia> SM = SMatrix{2,3}(reshape([1:6;],2,3))
2×3 SMatrix{2, 3, Int64, 6} with indices SOneTo(2)×SOneTo(3):
 1  3  5
 2  4  6

julia> D = Diagonal(fill(SM,2))
2×2 Diagonal{SMatrix{2, 3, Int64, 6}, Vector{SMatrix{2, 3, Int64, 6}}}:
 [1 3 5; 2 4 6]               
                [1 3 5; 2 4 6]

julia> A = fill(SM,2,2)'
2×2 adjoint(::Matrix{SMatrix{2, 3, Int64, 6}}) with eltype SMatrix{3, 2, Int64, 6}:
 [1 2; 3 4; 5 6]  [1 2; 3 4; 5 6]
 [1 2; 3 4; 5 6]  [1 2; 3 4; 5 6]

julia> D * A
2×2 Matrix{SMatrix{2, 2, Int64, 4}}:
 [35 44; 44 56]  [35 44; 44 56]
 [35 44; 44 56]  [35 44; 44 56]

julia> A * D
2×2 Matrix{SMatrix{3, 3, Int64, 9}}:
 [5 11 17; 11 25 39; 17 39 61]  [5 11 17; 11 25 39; 17 39 61]
 [5 11 17; 11 25 39; 17 39 61]  [5 11 17; 11 25 39; 17 39 61]

@jishnub jishnub added the linear algebra Linear algebra label Dec 4, 2023
@dkarrasch dkarrasch merged commit c1ca0d3 into master Dec 5, 2023
@dkarrasch dkarrasch deleted the jishnub/diagadjmul branch December 5, 2023 10:59
jishnub added a commit to JuliaLang/LinearAlgebra.jl that referenced this pull request Feb 20, 2025
This reinstates slightly altered versions of the methods that were
removed in JuliaLang/julia#52389. Sort of fixes
#1205, although this doesn't recover the full performance. However, this
version is more general, and works with the example presented in
JuliaLang/julia#52389. There's still a
performance regression, but the full performance may only be obtained
for mutable matrices, and we may not assume mutability in general.

Performance:
v1.10:

```julia
julia> n = 100
100

julia> A = adjoint(sparse(Float64, I, n, n));

julia> B = Diagonal(ones(n));

julia> @Btime $A * $B;
  837.119 ns (5 allocations: 2.59 KiB)
```

This PR
```julia
julia> @Btime $A * $B;
  1.106 μs (15 allocations: 5.56 KiB)
```
We need double the allocations here compared to earlier, as we firstly
materialize `D' * A'`, and then we again copy the adjoint of this
result. I wonder if this may be reduced.

---------

Co-authored-by: Daniel Karrasch <[email protected]>
dkarrasch added a commit to JuliaLang/LinearAlgebra.jl that referenced this pull request Feb 20, 2025
This reinstates slightly altered versions of the methods that were
removed in JuliaLang/julia#52389. Sort of fixes
#1205, although this doesn't recover the full performance. However, this
version is more general, and works with the example presented in
JuliaLang/julia#52389. There's still a
performance regression, but the full performance may only be obtained
for mutable matrices, and we may not assume mutability in general.

Performance:
v1.10:

```julia
julia> n = 100
100

julia> A = adjoint(sparse(Float64, I, n, n));

julia> B = Diagonal(ones(n));

julia> @Btime $A * $B;
  837.119 ns (5 allocations: 2.59 KiB)
```

This PR
```julia
julia> @Btime $A * $B;
  1.106 μs (15 allocations: 5.56 KiB)
```
We need double the allocations here compared to earlier, as we firstly
materialize `D' * A'`, and then we again copy the adjoint of this
result. I wonder if this may be reduced.

---------

Co-authored-by: Daniel Karrasch <[email protected]>
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.

2 participants