http://www.mathworks.com/matlabcentral/newsreader/view_thread/288428
MATLAB Central Newsreader  Order of operations in a vectorized function
Feed for thread: Order of operations in a vectorized function
enus
©19942015 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Tue, 03 Aug 2010 17:43:05 +0000
Order of operations in a vectorized function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288428#768108
Vladimir
I was wondering if there is any documentation on the internal workings of vectorized operations. Specifically, I have programmed an iterative solver for a system of PDEs, and the next iterations values are found using a line much like this one:<br>
<br>
grid(:,:)=(alpha.*(grid(:,:)+grid(:,:))gamma.*(grid(:,:)+grid(:,:))2*beta.*x_dd_zeta_etha)./(2*(alpha+gamma));<br>
<br>
all are matrices of the same size. I tried running it both like this, which I hoped would be a good replacement for running looped GaussSeidel iteration, and in Jacobi form  when the new values for grid are calculated based specifically on the previous iterations values. The convergence speed was different, but strangely enough, it was faster for the Jacobi version. Then I realized that for my code to really be GS iteration, it needs to calculate each value using the newest values available, but I have no way of knowing if that is the case, because .* and .^'s aren't transparent to the user. <br>
Help?

Tue, 03 Aug 2010 19:17:06 +0000
Re: Order of operations in a vectorized function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288428#768141
Roger Stafford
"Vladimir " <something@or.another> wrote in message <i39kf9$bt5$1@fred.mathworks.com>...<br>
> I was wondering if there is any documentation on the internal workings of vectorized operations. Specifically, I have programmed an iterative solver for a system of PDEs, and the next iterations values are found using a line much like this one:<br>
> <br>
> grid(:,:)=(alpha.*(grid(:,:)+grid(:,:))gamma.*(grid(:,:)+grid(:,:))2*beta.*x_dd_zeta_etha)./(2*(alpha+gamma));<br>
> <br>
> all are matrices of the same size. I tried running it both like this, which I hoped would be a good replacement for running looped GaussSeidel iteration, and in Jacobi form  when the new values for grid are calculated based specifically on the previous iterations values. The convergence speed was different, but strangely enough, it was faster for the Jacobi version. Then I realized that for my code to really be GS iteration, it needs to calculate each value using the newest values available, but I have no way of knowing if that is the case, because .* and .^'s aren't transparent to the user. <br>
> Help? <br>
          <br>
In the "vectorized" form you have written, it doesn't look like GaussSeidel iteration to me. For each element of grid, matlab will use the old values in grid for all its computations. It accomplishes this by temporarily storing two copies of the grid matrix in separate buffers, the old values and the new values. The computation on the right takes from the old value matrix and puts the results into the new value matrix. Imagine that you wrote "temp_grid" on the left side and after the computation was finished for all elements you then wrote grid = temp_grid to copy the new results back into grid. That would be equivalent to the vectorized action.<br>
<br>
Note that this is general different from the results you would get if you put it in nested forloops like this:<br>
<br>
for i1 = 1:n<br>
for i2 = 1:n<br>
grid(i1,i2) = the above expression in grid(i1,i2) ;<br>
end<br>
end<br>
<br>
A lot of beginners in matlab make the mistake of not distinguishing between these two different ways of proceeding.<br>
<br>
Roger Stafford

Tue, 03 Aug 2010 19:58:04 +0000
Re: Order of operations in a vectorized function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288428#768156
Vladimir
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i39pvi$7g8$1@fred.mathworks.com>...<br>
>           <br>
> In the "vectorized" form you have written, it doesn't look like GaussSeidel iteration to me. For each element of grid, matlab will use the old values in grid for all its computations. It accomplishes this by temporarily storing two copies of the grid matrix in separate buffers, the old values and the new values. The computation on the right takes from the old value matrix and puts the results into the new value matrix. Imagine that you wrote "temp_grid" on the left side and after the computation was finished for all elements you then wrote grid = temp_grid to copy the new results back into grid. That would be equivalent to the vectorized action.<br>
> <br>
> Note that this is general different from the results you would get if you put it in nested forloops like this:<br>
> <br>
> for i1 = 1:n<br>
> for i2 = 1:n<br>
> grid(i1,i2) = the above expression in grid(i1,i2) ;<br>
> end<br>
> end<br>
> <br>
> A lot of beginners in matlab make the mistake of not distinguishing between these two different ways of proceeding.<br>
> <br>
> Roger Stafford<br>
<br>
Roger, thank you. <br>
So, it seems there's no way of avoiding using loops in this case? "Vectorization" is only good when all the data are on the same 'time' level in the iteration, and if I specifically want to use the updated values, I need to loop? <br>
In general though, using .* and so on is substantially faster then using a loop for the same calculation?

Tue, 03 Aug 2010 23:16:42 +0000
Re: Order of operations in a vectorized function
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288428#768226
Roger Stafford
"Vladimir " <something@or.another> wrote in message <i39scc$h7k$1@fred.mathworks.com>...<br>
> Roger, thank you. <br>
> So, it seems there's no way of avoiding using loops in this case? "Vectorization" is only good when all the data are on the same 'time' level in the iteration, and if I specifically want to use the updated values, I need to loop? <br>
> In general though, using .* and so on is substantially faster then using a loop for the same calculation?<br>
          <br>
Yes, if you want some kind of elementbyelement updating, you would not want to vectorize it. Moreover it would have to be just the right kind of loops to match the sequencing you desire, whatever that is. By "elementbyelement updating" I mean having some elements within the single operation updated and using those updated results in other later parts of the same operation.<br>
<br>
As for the use of .* as opposed to multiplying elements one at a time, there should be no question of updating a portion of the results and thereby affecting later results, so there is no particular reason to do it with loops. That is because each multiplication is between corresponding elements and there is no interaction between different pairs in such an operation. The products of n different pairs of factors gives n different answers which are independent of each other.<br>
<br>
The same cannot be said for matrix multiplication, *, (without the 'dot'.) If you wrote<br>
<br>
X = A*X;<br>
<br>
where A is an n by n matrix and X is an n by 1 column vector, the results would be different from writing<br>
<br>
for i = 1:n<br>
X(i) = A(i,:)*X;<br>
end<br>
<br>
though<br>
<br>
Y = zeros(n,1);<br>
for i = 1:n<br>
Y(i) = A(i,:)*X;<br>
end<br>
X = Y;<br>
<br>
would give the same results as the vectorized operation. Try it on your computer and see. That is because the results of the changes to earlier X(i)'s would affect the computation of the later ones, whereas this wouldn't happen in the vectorized version. The one you would use in such a case depends on what kind of algorithm you wish to carry out.<br>
<br>
Roger Stafford