Fast Kronecker matrix multiplication

Fast Kronecker matrix multiplication for matrices of any size

Fast Kronecker matrix multiplication, for both full and sparse matrices
of any size. Never computes the actual Kronecker matrix and omits
multiplication by identity matrices.
y = kronm(Q,x) computes
y = (Q{k} kron ... Q{2} kron Q{1})*x
If Q contains only two matrices and x is a vector, the code uses the
identity
( Q{2} kron Q{1} )*vec(X) = vec(Q{1}*X*Q{2}'),
where vec(X)=x. If Q contains more than two matrices and/or if x has more
than one column, the algorithm uses a generalized form of this identity.
The idea of the algorithm is to see x as a multi-dimensional array and to
apply the linear maps Q{i} separately for each dimension i.
Acknowledgement:
This code follows the same idea as 'kronmult' by Paul G. Constantine &
David F. Gleich (Stanford, 2009). However, I avoid loops and allow for
non-square inputs Q{i}. I have also included the special treatment for
identity matrices.