indexing into a matrix to eliminate a loop

This file generates an ergodic density matrix. I would like to get rid of the inner loop (over ii). The first line in the outer loop - F1 = F0*P; - reshuffles columns of F0 to create F1, and it is as quick as it can get.

But depending on how the GPOL matrix looks like, the inner loop can be quite slow. Now this line of command (which I think was courtesy of Matt J,

   Zi(sub2ind(size(Zi),GPOL(ii,:),1:numel(GPOL(ii,:))))=F1(ii,:);

this is already pretty smart. It is reshuffling the rows of F1, one by one. here you can look at the entire code, and i will state the question below (you may be guessing it already)

function[F2,ijj,differ] = densityiter(F0,GPOL,P);

[nh ny]=size(F0);

niter = 400;

for ijj=1:niter

     F1 = F0*P;

     Z0 = 0*ones(nh,ny);

     F2  = 0*ones(nh,ny);

     for ii=1:nh

           Zi = Z0;

           Zi(sub2ind(size(Zi),GPOL(ii,:),1:numel(GPOL(ii,:))))=F1(ii,:);

           F2 = F2 + Zi;

     end

     differ = sum(sum(abs(F2-F0)));

     F0=F2;

     if differ<0.000001;break,else,end

end

The question then is: does anyone know how to replace the loop over ii to get the same effect only faster? to work on all rows of F1 at the same time? this is the piece

       for ii=1:nh

             Zi = Z0;

             Zi(sub2ind(size(Zi),GPOL(ii,:),1:numel(GPOL(ii,:))))=F1(ii,:);

             F2 = F2 + Zi;

       end

thank you. i hope this rephrasing on the problem helps.

[Matt J]

We can talk about it, but please highlight your code and apply the

formatting button first.

[Joao Ejarque]

Dear Matt

Fair points. I will do my best.

Answer 1

[Matt J]

 e=repmat(1:size(GPOL,2),nh,1);
 subs=[GPOL(:),e(:)];

 F2=accumarray(subs,F1(:),[nh,ny]);

[Joao Ejarque]

Dear Matt. Thank you. I will check both answers received. I will just add that I never cease to be amazed by how good you guys are with this. Regards, J.

[Joao Ejarque]

works brilliantly. i think i am going to pick your brain again on another "getting rid of a loop" problem soon. I promise not to make a habit of it though. regards. J.

Answer 2

[Amit]

I am assuming that GPOL and F1 are of same size. Then instead of:

for ii=1:nh
           Zi = Z0;
           Zi(sub2ind(size(Zi),GPOL(ii,:),1:numel(GPOL(ii,:))))=F1(ii,:);
           F2 = F2 + Zi;
     end

you can so something like this:

Zi = Z0;
[m n] = size(GPOL);
F1tmp = F1';
Zi(sub2ind(size(Zi),reshape(GPOL',1,m*n),repmat(1:numel(GPOL(1,:)),1,nh))) = F1tmp(1:m*n);

This might need some tweaking. If you post a sample data then I can verify this.

[Joao Ejarque]

Dear Amit. Thank you. I will check it and get back to you. J.

[Joao Ejarque]

Dear Amit. This does not work. But Matt J.´s answer below works so you can check it out. Regards. J.