Get rid of for loop

Hello,

I need to get rid of loops in my Matlab code. During this procedure, i encountered the following problem. Here is the part of the code

Nx=11; Ny=11; 
A=rand(Ny,Nx);
A1=A;
 
       for j=2:Ny-1      
       for i=2:Nx-1
          A(j,i)=A(j,i-1);
       end;
       end;
 
A1(2:Ny-1,2:Nx-1)=A1(2:Ny-1,1:Nx-2);
 
A2=A-A1;

I get that A is not equal to A1. The situation is the same with the j-1 index on the right side. I found out that in the case of

A1(2:Ny-1,2:Nx-1)=A1(2:Ny-1,1:Nx-2)

Matlab does only 1 operation instead of filling the entire matrix. However, if we take the index i + 1 or j + 1, then everything works fine. Is there any way to fix it without using the loops?

[Steven Lord]

I need to get rid of loops in my Matlab code.

Why? If the answer is "because loops are slow" sigh that has been the case at some times in the past, but MATLAB has improved quite a bit over the past two decades since I first heard that complaint.

[Al Ne]

Because I do parallel GPU computing.

Answer 1

[Walter Roberson]

Consider the difference between

for j=2:Ny-1      
   for i=2:Nx-1
      A(j,i) = A(j,i-1);
   end
end

and

old_A = A;
for j=2:Ny-1      
   for i=2:Nx-1
      A(j,i) = old_A(j,i-1);
   end
end

In the second of those, it is clear that none of the assignments into A can possibly affect the data that is to be copied. But in the first one, for each row, the first column, A(j,1) is copied into A(j,2), and then the next i iteration, A(j,2) is copied into A(j,3) -- but the A(j,2) is the one that was just copied from A(j,1). So the effect of your code is to copy the first entry of each row into all the other columns of the same row, except for the last column (because you stop at Nx-1 )

When you do a vectorized copy such as

A1(2:Ny-1,2:Nx-1)=A1(2:Ny-1,1:Nx-2);

then MATLAB is internally doing the equivalent of copying to a tempory array and doing the assignments out of that copy. Values are not propogated across. The internals are

temporary = A1(2:Ny-1,1:Nx-2);
A1(2:Ny-1,2:Nx-1) = temporary;

-- the act of indexing causes a temporary array with the selected data to be created, and data is copied out of that temporary into the destination locations.

If you wanted to copy the first element of each column across to all columns except the last, then

A(:,2:end-1) = A(:,1).*ones(1,size(A,2)-2);

[Al Ne]

Dear Walter,

Thank you so much for your clear explanation! My task is not just copping the first element. I need to compute the following relation

A= A(2:Ny-1,3:Nx)-A(2:Ny-1,1:Nx-2)

the first term in the right side works fine whereas the second one is not for the reason you described above. As you explained, it can be fixed if i create an additional array like

An=A;

A= An(2:Ny-1,3:Nx)-An(2:Ny-1,1:Nx-2)

However, this variant is not appropriate in my case due to the specific of the code. Is there any way to fix the relation without indtroducing the additional array?

[Walter Roberson]

A= A(2:Ny-1,3:Nx)-A(2:Ny-1,1:Nx-2)

In MATLAB, that is not equivalent to the for loop you had, for the reasons I described. It is not a matter of which of the two terms "works fine": MATLAB handles the expression as

temporary1 = Ny-1
temporary2 = 2:temporary1
(release temporary1->data, release temporary1->attributes, no actual name to remove)
temporary3 = 3:Nx
temporary4 = A(temporary2, temporary3)
(release temporary2->data, release temporary2->attributes, no actual name to remove)
(release temporary3->data, release temporary3->attributes, no actual name to remove)
temporary5 = Ny-1
temporary6 = 2:temporary1
(release temporary5->data, release temporary3->attributes, no actual name to remove)
temporary7 = Nx-2
temporary8 = 1:temporary7
(release temporary7->data, release temporary7->attributes, no actual name to remove)
temporary9 = A(temporary6, temporary8)
(release temporary6->data, release temporary6->attributes, no actual name to remove)
(release temporary8->data, release temporary8->attributes, no actual name to remove)
temporary10 = temporary4 - temporary5
(release temporary4->data, release temporary4->attributes, no actual name to remove)
(release temporary5->data, release temporary5->attributes, no actual name to remove)
(release A->data, release A->attributes;
 set A->attributes <- temporary10->attributes including data pointer)

In the above, = is short form that includes computation and creation of data attributes

There is no sub-expression there that is done "wrong": both are handled the same way, involving lots of generation of unnamed variables.

In MATLAB, the difference between an unnamed temporary variable and a named variable, is that a named variable is in the name table and "points" to what used to be the attribute block of the temporary -- the same pointer that would have been passed around before gets recorded in the name table. Conceptually the "global" and "persistent' attribute reside in the name table, but in reality the flags are part of the attributes block, and where I wrote

 (release A->data, release A->attributes;
 set A->attributes <- temporary10->attributes including data pointer)

there is really another step

 (set temporary10->attributes->flags <- A->attributes->flags;
  release A->data, release A->attributes;
  set A->attributes <- temporary10->attributes including data pointer)

so that global / persistent status get copied from the destination variable to the new space.

If this is not what you want to have happen, if you want A(2:Ny-1,3:Nx)-A(2:Ny-1,1:Nx-2) to be implemented in a way that propogates just-changed values, then I am going to need a reference implementation as for loops, and an explicit acknowledgement that you want the just-updated location to be the source for the next iteration's calculation.

[Al Ne]

Thank you so much for the clear explanation!