Mark.F
Mark.F

Reputation: 1694

How can I vectorize the following loops in matlab?

I am using an R2019a matlab. I have a matrix that I want to fill with values at certain positions, based on a calculation done on 2 additional vectors. Currently I'm doing it with 2 loop, but this doesn't take advantage of matlab's vectorization abilities.

How can the following script be performed in a vectorized manner:

C = zeros(size(vecB), size(vecA));

% Calculate face-vertix connectivity of face area values:
posMat = sparse(A, repmat((1:size(A,1))',1,3), ...
                1, size(B,1), size(A,1));

for i = 1:size(B,1)
    for j = 1:size(A,1)
        if posMat(i,j) == 1
            C(i,j) = vecA(j)/vecB(i)/3;
        end
    end
end

Sizes of variables in the script:

size(A)  =         5120           3
size(B)  =         2562           3
size(B)  =         2562        5120
size(posMat)  =    2562        5120
size(vecA) =       5120           1
size(vecB) =       2562           1

Upvotes: 0

Views: 93

Answers (1)

Will
Will

Reputation: 1880

Because the condition inside your loop is only satisfied by nonzero elements of posMat, the important thing for efficiency is to take advantage of the fact that this matrix is sparse in a vectorized implementation.

If you use multiplication by a sparse condition as the means of setting the elements that fail the condition to 0 (rather than initializing the matrix using zeros) then only the elements that pass the condition will actually be evaluated.

From R2016b onwards, vecA.'./vecB/3 is implicitly expanded to return the full 2562×5120 double array of vecA(j)/vecB(i)/3 for all values of i and j. But (posMat == 1) .* vecA.'./vecB/3 returns a 2562×5120 sparse double array where the division has only been evaluated for the elements where posMat == 1.

If a sparse value of C if acceptable, then

C = (posMat == 1) .* vecA.'./vecB/3;

will suffice. If the full-storage form is required this output can simply be passed to the full function.

Upvotes: 3

Related Questions