Sparse kronecker product

I have a piece of code that computes the kronecker product between sparse matrices in csr format in the following way

import scipy.sparse as sprs
import numpy as np

A          = np.random.rand( 4,4 )
nids_left = 10
nids_right= 10  

idm      = sprs.identity( 2**(nids_left) ,dtype=A.dtype, format="csr")           
A        = sprs.kron( idm, A, format="csr" )
idm      = sprs.identity( 2**(nids_right) ,dtype=A.dtype, format="csr")           
A        = sprs.kron( A, idm, format="csr" )

The nids_left and nids_right indicate the total number of 2x2 identities to the left and right of the a A matrix; in other words I want to compute IxIxI...xIxAxIx...xI in an efficient way, but I just realized that scipy.sparse uses a single thread when doing the kronecker product. Is there a more efficient way to do this? In particular, if nids_left and nids_right become larger, memory considerations aside.