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

Add 3-arg * methods #37898

Merged
merged 32 commits into from
Jun 7, 2021
Merged
Changes from 1 commit
Commits
Show all changes
32 commits
Select commit Hold shift + click to select a range
2ed51d6
add 3-arg star
Oct 6, 2020
14fa53e
also handle 3 matrices
Oct 6, 2020
b0c33f4
a few more cases, and tests
Oct 6, 2020
1981525
two scalars + one array
Oct 6, 2020
71dfbb0
avoid dot dispatch
Oct 6, 2020
be5eee3
add some array-of-array tests
Oct 6, 2020
8e9bcd5
fix order of multiplication
Oct 6, 2020
cfa28cb
add a docstring
Oct 6, 2020
93d03b8
preserve order in fallbacks
Oct 6, 2020
b04b3a5
opt-in, take 1
Oct 6, 2020
1a29b74
remove two-scalar methods
Oct 6, 2020
6c3c9c0
fewer adjoint = less chance of spurious dimensions
Oct 6, 2020
a1ddb42
Update stdlib/LinearAlgebra/src/matmul.jl
mcabbott Oct 7, 2020
ea98220
four-argument *, why not
Oct 8, 2020
d6954e0
one more case
Oct 8, 2020
ffbb4b8
rm ambiguity
Oct 8, 2020
f0b3648
use 3-arg dot only when length(x)<64
Oct 17, 2020
1c92294
more explicit name for _SafeMatrix
Oct 17, 2020
0a71d33
don't use 3-arg dot at all
Oct 18, 2020
fb3b7da
use RealOrComplex for eltypes, too
Oct 27, 2020
f005fb5
better fallback for mat_vec_scalar
Nov 17, 2020
eac95e7
remove a now-duplicate definition of StridedMaybeAdjOrTransMat
Nov 20, 2020
2f39e33
Apply suggestions from code review
mcabbott Dec 11, 2020
d9b456e
constrain adjoint eltypes
Dec 11, 2020
766368f
two tiny fixes
mcabbott May 14, 2021
15f9f21
optimise dot(x,A,y)-like cases based on whether A is transposed
mcabbott May 16, 2021
5ef2922
two tests
mcabbott May 28, 2021
f778a00
Update stdlib/LinearAlgebra/test/matmul.jl
mcabbott May 28, 2021
eacc11b
change zero(γ) -> false as suggested
mcabbott Jun 7, 2021
ebb52fb
better docstring
mcabbott Jun 7, 2021
e7aa4ad
wording, plus examples
mcabbott Jun 7, 2021
468d19b
docstring
mcabbott Jun 7, 2021
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
wording, plus examples
  • Loading branch information
mcabbott committed Jun 7, 2021
commit e7aa4ad38625cb60ac1fedf30284f9893eb6a28d
33 changes: 25 additions & 8 deletions stdlib/LinearAlgebra/src/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1089,22 +1089,39 @@ const RealOrComplex = Union{Real,Complex}
*(A, B::AbstractMatrix, C)
A * B * C * D

Chained multiplication of 3 or 4 matrices is done in the most efficient order,
i.e. the order which minimises the total number of scalar multiplications involved,
based on the sizes of the arrays.

If the last factor is a vector, or the first a transposed vector,
it is efficient to do these first. In particular `x' * B * y` means `(x' * B) * y`
for an ordinary colum-major `B`. This is usually faster than `dot(x, B, y)`,
Chained multiplication of 3 or 4 matrices is done in the most efficient sequence,
based on the sizes of the arrays. That is, the number of scalar multiplications needed
for `(A * B) * C` (with 3 dense matrices) is compared to that for `A * (B * C)`
to choose which of these to execute.

If the last factor is a vector, or the first a transposed vector, then it is efficient
to deal with these first. In particular `x' * B * y` means `(x' * B) * y`
for an ordinary colum-major `B::Matrix`. This is often equivalent to `dot(x, B, y)`,
Copy link
Member

Choose a reason for hiding this comment

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

colum -> column

Copy link
Member

Choose a reason for hiding this comment

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

Hm, why not be explicit: "For scalar eltypes, this is equivalent to dot(x, B, y)" or something in that direction.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Maybe "Unlike dot(), ..."? I don't want to dwell on pinning down exactly what types they agree or disagree on.

Copy link
Contributor Author

@mcabbott mcabbott Jun 7, 2021

Choose a reason for hiding this comment

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

As an aside, trying to puzzle out what recursive dot means in the 3-arg case, is this the intended behaviour?

julia> x = [rand(2,2) for _ in 1:3]; A = [rand(2,2) for _ in 1:3, _ in 1:3]; y = [rand(2,2) for _ in 1:3];

julia> dot(x,A,y)
7.411848453027886

julia> @tullio _ := x[i][b,a] * A[i,j][b,c] * y[j][c,a]
7.4118484530278845

julia> @tullio _ := tr(x[i]' * A[i,j] * y[j])
7.411848453027887

The action on the component matrices looks like a trace of a matrix product.

although it allocates an intermediate array.

If the first or last factor is a number, this will be fused with the matrix
multiplication, using 5-arg [`mul!`](@ref).

See also [`dot`](@ref), [`muladd`](@ref).
See also [`muladd`](@ref), [`dot`](@ref).

!!! compat "Julia 1.7"
These optimisations require at least Julia 1.7.

# Examples
```
julia> A, B, C = randn(100,10), randn(10,100), randn(100,10);

julia> using BenchmarkTools

julia> @btime ($A * $B) * $C; # slow
13.500 μs (3 allocations: 86.14 KiB)
mcabbott marked this conversation as resolved.
Show resolved Hide resolved

julia> @btime $A * ($B * $C); # fast
1.892 μs (2 allocations: 8.81 KiB)

julia> @btime $A * $B * $C; # automatic
1.896 μs (2 allocations: 8.81 KiB)
```
"""
*(A::AbstractMatrix, B::AbstractMatrix, x::AbstractVector) = A * (B*x)

Expand Down