# How to correctly vectorize?

1 view (last 30 days)

Show older comments

Hi,

I have this code (Matlab 2018b):

m=5;

f=[.1,-3,.05,70.56,110.32456];

c=zeros(m);

M=10;

for i=1:M

for j=1:m

fj=f(j);

for k=j:m

c(j,k)=c(j,k)+fj*f(k);

c(k,j)=c(j,k);

end

end

end

cc=M*(f'*f);

c-cc

If M=1 (or =2) the result is all zeros. If M=10, the result is not all zeros, but some. If M=100, the result is not zeros at all. I have plenty of this type of code and want to accelerate with vectorization, but I am confused about the results.

What is the correct vectorization of these kind of for loops? Why it is not zero all the times? I migt imagine that the result is around the minimum number of representation but here the difference is 1e-12 - 1e-17. It seems to me way too high.

So what should I do? Which is correct, vectorized or for loop? With for loops it works correctly.

Csaba

##### 4 Comments

Jan
on 20 Dec 2018

### Accepted Answer

Jan
on 20 Dec 2018

Edited: Jan
on 20 Dec 2018

The differences between the loop and the linear algebra implementation have the expected range. The matrix multiplication uses highly optimized BLAS routines. The dot products contain a sum and summing is numerical instable in general, see e.g. https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum

If you display the relative errors, you see that the deviations are in the magnitude of eps:

(c - cc) ./ c

This is the typical dimension of errors, e.g. caused by calculating the values one time with floating point commands and the other time with SSE/AVX. Both results are correct.

What do you consider as correct result of:

1 + 1e-17 - 1

##### 5 Comments

Jan
on 20 Dec 2018

Edited: Jan
on 21 Dec 2018

By the way, eps is 2.2e-16.

Maybe this thread clarifies some details: https://www.mathworks.com/matlabcentral/answers/57444-faq-why-is-0-3-0-2-0-1-not-equal-to-zero

0.3 - 0.2 - 0.1

This is not 0.0 also.

There are several methods for creating a sum with reduced rounding effects: https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!