JuliaArrays / StaticArrays.jl

Statically sized arrays for Julia
Other
761 stars 147 forks source link

All dispatches of static arrays with factorizations are missing #478

Open ChrisRackauckas opened 6 years ago

ChrisRackauckas commented 6 years ago

For example, while you can build the lu(A), every use of it errors:

using StaticArrays, LinearAlgebra
A = @SMatrix rand(3,3)
Alu = lu(A)
c = rand(3)

b = rand(3)
Alu \ b
Alu / b
ldiv!(c,Alu,b)
b = @SVector rand(3)
Alu \ b
Alu / b
ldiv!(c,Alu,b)
b = @MVector rand(3)
Alu \ b
Alu / b
ldiv!(c,Alu,b)

c = @MVector rand(3)

b = rand(3)
Alu \ b
Alu / b
ldiv!(c,Alu,b)
b = @SVector rand(3)
Alu \ b
Alu / b
ldiv!(c,Alu,b)
b = @MVector rand(3)
Alu \ b
Alu / b
ldiv!(c,Alu,b)

This is likely true for the other factorizations as well. Quick test:

Alu = qr(A)
c = rand(3)
b = rand(3)
Alu \ b

There should be a way to programmatically do this and copy Base.

xref: https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/339

andyferris commented 6 years ago

Right... this is clearly important but I haven't had time for investigating this (and since it isn't used at my workplace, it doesn't get prioritized in my personal time, sorry).

In theory, Julia is good at code composition... my opinion is that we shouldn't really need our own factorization containers, nor to define our own method specializations in order for this stuff to work. (Some cases might benefit from unrolled code for speed of course, but there shouldn't be errors). Unfortantely, Base and LinearAlgebra frequently make assumptions that are true for Array but not all methods work on a bare AbstractArray interface.

Would it be fair to say that this is something that should be done in tandem to LinearAlgebra so that the two codebases compose with each other more naturally?

cc @andreasnoack

andreasnoack commented 6 years ago

Some of the factorizations in LinearAlgebra are still too closely tied to LAPACK. In particular, it would be great to get rid of any hard-coded Arrays. However, there might still be situations where special matrix types will benefit from a special version of a factorization.

andyferris commented 6 years ago

Right, totally agreed.

I guess where I'm at is I'm happy to use the container types from LinearAlgebra (or else inherit from a LinearAlgebra.AbstractLU supertype, for example) and provide specialized methods here where they make a speedup. At the moment I'm stuck with basic questions like what is the interface guaranteed to be provided by a Factorization?