block diagonal interaction terms
Some models have a block diagonal structure in the interaction terms, similar to the already exploited N_fl structure. However, the non-interacting part may not be block diagonal such that N_fl cannot be used. This branch aims to implement this structure!
I would propose to redefine the operator type to contain a list of the "dense" blocks of the operators. I.e., P(:) --> P(:,b); O(:,:) --> O(:,:,b) and so on. b should be the index of the block and appear as the last index for memory efficiency as with n_fl. The rest of the code should then be updated accordingly.
It might also be convenient, to use the number of blocks as an optional argument when the operator is allocated, ideally in a way that the existing model without that structure can stay as they are. This way, the more experienced user can make use of this option whereas a "new" user doesn't have to care about this in the beginning.