Store the diagonal of Diagonal

[baggepinnen]

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.

[jonniedie]

ComponentArrays 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 ComponentArrays to only wrap StridedArray/DenseArrays 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 ComponentArrays 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 StaticArrays and CuArrays, which I've wanted to have back for a while now. I would still have to write manual LinearAlgebra methods for the NonStridedComponentArrays, 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.

[baggepinnen]

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.