Let $\rho_{ABCD}$ be a sparse matrix of 4 systems each in a $d$-dimensional Hilbert space.
For $d<7$ in a reasonable time (few seconds) I able to perform the partial trace $\rho_{AD}$ using the code proposed in http://www3.imperial.ac.uk/people/m.tame/research. I would need an efficient algorithm for calculating $\rho_{AD}$ where $d\geq7$. The algorithm of the site above does not exploit any properties of the matrix and it requires a lot of permutations and rearrangements.
Do you know an efficient algorithm for calculating the partial trace of qudit which uses the fact that the matrix is sparse? It would also be interesting if the algorithm can take advantage of parallel computation.
Thank you very much in advance for your answers.
Regards,
Silvio