Eliminate for loop: can this code be vectorized?

Hi all,

I am trying to optimize my code, my question is if it is possible to vectorize the following code:

for r = 1:2:(N-1)
  U(r,1)=U(r,1)+ (1./K(r,r)).*(F(r,1)-K(r,:)*U(:,1));
end

where U and F are Nx1 vectors and K is a NxN matrix. This part of the code takes a considerable amount of time when N becomes large. Maybe someone has a suggestion on how to vectorize this part?

Thank you!

Answer 1

[KSSV]

r = 1:2:(N-1) ;
U(r,1)=U(r,1)+ (1./K(r,r))*(F(r,1)-K(r,:)*U(:,1));

unchecked check the values before using.

[Quinten Rensen]

Thank you for your answer, that seems to be a nice solution!

Unfortunately I get the following error message:

"Unable to perform assignment because the size of the left side is 289-by-1 and the size of the right side is 289-by-289."

while N = 8450. I think it goes wrong for K(r,r). I will have a look at it.

[Jan]

Yes, K(r,r) is the problem. It is a matrix, but you need a vector here. Try:

K(sub2ind(size(K), r, r))

[Quinten Rensen]

That indeed seems a nice solution, however, I get the same error message. This might be due to the K(r,:) term, but when I try to use sub2in() for that term as well, it gives:

""sub2ind" previously appeared to be used as a function or command, conflicting with its use here as the name of a variable. A possible cause of this error is that you forgot to initialize the variable, or you have initialized it implicitly using load or eval."

[Jan]

Please post the piece of code, which reproduces the error. It seems like you have created a variable called "sub2ind".

[Quinten Rensen]

Hi Jan, thank you for your quick replies. Maybe it is useful to note that I have defined K as a sparse matrix, so K = [(indices) , value ]

Here is the original code:

function [U,res]=GS(U,F,nelx,nely,sweeps,K)

omega = 1;
count =1;
N = 2.*(nelx+1).*(nely+1);
while count <= sweeps      %Number of sweeps;

  for r = 1:2:(N-1)
      U(r,1)=U(r,1)+omega.*(1./K(r,r)).*(F(r,1)-K(r,:)*U(:,1));
  end

  for c = 2:2:N
         U(c,1)=U(c,1)+omega.*(1./K(c,c)).*(F(c,1)-K(c,:)*U(:,1));
  end

     count = count + 1;
  end
  %%Calculate final residual at all points
      res = zeros(2.*(nelx+1).*(nely+1),1); %Preallocate
      res = F - K*U;                %All residuals;
  end%End of GS function

And this is the code after I implemented the suggestions:

function [U,res]=GS(U,F,nelx,nely,sweeps,K)

omega = 1;
count =1;
N = 2.*(nelx+1).*(nely+1);
while count <= sweeps      %Number of sweeps;
  r = 1:2:(N-1)
    U(r,1)=U(r,1)+omega.*(1./K(sub2ind(size(K), r, r))).*(F(r,1)-K(sub2ind(size(K), r, :))*U(:,1));
  
  c = 2:2:N
       U(c,1)=U(c,1)+omega.*(1./K(sub2ind(size(K), c, c))).*(F(c,1)-K(sub2ind(size(K), c, :))*U(:,1));
  
     count = count + 1;
  end
  %%Calculate final residual at all points
      res = zeros(2.*(nelx+1).*(nely+1),1); %Preallocate
      res = F - K*U;                %All residuals;
  end%End of GS function

[Quinten Rensen]

Ok, I found what was wrong, I had to use the transpose of (1./K(sub2ind(size(K), r, r))). Then the code becomes

r = 1:2:(N-1)
    U(r,1)=U(r,1)+omega.*(1./K(sub2ind(size(K), r, r)))'.*(F(r,1)-K(r,:)*U(:,1));

Thank you for your help!