Fastest way to implement the given lines of code.

Let A, B, C be any random matrices of order mxn; m,n would be large enough. What is the fastest way to implement the following lines:

A = rand(m,n) ;
B = rand(m,n) ; 
C = rand(m,n) ; 
t1 = tic ;
for j = 2:n
    for i = 2:m
        H = 0.5*(A(i,j) + A(i+1,j)) ;
        if (H > 0.0)
            C(i,j) = C(i,j)-H*(B(i+1,j)-B(i,j)) ;
        end 
    end 
end 
toc(t1)

cheers...

Sreenivas

[James Tursa]

As currently written, I think the fastest way is this:

error('Index exceeds matrix dimensions');

because when i=m, A(i+1,j) will generate an error. Similar problem with B(i+1,j).

Did you mean i=1:m-1? Please double check your code and let us know what the real algorithm is.

[KSSV]

Dear James Tursa

Yes you are right...i = 1:m-1. I apologize for the mistake in the loop.

Answer 1

[Matt J]

I suspect you're looking for something like the following,

      A = rand(m,n) ;
      B = rand(m,n) ; 
      C = rand(m,n) ; 

   H=0.5*conv2(A,[1;1],'same');    %EDITED
   deltaB=conv2(B,[1;-1],'same');

   idx=H>0;
   C(idx)=C(idx)-H(idx).*deltaB(idx);

[James Tursa]

Missing the 0.5 factor on H

[Matt J]

Fixed!

[Matt J]

As James Tursa said, code 1. is faster compared to other three in the present context of matrices.

I don't see where James said this. In any case, here is a modification to my original proposal. I find it to be 4 times as fast as the originally posted for-loop implementation, when testing with m=n=5000.

       H=conv2(A,0.5*[1;1],'same');   
       deltaB=conv2(B,[1;-1],'same');

       H(H<=0)=0;

       C=C-H.*deltaB;

Answer 2

[KSSV]

Dear James Tursa and Matt

Thanks a lot for the reply. I have tried four different methods: 1. The way the code given above 2. Vectorising the code 3. Code suggested by Matt before editing 4. code suggested by Matt after editing.

As James Tursa said, code 1. is faster compared to other three in the present context of matrices.

Thanking you..

Sreenivas