Thread Subject:

help on vectorizing this

Hello,
I need some help on vectorizing the following loop

for t=1:num_probes
    for i=1:num_ebc
        f(node_ebc(i),t)=prop_ebc(i,t);
        k(node_ebc(i),node_ebc(i))=1;
        for j=1:num_nodes
if node_ebc(i)~=j
                  f(j,t)=f(j,t)-k(j,node_ebc(i))*prop_ebc(i,t);
                  k(node_ebc(i),j)=0;
                  k(j,node_ebc(i))=0;
              end
        end
    end
end

What I am thinking is like this
Rearrange oder so I can avoid the
"if node_ebc(i)~=j"


f=f-k(:,node_ebc).*prop;
f(node_ebc,:)=prop;


k(cm.node_ebc,:)=0;
k(:,cm.node_ebc)=0;
k(sub2ind(size(k),cm.node_ebc,cm.node_ebc))=1;


What do you think?

Marios Karaoulis <marios.karaoulis@gmail.com> wrote in message <d6bed29a-31cb-400b-a836-aed0e490610d@c13g2000vbr.googlegroups.com>...
> What I am thinking is like this
> Rearrange oder so I can avoid the
> "if node_ebc(i)~=j"
>
>
> f=f-k(:,node_ebc).*prop;
> f(node_ebc,:)=prop;
>
>
> k(cm.node_ebc,:)=0;
> k(:,cm.node_ebc)=0;
> k(sub2ind(size(k),cm.node_ebc,cm.node_ebc))=1;
>
>
> What do you think?

I think I must add that
f=[num_nodesXnum_probes];
K=[num_nodes,num_nodes];

Marios Karaoulis <marios.karaoulis@gmail.com> wrote in message <73f25889-8b76-4aa7-abfe-a6bb1923b88e@t11g2000vbc.googlegroups.com>...
> Hello,
> I need some help on vectorizing the following loop
>
> for t=1:num_probes
> for i=1:num_ebc
> f(node_ebc(i),t)=prop_ebc(i,t);
> k(node_ebc(i),node_ebc(i))=1;
> for j=1:num_nodes
> if node_ebc(i)~=j
> f(j,t)=f(j,t)-k(j,node_ebc(i))*prop_ebc(i,t);
> k(node_ebc(i),j)=0;
> k(j,node_ebc(i))=0;
> end
> end
> end
> end
>

Let me try and explain what I see, to insure it's correct:
You have a number of probes you're going to loop over. From all of these probes you are generating to matrices f,k which are dependent on none_ebc. You want the resulting f,k to be a 2-dimensional matrix of some size generated by the inner computation from all positions of all probes?

Once again I think this would be easiest if we had some type of sample data or at least sizes to play with.

Anyway:
Since the k matrix isn't preallocated it's being filled with zeros as you increment. I would recommend preallocating it with zeros, and then just adding ones where you have to. Since it appears the ones are all on the diagonal this would be very easy and I think you came close with your try in the second post:
k = zeros(the_final_k_size);
k(sub2ind(size(k),node_ebc,node_ebc))=1;

As for f:
It's doable with a call to meshgrid()

Thanks (again!) for your reply. After some thinking this is what I did
so far. Although I do not write it in the forum, all of my matrices
are pre-allocated
M=num_nodes % like 20000;
N=num_probes; % 200
K=num_ebc; % 2000

k=zeros(M,M);
f=zeros(M,N);
prop_ebc=zeros(K,N);
node_ebc=zeros(K,1); %This is matrix that has the index on which
nodes I need to do some math.
...... %rest of code where k, node_ebc, prop_ebc are no longer zeros


for t=1:N
    for i=1:K
        f(:,t)=f(:,t)-k(:,node_ebc(i)).*prop_ebc(i,t);
    end
end


f(node_ebc,:)=prop_ebc;
k(node_ebc,:)=0;
k(:,node_ebc)=0;
k(sub2ind(size(k),node_ebc,node_ebc))=1;


The only loop I have now, is very time-consuming.

Marios Karaoulis <marios.karaoulis@gmail.com> wrote in message <05517269-cbbf-46db-912b-73ec53b4c08b@37g2000prx.googlegroups.com>...
> Thanks (again!) for your reply. After some thinking this is what I did
> so far. Although I do not write it in the forum, all of my matrices
> are pre-allocated
> M=num_nodes % like 20000;
> N=num_probes; % 200
> K=num_ebc; % 2000
>
> k=zeros(M,M);
> f=zeros(M,N);
> prop_ebc=zeros(K,N);
> node_ebc=zeros(K,1); %This is matrix that has the index on which
> nodes I need to do some math.
> ...... %rest of code where k, node_ebc, prop_ebc are no longer zeros
>
>
> for t=1:N
> for i=1:K
> f(:,t)=f(:,t)-k(:,node_ebc(i)).*prop_ebc(i,t);
> end
> end
>
>
> f(node_ebc,:)=prop_ebc;
> k(node_ebc,:)=0;
> k(:,node_ebc)=0;
> k(sub2ind(size(k),node_ebc,node_ebc))=1;
>
>
> The only loop I have now, is very time-consuming.

Actually the loop isn't bad compared to a bsxfun approach on the whole thing. However, you have no reason to operate on the whole of (k(:,nodes)) because it's VERY sparse and there's only one value per row.

Here's three approaches with results:

M=2000; %num_nodes % like 20000;
N=200;%num_probes; % 200
K=2000;%num_ebc; % 2000

k=eye(M); %zeros(M,M);
f=zeros(M,N);
prop_ebc=rand(K,N);
node_ebc=[1:K].'; %This is matrix that has the index on which
%nodes I need to do some math.
%...... %rest of code where k, node_ebc, prop_ebc are no longer zeros

%% Yours
tic
for t=1:N
    for i=1:K
        f(:,t)=f(:,t)-k(:,node_ebc(i)).*prop_ebc(i,t);
    end
end
toc

%% bsxfun on whole k matrix
tic
f2 = zeros(size(f));
kpart = -k(:,node_ebc);
for ii = 1:N
    f2(:,ii) = sum(bsxfun(@times,kpart,prop_ebc(:,ii)),2);
end
toc

%% bsxfun on extracted parts of k.
tic
kpart = -k(:,node_ebc);
idk = find(kpart);
f3 = bsxfun(@times,prop_ebc,kpart(idk));
toc
isequal(f,f2,f3)

%{
Elapsed time is 2.750416 seconds.
Elapsed time is 10.143211 seconds.
Elapsed time is 0.106752 seconds.

ans =

     1
}

This versions works only if M==K, otherwise it says
??? Error using ==> bsxfun
Non-singleton dimensions of the two input arrays must match each other.


Let me rephrase the problem with some initial data.


clear all
num_nodes=100;
num_probes=10;
num_ebc=20;
prop_ebc=rand(num_ebc,num_probes);

k=rand(num_nodes,num_nodes);
% k=eye(num_nodes);
f=zeros(num_nodes,num_probes);
prop_ebc=rand(num_ebc,num_probes);
node_ebc=randi([1 num_nodes],num_ebc,1);

k2=k;
f2=f;

k3=k;
f3=f;

%%--------------------First Attempt ----------------------------------------
tic
for t=1:num_probes
% k_tmp=k2;
    for i=1:num_ebc
        f(node_ebc(i),t)=prop_ebc(i,t);
        k(node_ebc(i),node_ebc(i))=1;
        for j=1:num_nodes
if node_ebc(i)~=j
                  f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
                  k(node_ebc(i),j)=0;
                  k(j,node_ebc(i))=0;
              end
        end
    end
end
toc



%%--------------------Second Attempt ---------------------------------------
tic
for t=1:num_probes
    for i=1:num_ebc
        f2(:,t)=f2(:,t)-k2(:,node_ebc(i)).*prop_ebc(i,t);
    end
end
   

f2(node_ebc,:)=prop_ebc;
k2(node_ebc,:)=0;
k2(:,node_ebc)=0;
k2(sub2ind(size(k2),node_ebc,node_ebc))=1;
toc

isequal(k,k2)
isequal(f,f2)



%%---------------------Third Attempt ---------------------------------------

tic
kpart = -k3(:,node_ebc);
idk = find(kpart);


f3 = ( bsxfun(@times,prop_ebc,kpart(idk)) );
toc




isequal(f2,f3)




With my version (2nd attempt) I k==k2, but f~=f2.
I thought it worked, cause in my case
k2(j,node_ebc(i))*prop_ebc(i,t)=0
so it doesn't effect me. But in some other case, this might not be zero.

"Marios Karaoulis" <marios_hellas_23@yahoo.gr> wrote in message <i9a3r7$jfk$1@fred.mathworks.com>...
> This versions works only if M==K, otherwise it says
> ??? Error using ==> bsxfun
> Non-singleton dimensions of the two input arrays must match each other.
>
>
> Let me rephrase the problem with some initial data.
>
>
> clear all
> num_nodes=100;
> num_probes=10;
> num_ebc=20;
> prop_ebc=rand(num_ebc,num_probes);
>
> k=rand(num_nodes,num_nodes);
> % k=eye(num_nodes);
> f=zeros(num_nodes,num_probes);
> prop_ebc=rand(num_ebc,num_probes);
> node_ebc=randi([1 num_nodes],num_ebc,1);
>
> k2=k;
> f2=f;
>
> k3=k;
> f3=f;
>
> %%--------------------First Attempt ----------------------------------------
> tic
> for t=1:num_probes
> % k_tmp=k2;
> for i=1:num_ebc
> f(node_ebc(i),t)=prop_ebc(i,t);
> k(node_ebc(i),node_ebc(i))=1;
> for j=1:num_nodes
> if node_ebc(i)~=j
> f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
> k(node_ebc(i),j)=0;
> k(j,node_ebc(i))=0;
> end
> end
> end
> end
> toc
>
>
>
> %%--------------------Second Attempt ---------------------------------------
> tic
> for t=1:num_probes
> for i=1:num_ebc
> f2(:,t)=f2(:,t)-k2(:,node_ebc(i)).*prop_ebc(i,t);
> end
> end
>
>
> f2(node_ebc,:)=prop_ebc;
> k2(node_ebc,:)=0;
> k2(:,node_ebc)=0;
> k2(sub2ind(size(k2),node_ebc,node_ebc))=1;
> toc
>
> isequal(k,k2)
> isequal(f,f2)
>
>
>
> %%---------------------Third Attempt ---------------------------------------
>
> tic
> kpart = -k3(:,node_ebc);
> idk = find(kpart);
>
>
> f3 = ( bsxfun(@times,prop_ebc,kpart(idk)) );
> toc
>
>
>
>
> isequal(f2,f3)
>
>
>
>
> With my version (2nd attempt) I k==k2, but f~=f2.
> I thought it worked, cause in my case
> k2(j,node_ebc(i))*prop_ebc(i,t)=0
> so it doesn't effect me. But in some other case, this might not be zero.

Oh! I thought k was just a sparse identity matrix not a boundary applied stiffness matrix. Are you sure 'f' is right? It seems to me that f2 should be right. If you inspect f,f2, f has values everywhere, where f2 is zeros on the rows of k that had been zeroed out.

> Oh! I thought k was just a sparse identity matrix not a boundary applied stiffness matrix. Are you sure 'f' is right? It seems to me that f2 should be right. If you inspect f,f2, f has values everywhere, where f2 is zeros on the rows of k that had been zeroed out.


Not it is not right. f and f2 are equal only to the (node_ebc) indexing, since i explicit assign values there.

In all other indexes, it is not.

>
> %%--------------------Second Attempt ---------------------------------------
> tic
> for t=1:num_probes
> for i=1:num_ebc
> f2(:,t)=f2(:,t)-k2(:,node_ebc(i)).*prop_ebc(i,t);
> end
> end
>

FYI, the above code is usually called *matrix multiplication*

f2 = -k2(:,node_ebc)*prop_ebc

Bruno

> %%--------------------First Attempt ----------------------------------------
> tic
> for t=1:num_probes
> % k_tmp=k2;
> for i=1:num_ebc
> f(node_ebc(i),t)=prop_ebc(i,t);
> k(node_ebc(i),node_ebc(i))=1;
> for j=1:num_nodes
> if node_ebc(i)~=j
> f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
> k(node_ebc(i),j)=0;
> k(j,node_ebc(i))=0;
> end
> end
> end
> end
> toc

I don't think this is doing what you want. When i =1, that node_ebc position in f, is not overwritten; however, all of the others are. What you would want I think is:
> for j=1:num_nodes
> if all(node_ebc~=j)
> f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
> en
> end

And thus it checks every place in node_ebc, not just the current one.

Let me know if this is correct.

> I don't think this is doing what you want. When i =1, that node_ebc position in f, is not overwritten; however, all of the others are. What you would want I think is:
> > for j=1:num_nodes
> > if all(node_ebc~=j)
> > f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
> > en
> > end
>
> And thus it checks every place in node_ebc, not just the current one.
>
> Let me know if this is correct.

I want to f(node_ebc(1)) to have a specific value, it is the Dirichlet boundaries conditions. If this node is not on the boundary (eg do not belong to node_ebc), make this calculation.

The initial part of the code is correct (1st attempt). It is not my code, is a typical FEM code.

"Marios Karaoulis" <marios_hellas_23@yahoo.gr> wrote in message <i9aa9f$s2s$1@fred.mathworks.com>...
> > I don't think this is doing what you want. When i =1, that node_ebc position in f, is not overwritten; however, all of the others are. What you would want I think is:
> > > for j=1:num_nodes
> > > if all(node_ebc~=j)
> > > f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
> > > en
> > > end
> >
> > And thus it checks every place in node_ebc, not just the current one.
> >
> > Let me know if this is correct.
>
> I want to f(node_ebc(1)) to have a specific value, it is the Dirichlet boundaries conditions. If this node is not on the boundary (eg do not belong to node_ebc), make this calculation.
>
> The initial part of the code is correct (1st attempt). It is not my code, is a typical FEM code.

It is not right, and your boundary conditions are not maintained. You set the boundary conditions in this line:
   f(node_ebc(i),t)=prop_ebc(i,t);
and then you overwrite them in this line:
   f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);

The reason is because you check to see if it's a boundary point against ONE of the node_ebc, the one in the ith position. Thus if it's in ANY other position, so say 19 of your 20 node_ebcs the following line overwrites it. Simply look at the example data you provided earlier: look at what it should be (prob_ebc) and you'll see that NONE of the boundary points are equal to these.

Using this as your conditional statement:
if all(node_ebc~=j)
...
end
fixes this problem. Try it and see, you'll notice your BCs are maintained in every coordinate position they should be.

-Sean

> Using this as your conditional statement:
> if all(node_ebc~=j)
> ...
> end
> fixes this problem. Try it and see, you'll notice your BCs are maintained in every coordinate position they should be.

Working off this assumption, this will speed the code up a ton. It removes the inner loop and uses just the outer two. Here's a bench test:
%%%%
num_nodes=1000;
num_probes=100;
num_ebc=200;

k=rand(num_nodes,num_nodes);
f=zeros(num_nodes,num_probes);
prop_ebc=rand(num_ebc,num_probes);
node_ebc=randi([1 num_nodes],num_ebc,1);

f2=f;

%% --------------------First Attempt ----------------------------------------
tic
for t=1:num_probes
    for i=1:num_ebc
        f(node_ebc(i),t)=prop_ebc(i,t);
        for j=1:num_nodes
            if all(node_ebc~=j)
                  f(j,t)=f(j,t)-k(j,node_ebc(i))*prop_ebc(i,t);
            end
        end
    end
end
toc



%% --------------------Second Attempt ---------------------------------------
tic
[nn tt] = meshgrid(1:num_probes,node_ebc);
f2(sub2ind(size(f2),tt,nn)) = prop_ebc;
jlist = setdiff(1:num_nodes,node_ebc);

for t=1:num_probes
    for i=1:num_ebc
       f2(jlist,t)= f2(jlist,t)-sum(k(jlist,node_ebc(i)).*prop_ebc(i,t),2);
    end
end


k(node_ebc,:)=0;
k(:,node_ebc)=0;
k(sub2ind(size(k),node_ebc,node_ebc))=1;

toc

isequal(f,f2)

%{
Elapsed time is 60.797902 seconds.
Elapsed time is 0.604357 seconds.

ans =

     1
%}

> %% --------------------First Attempt ----------------------------------------
> tic
> for t=1:num_probes
> for i=1:num_ebc
> f(node_ebc(i),t)=prop_ebc(i,t);
> for j=1:num_nodes
> if all(node_ebc~=j)
> f(j,t)=f(j,t)-k(j,node_ebc(i))*prop_ebc(i,t);
> end
> end
> end
> end
> toc
>


 
I see what do you mean.
This code was written from a fortan code which goes like this

do 9 t=1,num_probes
     do 8 i=1,num_ebc
           f(node_ebc(i),t)=prop_ebc(i,t);
              do 7 j=1,num_nodes
                  if (j.eq.node_ebc(i)) go to 7
                     f(j)=f(j)-k(j,node_ebc(i))*prop_ebc(i)
7 continue
8 continue
9 continue

 I didn't want to use go to 7, so i replace it with the if (node_ebc(j)~=j) ,

Obviously the end statement in my code is inside the loop.

I would never have found it, it war pure luck that in my case
k(j,node_ebc(i))*prop_ebc(i,t)=0
in all times, so it didn't effect my data.

Thank you.

> if all(node_ebc~=j)


By adding this, now k is not correct with the one I used
k(node_ebc,:)=0;
k(:,node_ebc)=0;

It needs to take into account this.

"Marios Karaoulis" <marios_hellas_23@yahoo.gr> wrote in message <i9aem4$hnp$1@fred.mathworks.com>...
> > if all(node_ebc~=j)
>
>
> By adding this, now k is not correct with the one I used
> k(node_ebc,:)=0;
> k(:,node_ebc)=0;
>
> It needs to take into account this.

I found it.

tmp=[1:1:num_nodes];
tmp2=setdiff(tmp,node_ebc);
k2(node_ebc,tmp2)=0;
k2(tmp2,node_ebc)=0;

>
> I found it.
>
> tmp=[1:1:num_nodes];
> tmp2=setdiff(tmp,node_ebc);
> k2(node_ebc,tmp2)=0;
> k2(tmp2,node_ebc)=0;



Now this part is a lot slower when I compare it to the original one

k2(node_ebc,:)=0
k2(:,node_ebc)=0;

Any more ideas?

Marios Karaoulis <marios.karaoulis@gmail.com> wrote in message <73f25889-8b76-4aa7-abfe-a6bb1923b88e@t11g2000vbc.googlegroups.com>...
> Hello,
> I need some help on vectorizing the following loop
>
> for t=1:num_probes
> for i=1:num_ebc
> f(node_ebc(i),t)=prop_ebc(i,t);
> k(node_ebc(i),node_ebc(i))=1;
> for j=1:num_nodes
> if node_ebc(i)~=j
> f(j,t)=f(j,t)-k(j,node_ebc(i))*prop_ebc(i,t);
> k(node_ebc(i),j)=0;
> k(j,node_ebc(i))=0;
> end
> end
> end
> end
>

in such scenario, oftentimes a well designed for-loop outperforms the best vectorization.

Hi again,
I came back again after some time. From my research the original code was correct, so I come back a few steps.

%some random data
clear all
num_nodes=1000;
num_probes=15;
num_ebc=50;
k=rand(num_nodes,num_nodes);
f=zeros(num_nodes,num_probes);
prop_ebc=rand(num_ebc,num_probes);
node_ebc=randi([1 num_nodes],num_ebc,1);
node_ebc=unique(node_ebc);
num_ebc=length(node_ebc);
prop_ebc=rand(num_ebc,num_probes);


k2=k;
f2=f;

k3=k;
f3=f;

k4=k;
f4=f;

% End of random data



k2=k; %keep a backup of k, to use it inside loop


% Initial code
tic
for t=1:num_probes
     for i=1:num_ebc
        f2(node_ebc(i),t)=prop_ebc(i,t);
         k3(node_ebc(i),node_ebc(i))=1;
         for j=1:num_nodes
  if (node_ebc(i)==j); continue; end
                          f2(j,t)=f2(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
                          k3(node_ebc(i),j)=0;
                          k3(j,node_ebc(i))=0;
               
                    end
    end
end
toc


%And the matlab code is this

tic
for t=1:num_probes
     for i=1:num_ebc
        f3(node_ebc(i),t)=prop_ebc(i,t);
        jlist = setdiff(1:num_nodes,node_ebc(i));
        f3(jlist,t)= f3(jlist,t)-sum(k2(jlist,node_ebc(i)).*prop_ebc(i,t),2);
   end
end

k4(node_ebc,:)=0;
k4(:,node_ebc)=0;
k4(sub2ind(size(k4),node_ebc,node_ebc))=1;
toc


isequal(f2,f3)
isequal(k4,k3)

The loop is really slow. How can I speed it up?

Or to be more accurate (forget the previous version)





clear all
num_nodes=10;
num_probes=1;
num_ebc=4;


k=rand(num_nodes,num_nodes);
% k=eye(num_nodes);
f=zeros(num_nodes,num_probes);
prop_ebc=rand(num_ebc,num_probes);
node_ebc=randi([1 num_nodes],num_ebc,1);
node_ebc=unique(node_ebc);
num_ebc=length(node_ebc);
prop_ebc=rand(num_ebc,num_probes);


k2=k;
f2=f;

k3=k;
f3=f;

k4=k;
f4=f;
%--------------------First Attempt ----------------------------------------
% % tic
% % for t=1:num_probes
% % for i=1:num_ebc
% % f(node_ebc(i),t)=prop_ebc(i,t);
% % k(node_ebc(i),node_ebc(i))=1;
% % for j=1:num_nodes
% % if (node_ebc(i)~=j)
% % % if all(node_ebc~=j)
% % f(j,t)=f(j,t)-k2(j,node_ebc(i))*prop_ebc(i,t);
% % k(node_ebc(i),j)=0;
% % k(j,node_ebc(i))=0;
% % end
% % end
% % end
% % end
% % toc

%--------------------Second Attempt ---------------------------------------

for t=1:num_probes
    k3=k2;
    for i=1:num_ebc
       f2(node_ebc(i),t)=prop_ebc(i,t);
       k3(node_ebc(i),node_ebc(i))=1;
        for j=1:num_nodes
if (node_ebc(i)==j); continue; end
                  f2(j,t)=f2(j,t)-k3(j,node_ebc(i))*prop_ebc(i,t);
                  k3(node_ebc(i),j)=0;
                  k3(j,node_ebc(i))=0;
              
        end
    end
end










%---------------------Third Attempt ---------------------------------------

tic
% [nn tt] = meshgrid(1:num_probes,node_ebc);
% f4(sub2ind(size(f2),tt,nn)) = prop_ebc;
% % f2(node_ebc,:)=prop_ebc;
% jlist = setdiff(1:num_nodes,node_ebc);

for t=1:num_probes
% k5=k;
    for i=1:num_ebc
% f4(node_ebc(i),t)=prop_ebc(i,t);
% k5(node_ebc(i),node_ebc(i))=1;
       jlist = setdiff(1:num_nodes,node_ebc(i));
       f4(jlist,t)= f4(jlist,t)-sum(k2(jlist,node_ebc(i)).*prop_ebc(i,t),2);
    end
end

f4(node_ebc,:)=prop_ebc;


k4(node_ebc,:)=0;
k4(:,node_ebc)=0;
k4(sub2ind(size(k2),node_ebc,node_ebc))=1;
toc

isequal(k3,k4)
isequal(f2,f4)