MATLAB: How to make this function run faster

tensor

Hello,
I have this function inside a scheme which takes 70% of the total time of running the script. The function decomposes a 4th order tensor such that
% INPUT 
% F = 2x2 matrix; D = 3x3 matrix;
%








% OUTPUT
F = eff;
 % =========================== Calculate D_UL
J = det(F);
% ================================================================
% global i j k l
%%%   ========= C_ijkl = J^-1 * F_im * F_jn * D_mnpq * F_kp * F_lq ========
D1111 = D(1,1); D1112 = D(1,3); D1121 = D(3,1); D1122 = D(1,2);
D1211 = D(3,1); D1212 = D(3,3); D1221 = D(3,3); D1222 = D(3,2);
D2221 = D(2,3); D2222 = D(2,2); D2121 = D(3,3); D2122 = D1222 ;
D2211 = D(2,1); D2212 = D(2,3); D2111 = D(3,1); D2112 = D(3,3);
% ----------------------------

for i = 1:2
    for j = 1:2
        for k = 1:2
            for l = 1:2
%
AA = inv(J) * F(i,1) * F(j,i) * D1111 * F(k,1) * F(l,1) + inv(J) * F(i,1) * F(j,1) * D1112 * F(k,1) * F(l,2);
%
BB = inv(J) * F(i,1) * F(j,i) * D1121 * F(k,2) * F(l,1) + inv(J) * F(i,1) * F(j,i) * D1122 * F(k,2) * F(l,2);
%
CC = inv(J) * F(i,1) * F(j,2) * D1211 * F(k,1) * F(l,1) + inv(J) * F(i,1) * F(j,2) * D1212 * F(k,1) * F(l,2);
%
DD = inv(J) * F(i,1) * F(j,2) * D1221 * F(k,2) * F(l,1) + inv(J) * F(i,1) * F(j,2) * D1222 * F(k,2) * F(l,2);
%
EE = inv(J) * F(i,2) * F(j,1) * D2111 * F(k,1) * F(l,1) + inv(J) * F(i,2) * F(j,1) * D2112 * F(k,1) * F(l,2);
%
FF = inv(J) * F(i,2) * F(j,1) * D2121 * F(k,2) * F(l,1) + inv(J) * F(i,2) * F(j,1) * D2122 * F(k,2) * F(l,2);
%
GG = inv(J) * F(i,2) * F(j,2) * D2211 * F(k,1) * F(l,1) + inv(J) * F(i,2) * F(j,2) * D2212 * F(k,1) * F(l,2);
%
HH = inv(J) * F(i,2) * F(j,2) * D2221 * F(k,2) * F(l,1) + inv(J) * F(i,2) * F(j,2) * D2222 * F(k,2) * F(l,2);
% ============
C(i,j,k,l)=C(i,j,k,l)+ AA+BB+CC+DD+EE+FF+GG+HH;
% ----------------------------
            end
        end
    end
end
D_UL = [C(1,1,1,1) C(1,1,2,2) C(1,1,1,2);
        C(2,2,1,1) C(2,2,2,2) C(2,2,1,2);
        C(1,2,1,1) C(1,2,2,2) C(1,2,1,2)];

Best Answer

For starters, try...
...
IJ=inv(J);
for i = 1:2
  for j = 1:2
    for k = 1:2
      for l = 1:2
        AA = IJ*F(i,1)*F(j,i)*F(k,1)*[D1111*F(l,1) + D1112*F(l,2)];
        BB = IJ*F(i,1)*F(j,i)*F(k,2)*[D1121*F(l,1) + D1122*F(l,2)];
        CC = IJ*F(i,1)*F(k,1)*F(j,2)*[D1211*F(l,1) + D1212*F(l,2)];
        DD = IJ*F(i,1)*F(j,2)*F(k,2)*[D1221*F(l,1) + D1222*F(l,2)];
        EE = IJ*F(i,2)*F(j,1)*F(k,1)*[D2111*F(l,1) + D2112*F(l,2)];
        FF = IJ*F(i,2)*F(j,1)*F(k,2)*[D2121*F(l,1) + D2122*F(l,2)];
        GG = IJ*F(i,2)*F(j,2)*F(k,1)*[D2211*F(l,1) + D2212*F(l,2)];
        HH = IJ*F(i,2)*F(j,2)*F(k,2)*[D2221*F(l,1) + D2222*F(l,2)];
        C(i,j,k,l)=C(i,j,k,l)+ AA+BB+CC+DD+EE+FF+GG+HH;
      end
    end
  end
end
...
The next step is to move the partial products that are invariant in the inner loops out to the outer loop...for example, the term IJ*F(i,1) is invariant inside the loops over j,k,l and so could be computed in the i loop and made a temporary there. Then, that term times F(j,I) is invariant inside loops over k,l so could be done at the j loop level and that temporary carried through. Etc., etc., etc., ...
Just moving the inv(J) outside all the loops should be noticeable, however, factoring out the common terms should also be of some benefit.
You'll want to check carefully I didn't get eyes crossed, but I think I have the correct factors in each.
Also, of course, while it's not shown here, the array C should definitely be preallocted.
Related Question