Generalized pair couplings

[kbarros]

This PR provides initial support for #63, with some differences:

1. It uses the built-in LinearAlgebra.kron as the tensor product operation.
2. Given a general interaction of the form $\sum_i A_i \otimes B_i$, this PR will "extract" the bilinear spin-spin exchange, and the scalar biquadratic parts, if present. This is done to reduce the number of terms present in the sum of tensor products, which is empirically the slowest part of a simulation.
3. The generalized pair couplings work in the classical SU(N) dynamics but not yet spin wave theory. This will be addressed in a future PR.
4. Currently :dipole mode simulations cannot handle anything but scalar biquadratic. That will be addressed in another PR.

[Lazersmoke]

How many terms in the SVD sum do you need in the best/worst/typical case?

[kbarros]

This illustrates the behavior on the this branch. The user interface is subject to change, and that's the current blocker for this PR.

using Sunny

S1 = spin_matrices(; N=3)
S2 = spin_matrices(; N=4)
Si, Sj = Sunny.to_product_space(S1, S2)

D = Si'*randn(3,3)*Sj
@assert length(Sunny.svd_tensor_expansion(D, 3, 4)) == 3

D = (Si'*randn(3,3)*Sj)^2
@assert length(Sunny.svd_tensor_expansion(D, 3, 4)) == 9

# In practice, though, parameters to `set_pair_coupling!` will pass through this
# decomposition function:
bilin, tensordec = Sunny.decompose_general_coupling(D, S1, S2; fast=true)
@assert length(tensordec.data) == 6  # Just 6 terms left after extracting 3x3 dipole-dipole part

# As an optimization, `decompose_general_coupling`, could also extract the 5x5
# quadrupole-quadrupole part. For this particular example, there would be no terms
# left in the "genereal" expansion. (Not yet implemented.)

# bilin, biquad, tensordec = Sunny.decompose_general_coupling(D, S1, S2; fast=true)