Benchmark performance of starry rotation matrices

[dfm]

The rotation matrices are block diagonal and in the implementation in #17, the full dense matrix is constructed. This probably isn't what we want in all cases: for larger systems we might want to manually implement the matrix product. That being said, if we're using GPU acceleration, we may want to use the direct solve up to fairly high order. We should benchmark this on several different systems to select sensible defaults and give advice to users.