How can I exploit parallelism when defining values in a sparse matrix?

Question

The following MATLAB code loops through all elements of a matrix with size 2IJ x 2IJ.

    for i=1:(I-2)
        for j=1:(J-2)
            ij1 = i*J+j+1; % row
            ij2 = i*J+j+1 + I*J; % col 
            D1(ij1,ij1) = 2;
            D1(ij1,ij2) = -1;
        end
    end

Is there any way I can parallelize it use MATLAB's parfor command? You can assume any element not defined is 0. So this matrix ends up being sparse (mostly 0s).

What are the values of I and J? How long does this take to run? You can likely easily vectorize this, but I don’t think you can use parfor here. — Cris Luengo
If this takes so long that you need to improve the speed, you don't want to be creating a full matrix in the first place (or probably even a matrix at all). To efficiently make large sparse matrices like this use spdiags. — David
@CrisLuengo - I and J can be a couple hundred, anywhere between 200-300. How do I vectorize the above code and use spdiags to make it more efficient? — user1068636

rahnema1 rahnema1 · Accepted Answer · 2020-10-26T04:36:01

Before using parfor it is recommended to read the guidelines related to decide when to use parfor. Specially this:

Generally, if you want to make code run faster, first try to vectorize it.

Here vectorization can be used effectively to compute indices of the nonzero elements. Those indices are used in function sparse. For it you need to define one of i or j to be a column vector and another a row vector. Implicit expansion takes effect and indices are computed.

I = 300;
J = 300;
i = (1:I-2).';
j = 1:J-2;
ij1 = i*J+j+1;
ij2 = i*J+j+1 + I*J;

D1 = sparse(ij1, ij1, 2, 2*I*J, 2*I*J) +  sparse(ij1, ij2, -1, 2*I*J, 2*I*J);

However for the comparison this can be a way of using parfor (not tested):

D1 = sparse (2*I*J, 2*I*J);
parfor i=1:(I-2)
    for j=1:(J-2)
        ij1 = i*J+j+1; 
        ij2 = i*J+j+1 + I*J;  
        D1 = D1 + sparse([ij1;ij1], [ij1;ij2], [2;-1], 2*I*J, 2*I*J) ;
    end
end

Here D1 used as reduction variable.

How can I exploit parallelism when defining values in a sparse matrix?

1 Answers