You can’t really use those implementations at scale. And they traverse the matrices in non-local order.
The trick to working with Kronecker products is to never expand them out. You use the fact that things like inverses and Cholesky factors distribute, e.g.,
\text{cholesky}(A \otimes B) = \text{cholesky}(A) \otimes \text{cholesky}(B).
This is why we didn’t put them into Stan, because Stan always creates the full matrices it returns, so if we had a Kronecker product, it’d just be brute force like the implementations linked above.
Stan does allow that—the result is a submatrix of the original matrix. But it copies everything, so it’s inefficient.