Fast Kronecker matrix multiplication

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.