ComponentArrays.jl
ComponentArrays.jl copied to clipboard
Store the diagonal of Diagonal
I was wondering whether or not it would be a good idea to special-case the type Diagonal
to only store the internal vector representing the diagonal of the matrix? I'm not familiar enough with the internals of ComponentArrays to determine whether or not this would be a straightforward addition or not. Of course, this opens up for the question whether or not all structured matrices should have a special representation, and if there would be a way to handle them all uniformly.
ComponentArray
s used to be able to hold any type of internal array, including Diagonal
with its internal storage of a single vector, but I had to change this to get sane linear algebra method handling.
Since the LinearAlgebra
methods only dispatch on whether something is a subtype of StridedArray
/DenseArray
, you have to either restrict array wrapper types like ComponentArray
s to only wrap StridedArray
/DenseArray
s so they can be passed to BLAS routines by pointer or manually overload all of the LinearAlgebra methods. I used to have it set up where I wrote all of the linear algebra methods by hand, but this led to a method ambiguity hell. Every new type of AbstractArray
that someone wanted to multiply/divide/whatever with ComponentArray
s needed a special method overload. After a few issues were filed, I decided to give up on that method and just force the inner arrays to be strided/dense.
Another approach I tried was making ComponentArray
be a union of StridedComponentArray
and NonStridedComponentArray
concrete realizations. This way it would just dispatch to the correct one based on what the underlying array type was. Unfortunately, this made getproperty
/getindex
unable to be inferred and performance dropped off a cliff. I basically have no idea what makes constant propagation work or not and a lot of the new features I've tried to implement have failed because of this. I still think there might be a way to handle the Strided
/NonStridedComponentArray
thing, by writing separate methods for getindex
/getproperty
for each type and not trying to dispatch in the return of the @generated
function.
As far as features I'd like to implement, this is at the top of my list since it would allow for StaticArray
s and CuArray
s, which I've wanted to have back for a while now. I would still have to write manual LinearAlgebra methods for the NonStridedComponentArray
s, but I'm okay with that because at least the "normal" version wouldn't need them and I'd probably get fewer issues coming in.
But the bigger issue is that we need a version of LinearAlgebra
that dispatches on traits rather than supertypes. If you look through all of the other array wrapper type libraries, they all tend to have the same issues filed for method ambiguities. There are a few issues filed in the std repo about it, but they seemed to have lost traction. Maybe I'll try to revive that discussion.
Thanks for your detailed explanation. I recognize the issues you're describing with method ambiguities, they appear all too often when implementing a type that behaves like another, standard type.
Thanks for this package, it makes replacing simulink/modelica a breeze while keeping the model as code rather than a bunch of boxes and lines in a gui.