Prev: cell, indicator, replacing
Next: MZDDE toolbox question
From: Vladimir on 3 Aug 2010 13:43
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:

grid(:,:)=(-alpha.*(grid(:,:)+grid(:,:))-gamma.*(grid(:,:)+grid(:,:))-2*beta.*x_dd_zeta_etha)./(-2*(alpha+gamma));

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 Gauss-Seidel 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.
Help?
From: Roger Stafford on 3 Aug 2010 15:17
"Vladimir " <something(a)or.another> wrote in message <i39kf9$bt5$1(a)fred.mathworks.com>...
> 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:
>
> grid(:,:)=(-alpha.*(grid(:,:)+grid(:,:))-gamma.*(grid(:,:)+grid(:,:))-2*beta.*x_dd_zeta_etha)./(-2*(alpha+gamma));
>
> 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 Gauss-Seidel 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.
> Help?
- - - - - - - - - - -
In the "vectorized" form you have written, it doesn't look like Gauss-Seidel 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.

Note that this is general different from the results you would get if you put it in nested for-loops like this:

for i1 = 1:n
for i2 = 1:n
grid(i1,i2) = the above expression in grid(i1,i2) ;
end
end

A lot of beginners in matlab make the mistake of not distinguishing between these two different ways of proceeding.

Roger Stafford
From: Vladimir on 3 Aug 2010 15:58
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <i39pvi$7g8$1(a)fred.mathworks.com>...
> - - - - - - - - - - -
> In the "vectorized" form you have written, it doesn't look like Gauss-Seidel 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.
>
> Note that this is general different from the results you would get if you put it in nested for-loops like this:
>
> for i1 = 1:n
> for i2 = 1:n
> grid(i1,i2) = the above expression in grid(i1,i2) ;
> end
> end
>
> A lot of beginners in matlab make the mistake of not distinguishing between these two different ways of proceeding.
>
> Roger Stafford

Roger, thank you.
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?
In general though, using .* and so on is substantially faster then using a loop for the same calculation?
From: Roger Stafford on 3 Aug 2010 19:16
"Vladimir " <something(a)or.another> wrote in message <i39scc$h7k$1(a)fred.mathworks.com>...
> Roger, thank you.
> 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?
> In general though, using .* and so on is substantially faster then using a loop for the same calculation?
- - - - - - - - - - -
Yes, if you want some kind of element-by-element 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 "element-by-element updating" I mean having some elements within the single operation updated and using those updated results in other later parts of the same operation.

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.

The same cannot be said for matrix multiplication, *, (without the 'dot'.) If you wrote

X = A*X;

where A is an n by n matrix and X is an n by 1 column vector, the results would be different from writing

for i = 1:n
X(i) = A(i,:)*X;
end

though

Y = zeros(n,1);
for i = 1:n
Y(i) = A(i,:)*X;
end
X = Y;

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.

Roger Stafford
 | 
Pages: 1
Prev: cell, indicator, replacing
Next: MZDDE toolbox question