Partial trace of matrix

For a given density matrix which originates from a vector which contains four systems, i want to compute the reduced density matrix for each of the four subsystem. (i.e. I have a density matrix which has constructed from four subsystem outer products and I want to compute the partial trace of this matrix for each one of those subsystems).

I. I was wondering if there is an efficient way to do this

Answer 1

[Jean-Daniel]

http://www.dr-qubit.org/matlab.php contains all that.

Answer 2

[Nikolaos]

thx a lot.. but the problem is that this is not working for all cases. For example if you have tri-partite or four-partite systems. This is my way of making a partial trace of the initial matrix pall of dimensions 100x100 into the partial matrix pred of dimensions 10x10 according to pred(i,j)= Sum_ijkl pall(i,j,k,l)

red=zeros(10)
pd= reshape(pall, 10,10,10,10)
for ii=1:10
for jj=1:10
for mm=1:10
pred(ii,jj)=pred(ii,jj)+pd(ii,mm,jj,mm)
end
end
end