* has no method matching *(::Array{Float64,2}, ::Tridiagonal{Float64})

[Ken-B]

To replicate:

A = Tridiagonal(rand(4), rand(5), rand(4))
B = rand(5, 5)
A * B #works
B * A #MethodError
B * full(A) #works

Seems to me some fallback is missing. I'm on v0.4-rc2, haven't tried on master.

[Ken-B]

Also, more missing method for Bidiagonal:

A = Bidiagonal(rand(5), rand(4), false)
B = rand(5, 5)
A * B #MethodError
B * A #MethodError
B * full(A) #works
full(A) * B #works

[jw3126]

There are also performance issues:

n = 100

A = Tridiagonal(randn(n-1), randn(n), randn(n-1))
Ad = full(A)

B = Tridiagonal(randn(n-1), randn(n), randn(n-1))
Bd = full(B)

@time A*B  #  0.006436 seconds (13 allocations: 234.859 KB))
@time A*Bd  #  0.000037 seconds (7 allocations: 78.391 KB)
@time Ad*Bd  #  0.003669 seconds (7 allocations: 78.391 KB)

e.g. multiplication of two tridiagonal matrices is fastest if we are only given the info that the second one is tridiagonal and worst if we are informed about both matrices...

Anyway I am on it.

[timholy]

All the MethodErrors were resolved (presumably by JuliaLang/julia#15367), but the performance issues were not. Would be a great contribution!

[pkofod]

@jw3126 Great. Don't be shy to ask if you get stuck at something! And remember to work on a commit after the one Tim mentioned :)

[jw3126]

@timholy @pkofod Thanks! I have already a fast kernel for the right multiplication by something tridiagonal.

n = 100

A = randn(n,n)
B = Tridiagonal(randn(n-1), randn(n), randn(n-1))
fullB = full(B)
C = randn(n,n)

@time A_mul_B!(C, A, B)  #  0.000022 seconds (4 allocations: 160 bytes)
@time A_mul_B!(C, A, fullB)  # 0.002359 seconds (4 allocations: 160 `bytes)

However multiplication is not yet dispatched in a way that this kernel is called.
The issue is that there are currently two patterns of multiplication function around:

1. A * B calls A_mul_B!(C, A, B)
2. A * B calls full(A) * full(B)

I could fix that problem by overloading * for combinations where the right matrix is tridiagonal.

However I think this is a more general issue. I think we should try to get rid of 2) completely and substitute it by definitions like

A_mul_B!(C, A, B) = A_mul_B!(C, full(A), full(B))

Otherwise I it is easy to introduce performance bugs by overloading * and thereby accidentally preventing calls to some fast A_mul_B! for methods.
In other words we have two alternatives:

1. Have A_mul_B! methods for various type combinations
2. Have lots of * and lots of A_mul_B! methods for various type combinations

I think 2) is code duplication and should be omitted. What do you think? Is this reasonable and possible?

A related question: Is there a way to test that a call of some method (like *) results in a call of some other method (like A_mul_B!) ?

[timholy]

I'm not sure we need to convert anything to full now, if you write everything in terms of AbstractMatrix. For optimized implementations, there could be a "priority" issue (which matrix is "in charge" when you can solve the problem using sparse iteration?) that could lead to ambiguities.

JuliaLang/julia#15459 sketches one potential way out of this problem. It is a work-in-progress, and the API will likely change. (Any good ideas are of course welcome!) But the idea might be that we could dispatch to a single method when either A or B (or both) exhibit some kind of sparsity.

[pkofod]

I think @jw3126 was referring to https://github.com/JuliaLang/julia/blob/master/base/linalg/bidiag.jl#L228

[tkelman]

a tip @pkofod - when linking to lines of code, hit y and github will replace the link with one that has the commit sha in it like so https://github.com/JuliaLang/julia/blob/e588b4b456df4db9c34bec67580e43c3373a253c/base/linalg/bidiag.jl#L228 - that way when looking at this issue in a few months your link will still make sense even if the code has moved around

[pkofod]

a tip @pkofod - when linking to lines of code, hit y and github will replace the link with one that has the commit sha in it like so https://github.com/JuliaLang/julia/blob/e588b4b456df4db9c34bec67580e43c3373a253c/base/linalg/bidiag.jl#L228 - that way when looking at this issue in a few months your link will still make sense even if the code has moved around

Will do! Neat trick, thanks.

[jw3126]

@pkofod is right I had that line of code in mind. @timholy so you say for various sparse cases it does not make sense to define * by calling A_mul_B!(C, A, B) anyway?

[pkofod]

Before I realized you were working on this, I poked around, that why I figured that was the line. I think you should just remove Tridiagonal from Special, and then be sure that you have covered the different cases of multiplication using Tridiagonals?

[timholy]

@jw3126, sorry, I was not very clear. I think having * call A_mul_B! is fine (I understand there's long-term interest in changing that API, but that's really a separate issue). All I meant was that I hope we can soon get to a point where we might be able to write a single

function A_mul_B!(C::AbstractMatrix, A::AbstractSparse, B::AbstractSparse)

and never have to convert anything to full.

Right now, whether to convert to full seems to be a somewhat difficult decision: the generic fallback means you don't need to make that conversion (the generic algorithm should work for sparse inputs), but of course the performance is poorer than you'd get from using BLAS, which is an argument in favor of performing the conversion.

[andreasnoack]

Fixed by JuliaLang/julia#15505