Add explicit arithmetics again to LowerTriangularArrays

[maximilian-gelbrecht]

See discussion in #503

[milankl]

List of in-place operations we would want to define manually

• A += B here https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/6860c05ce1fe3c15d1e1e506f2af28ec574c5746/src/dynamics/tendencies.jl#L243 and here https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/6860c05ce1fe3c15d1e1e506f2af28ec574c5746/src/dynamics/tendencies.jl#L703

• A = -A - B here https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/6860c05ce1fe3c15d1e1e506f2af28ec574c5746/src/dynamics/tendencies.jl#L279

• A += c*B here https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/6860c05ce1fe3c15d1e1e506f2af28ec574c5746/src/dynamics/tendencies.jl#L733 and here https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/6860c05ce1fe3c15d1e1e506f2af28ec574c5746/src/dynamics/geopotential.jl#L125

• A .*= c here (and other places in same file) https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/6860c05ce1fe3c15d1e1e506f2af28ec574c5746/src/SpeedyTransforms/spectral_gradients.jl#L123

[maximilian-gelbrecht]

I added the add!(C,A,B) and add!(C,A) that does C.data .+= A.data .+ B.data

Sorry for repeating this discussion again, but to be honest, the more I think about it, the less I like the non-performant broadcast, just because they are so many special cases that we'd have to cover to really have good performance everywhere.

I don't have much time this week (public holiday in Germany and so on), but I could be working on it in some way next week.

There are a few possibilities:

• The kernel one: Do we actually want to do use a broadcast here later anyway? Maybe we just put those operations in a KernelAbstractions kernel anyway, then we don't need to think about broadcasting at all
• The pragmatic one: Just use broadcasts with .data in the dynamical core where we don't have an explicit function defined
• The fancy pragmatic one: Write a macro that does the thing above automatically (@L. \$L += \$A -> @. L.data += A.data) (macros don't know the types of variables, so the \$ is used as a marker, similar to how you can escape functions in the @. macro)
• Make the broadcast more performant: I was thinking about if there's a way to write a more performant copyto function using linear indexing, but I don't really know yet if that's possible. Getting to equal performance as just using .data will be quite hard though
• Make it <:AbstractArray{T,N-1} (but keep both indexing styles), solves the broadcasting performance issue, but makes all kind of operations with regular N-dim matrices/arrays quite messy, would need a lot of extra definitions/functions there. However, those are not used in the dynamical core, this is more an interfacing issue. Probably it would be fine to drop all kind of arithmetics with regular N-dim arrays and just provide a constructor for conversion.

As before, I don't really have a strong opinion on this, but tend to "the pragmatic one". Is your list above actually exhaustive? Because if it really is only used at these few occasion, some kind of pragmatic solution is probably the most time-efficient solution for us.

But, if you have a strong opinion on this, I am happy to go with your choice.

[milankl]

Sorry for repeating this discussion again, but to be honest, the more I think about it, the less I like the non-performant broadcast, just because they are so many special cases that we'd have to cover to really have good performance everywhere.

I agree with you! But I'm actually in favour of this one

Make it <:AbstractArray{T,N-1} (but keep both indexing styles), solves the broadcasting performance issue, but makes all kind of operations with regular N-dim matrices/arrays quite messy, would need a lot of extra definitions/functions there. However, those are not used in the dynamical core, this is more an interfacing issue. Probably it would be fine to drop all kind of arithmetics with regular N-dim arrays and just provide a constructor for conversion.

because then at least from an internal perspective it's all consistent and we're not abusing any subtyping (like a vector that's a matrix as we currently do). We'd then only need to define

• [l,m] indexing as additional method for convience
• and maybe some matrix-like printing for the matrix case to see the triangular structure, this might come in handy for testing too!

Is your list above actually exhaustive?

I believe so! Checked everything, maybe I missed one occasion, but it's not just examples it's supposed to be all cases.

[maximilian-gelbrecht]

Yeah, I think that's good as well.

I have some time tomorrow and can start working on a draft. The basic change only takes little work, it's more about adjusting tests, some interfacing etc.

But would you agree that we don't support for any kind of operations with regular AbstractArray{T,N} then, and just provide constructors/conversion Array/`full``... for it?

[milankl]

But would you agree that we don't support for any kind of operations with regular AbstractArray{T,N} then, and just provide constructors/conversion Array/`full``... for it?

Yeah I think that's fine. Operations between LTAs should be possible, but I don't see much of a new for LTA + Array or LTA * Vector for example.

Do you want to merge this pull request first though? Looks like it wouldn't hurt to have them?

[maximilian-gelbrecht]

Yeah, we can merge this. The next PR will partially revert it again, but maybe I need more time than I think right now. So merging it now makes some sense and won't hurt.