Partial trace algorithm (original)

In general, is there a partial trace algorithm (ideally for systems of any size) that can be coded using basic matrix operations found in software like Mathematica or Maple? All of the methods I'm aware of seem to be suited for humans, and I am struggling to think of a way of implementing them in any programming language I know. I would appreciate any kind of help with this.

Thank you.

• This might be better suited at scicomp.stackexchange than here. Oct 24 '13 at 20:06
• @KyleKanos: Oh, thanks. It's usually something that physicists deal with so I thought that asking here would be appropriate.
– sps
Oct 24 '13 at 20:15

$\newcommand{\ket}{\lvert#1\rangle}$ I don't know why the methods are more suited for humans than for computer, and you can easily find a MatLab implementation by googling, for example on Toby Cubitt's webpage. (I haven't checked the algorithm, but I trust Toby to get such a thing correctly.)
I'll try below to give you an explanation “suited for computers”. Let $\rho_{AB}$ be the density matrix you're interested in, with the $d_A d_B$ base vectors in the following order : $$\ket{1,1},\ket{1,2},\dots,\ket{1,d_B},\ket{2,1},\dots,\ket{d_A,d_B}$$
To obtain $\rho_B=\text{Tr}_A(\rho_{AB})$, you just have to add the $d_A$ diagonal submatrices of $\rho_AB$. Something like $\rho_B=\sum_{j=0}^{d_A-1}\rho_{AB}[jd_B+1\dots jd_B+d_A;jd_B+1\dots jd_B+d_A]$.
To get $\rho_A$, you can reorder the indices to go back to the previous algorithm. Or you directly take the relevant submatrices, but they are somehow scattered. $\rho_A[k;l]=\sum_{i}^{d_B}\rho_{AB}[(k-1)d_B+i;(l-1)d_B+i]$