Skip to content

Commit 7a645dd

Browse files
jishnubdkarrasch
andauthored
Avoid materializing arrays in bidiag matmul (#55450)
Currently, small `Bidiagonal`/`Tridiagonal` matrices are materialized in matrix multiplications, but this is wasteful and unnecessary. This PR changes this to use a naive matrix multiplication for small matrices, and fall back to the banded multiplication for larger ones. Multiplication by a `Bidiagonal` falls back to a banded matrix multiplication for all sizes in the current implementation, and iterates in a cache-friendly manner for the non-`Bidiagonal` matrix. In certain cases, the matrices were being materialized if the non-structured matrix was small, even if the structured matrix was large. This is changed as well in this PR. Some improvements in performance: ```julia julia> B = Bidiagonal(rand(3), rand(2), :U); A = rand(size(B)...); C = similar(A); julia> @Btime mul!($C, $A, $B); 193.152 ns (6 allocations: 352 bytes) # nightly v"1.12.0-DEV.1034" 18.826 ns (0 allocations: 0 bytes) # This PR julia> T = Tridiagonal(rand(99), rand(100), rand(99)); A = rand(2, size(T,2)); C = similar(A); julia> @Btime mul!($C, $A, $T); 9.398 μs (8 allocations: 79.94 KiB) # nightly 416.407 ns (0 allocations: 0 bytes) # This PR julia> B = Bidiagonal(rand(300), rand(299), :U); A = rand(20000, size(B,2)); C = similar(A); julia> @Btime mul!($C, $A, $B); 33.395 ms (0 allocations: 0 bytes) # nightly 6.695 ms (0 allocations: 0 bytes) # This PR (cache-friendly) ``` Closes #55414 --------- Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>
1 parent 99b8868 commit 7a645dd

File tree

4 files changed

+422
-68
lines changed

4 files changed

+422
-68
lines changed

stdlib/LinearAlgebra/src/LinearAlgebra.jl

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -673,7 +673,9 @@ matprod_dest(A::Diagonal, B::Diagonal, TS) = _matprod_dest_diag(B, TS)
673673
_matprod_dest_diag(A, TS) = similar(A, TS)
674674
function _matprod_dest_diag(A::SymTridiagonal, TS)
675675
n = size(A, 1)
676-
Tridiagonal(similar(A, TS, n-1), similar(A, TS, n), similar(A, TS, n-1))
676+
ev = similar(A, TS, max(0, n-1))
677+
dv = similar(A, TS, n)
678+
Tridiagonal(ev, dv, similar(ev))
677679
end
678680

679681
# Special handling for adj/trans vec

0 commit comments

Comments
 (0)