matrix vectorization finite difference avoid for loop

hi,

I am one of the winners of the recent simulink student challenge, as you can see the work in coding never concludes ;-)

I am working lately with the vectorized code to improve performance and readability, more precicely I am replacing my 2 nested for loop with a single vectorized matrix, using a finite difference method for solving a PDE.

Unfortunately I have been able to replace only the first loop because I solve the temperature field along rows (vectorized form) and I update it during time (columns) with the for loop that I'd like to avoid.

To simplify my question now it looks like:

for kk = 1:l_t-1
T(2:end,kk+1) = Dt./(1+(m_is*cp)./(rho_oil_f(T(2:end,kk)).*cp_oil_f(T(2:end,kk))*A)).*...
    ((-kx_f(T(2:end,kk),Tinf_e,lambdaV,V_HTF_f,l_x) + 1./Dt - V_HTF_f(T(2:end,kk))./Dx + (m_is*cp)./(rho_oil_f(T(2:end,kk)).*cp_oil_f(T(2:end,kk))*A*Dt)).*T(2:end,kk) +...
    (V_HTF_f(T(2:end,kk))./Dx).*T(1:end-1,kk) + kx_f(T(2:end,kk),Tinf_e,lambdaV,V_HTF_f,l_x)*Tinf_e);
end

and I'd like to have:

T(2:end,2:end) = Dt./(1+(m_is*cp)./(rho_oil_f(T(2:end,1:end-1)).*cp_oil_f(T(2:end,1:end-1))*A)).*...
    ((-kx_f(T(2:end,1:end-1),Tinf_e,lambdaV,V_HTF_f,l_x) + 1./Dt - V_HTF_f(T(2:end,1:end-1))./Dx + (m_is*cp)./(rho_oil_f(T(2:end,1:end-1)).*cp_oil_f(T(2:end,1:end-1))*A*Dt)).*T(2:end,1:end-1) +...
    (V_HTF_f(T(2:end,1:end-1))./Dx).*T(1:end-1,1:end-1) + kx_f(T(2:end,1:end-1),Tinf_e,lambdaV,V_HTF_f,l_x)*Tinf_e);

the solution is obviously incorrect, thus I have tried to reproduce the computation with two simple schemes. For example it's clear how it works when I write:

clc;clear;
dati = ones(5);
dati(:,2) = dati(:,1) + dati(:,2)*2;
dati(:,3) = dati(:,2) + dati(:,3)*2;
dati(:,4) = dati(:,3) + dati(:,4)*2;
dati(:,5) = dati(:,4) + dati(:,5)*2;
dati

But the same cannot be said for this other example which I would like to implement:

dati2 = ones(5);
dati2(:,2:end) = dati2(:,1:end-1) + dati2(:,2:end)*2

I cannot understand how can I obtain the same procedure of the previous example in the complete vectorized form.

many thanks !

Roberto

[Joel Handy]

dati2 = ones(5);
dati2(:,2:end) = dati2(:,1:end-1) + dati2(:,2:end)*2

This code is doing elementwise operations. . .mostly. Technically dati2(:,2:end)*2 is a matrix multiplication, but it looks like you know the difference between the '*' and '.*' operators based on the code above.

You seem to have a pretty good grasp of element wise operations baed on your first two code examples. Your last two code examples are directly equivalent. The latter example just has more compact indexing notation. In both, You are taking a 4 element sub array out of dati2 and adding it elementwise to two times a different 4 element sub array of dati2. A third way you could write this examples is:

dati3 = ones(5);
subA = dati2(:,1:end-1);
subB = dati2(:,2:end);
newSub = subA + subB*2;
dati2(:,2:end) = newSub

Does that make any sense? I'm a little worried I have missed the mark on your question since it seems like you've already vectorized am uch more complicated equation.

[Joel Handy]

Rereading your question and looking back at your first two examples, it looks like each column of T is dependant on the values in the previous column. I dont think you will be able to vectorize this.

Matlab's original purpose is to perform matrix calculations. the underlying code is optimzed for these sorts of cacluations. Therefore, if you can formulate your problem as a matrix operation, it will run more efficiently. Vectorization is that process of formulating a problem into a matrix calculation.

Your calculations appear to be, by necessity, sequential dictating a for loop.

[Roberto Tascioni]

thank you Joel Handy

you are right, this process is intrinsically sequential and unfortunately it cannot benefit by vectorization.

Answer 1

[Pravin Jagtap]

From the code, I am assuming

• The matrix used for computation represents the ‘Temperature’ along row-wise direction and ‘Time’ along column-wise direction.
• Explicit Finite Difference solver.

In order to vectorize the matrix operation, one must make sure the correctness of output. The sample code of Row-wise vectorized v/s Row and Column-wise vectorized version (given by you) gives different results.

Row-wise vectorized version:

clc;clear;

dati = ones(5);

dati(:,2) = dati(:,1) + dati(:,2)*2;

dati(:,3) = dati(:,2) + dati(:,3)*2;

dati(:,4) = dati(:,3) + dati(:,4)*2;

dati(:,5) = dati(:,4) + dati(:,5)*2;

dati;

In the above code, you are updating matrix column-wise (in time-direction) and all column operations are vectorized at every ‘Time’-step. Important thing here is, update of each column depends on values of previous column and current column.

Row and Column-wise vectorized version:

dati2 = ones(5);

dati2(:,2:end) = dati2(:,1:end-1) + dati2(:,2:end)*2;

In the above code, you are intending to do vectorization in both directions. The second line of above code is doing element-wise matrix addition which vectorized (but performing different operations compared with the first version).

(Note: Check the output of ‘dati’ and ‘dati2’ and verify the correctness before going for optimized code.)

Finite-Difference Code:

The code given by you cannot be vectorized in both the direction (‘Temperature’ and ‘Time’) because of data dependencies (updating ‘Temperature’ columns is iterative process in ‘time’ direction). Considering the explicit solver of finite-difference you can certainly vectorize the operation of updating all elements of column (‘Temperature’) at each timestep.

For more details on Vectorization:

https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html?s_tid=srchtitle