Fix initialization of W for static arrays

[andreasnoack]

To trigger this, I'm using LabelledArrays in the AD tests I recently added. I'm not sure the ArrayInterface.lu_instance approach can be made completely robust and, at least, it isn't right for LabelledArrays right now. The setup costs of computing the LU once should be small. Generally, I think it would be useful if initialization was based on the actual values since it's so easy to get the types wrong when you do it heuristically instead of using the instances.

[andreasnoack]

I should add that a tricky complication here is that the definition in https://github.com/JuliaDiff/ForwardDiff.jl/blob/c4487d5620983581f6a63284cca3599b8dacf203/src/jacobian.jl#L97-L105 isn't sound so the issue here is only triggered by the tests when LabelledArrays is loaded before ForwardDiff. Ref JuliaDiff/ForwardDiff.jl#468

[ChrisRackauckas]

The setup costs of computing the LU once should be small.

No. With auto switch methods (the default in DiffEq), this can be an order of magnitude more than the entire computation if it's a large non-stiff system. lu_instance was made to fix these performance issues, so I don't think we should just revert it given the known history on this subject unless something very major changed. It would be best to fix lu_instance for what we can and make it a primitive that's supported by every possible array type (hence it's part of ArrayInterface.jl).

at least, it isn't right for LabelledArrays right now

How so? The Jacobian is just a static array and static arrays are fine?

[andreasnoack]

With auto switch methods (the default in DiffEq), this can be an order of magnitude more than the entire computation if it's a large non-stiff system.

Fair enough. I'll try to patch the bug in a different way.

It would be best to fix lu_instance for what we can and make it a primitive that's supported by every possible array type

I don't think is a good approach. It's too easy to make errors when you try to guess the return types like this. The only robust way of doing this is to base it on the instances. If computing the instances of the stiff solver is too costly for non-stiff problems then the initialization of the stiff solver should probably be made lazy somehow.

The Jacobian is just a static array and static arrays are fine?

The Jacobian of an SLVector valued ODE is an SLMatrix, not an SMatrix (with the caveat mentioned in #1269 (comment)).

[ChrisRackauckas]

I don't think is a good approach. It's too easy to make errors when you try to guess the return types like this. The only robust way of doing this is to base it on the instances. If computing the instances of the stiff solver is too costly for non-stiff problems then the initialization of the stiff solver should probably be made lazy somehow.

Making the cache construction lazy would be really cool since it would allow non-stiff ODEs with auto-switch to essentially have no overhead. Right now they always build the caches for the Jacobians and such even if they aren't used.

The Jacobian of an SLVector valued ODE is an SLMatrix, not an SMatrix (with the caveat mentioned in #1269 (comment)).

That's odd?

using LabelledArrays
x = SLVector(a=1,b=2,c=3)
x.*x'


3×3 SArray{Tuple{3,3},Int64,2,9} with indices SOneTo(3)×SOneTo(3):
 1  2  3
 2  4  6
 3  6  9

[andreasnoack]

Just for the reference, I figured out what was going on here. For all but the 1x1 case, the Jacobian of the ODE system will be a SMatrix when the ODE is SLVector values but because of the branching in https://github.com/SciML/LabelledArrays.jl/blob/0998b0b881bf994645a97aea675993d91614b05d/src/slarray.jl#L66-L70 the 1x1 case ends up as an SLArray.