Assigning to multiple sub-matrices of a matrix simultaneously. Possible optimization by vectorized indices

Question

Is there a clever way to vectorize a for-loop which assigns elements to submatrices of a matrix?
Initially, I had two for-loops:

U=zeros(6*(M-2),M-2);
for k=2:M-3  
    i=(k-1)*6+1; 
    for j=2:M-3
        U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
    end
end

Then I vectorized the inner loop, such that the code now reads

U=zeros(6*(M-2),M-2);
j=2:M-2;
for k=2:M-3
    i=(k-1)*6+1;
    U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
end

This has reduced my CPU time by more than 90%, so I wondered if I could do the same with the outer loop, but it seems a bit tricky, since I assign to (6x1)-matrices within the U matrix. I tried

U=zeros(6*(M-2),M-2);
k=2:M-3;
i=(k-1)*6+1;
j=2:M-2;
U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

but this fails, since i:i+5 only takes out the first 6 indices I want.

I have also tried to use the reshape() function to convert the matrix into a vector, but it still seems difficult to assign to several blocks of elements at once. There are in total three such for-loops in the code, so I guess an alternative optimization is to parallelize them somehow. However, without access to the parallel toolbox, vectorization seems to me as a good solution if possible.

The code is part of a subroutine in a numerical finite difference method for solving a system of 6 equations on a grid, so this question could be relevant for anyone working with matrix calculations on systems of equations, particularly PDEs. Suggestions to optimizing the code would be greatly appreciated!

Answer 1

To understand how you can write the assignment in one line without loops, it may help to draw the array temp as a rectangle. Then, the different summands that will combine to U are nothing but sub-rectangles of temp (or sub-grids, if you want to keep track of individual elements in temp that will result in a specific element of U) that are shifted to the left, right, top, bottom, respectively.

%# define row, column shifts
rowShift = 6;
colShift = 1;

%# That's how we'd like to shift 
%# U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+
%# D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

%# assign U
U = A * temp(rowShift+1 : end-rowShift, colShift+1 : end-colShift) +... 
    B * temp(rowShift+1 : end-rowShift, 1 : end-2*colShift) + ...
    C * temp(rowShift+1 : end-rowShift, 2*colShift+1 : end) + ...
    D * temp(1 : end-2*rowShift, colShift+1 : end-colShift) + ...
    E * temp(2*rowShift+1 : end, colShift+1 : end-colShift);

Answer 2

In order to select a non-rectangular portion of a matrix, you need to use linear indices: in a 3x3 matrix A, A(3,3)==A(9), and A([1 3 5 7 9]) is a vector that can't be achieved through the row/column indexing method.

The sub2ind function converts row/column indices to linear indices, so you can use it in the form sub2ind(size(U),i:i+5,j) to get the linear indices of one block of U. Change your loop to do only the work of collecting the linear indices, and then you can say outside the loop:

U(ind_U) = A*temp(ind_A) + B*temp(ind_B) ...

Also, any time you are are dealing with FDM or FEM, consider whether you should be using sparse matrices.