A VERY slow Matlab code!

A part of my MATLAB code is extremely slow and I would appreciate any help!

Here is the slow function that I have:

function du = odefungPC(u,Dijk,Psi2Norm,Size)
du=zeros(3*Size,1);
for i=1:Size
    for j=1:Size
        for k=1:Size
            du(i)=du(i)+1/Psi2Norm(i)*u(Size+j)*u(2*Size+k)*Dijk{i}(j,k);
            du(Size+i)=du(Size+i)+1/Psi2Norm(i)*u(j)*u(2*Size+k)*Dijk{i}(j,k);
            du(2*Size+i)=du(2*Size+i)-2/Psi2Norm(i)*u(j)*u(Size+k)*Dijk{i}(j,k);
        end
    end
end
end

The dimension of the inputs are as follows:

Size is a scalar equal to 120.

u is of size (3*120) x 1.

Dijk is a cell array of size 120 x 1 and every element of the cell array is a 2D array of size 120 x 120.

Psi2Norm is a vector of size 120 x 1.

This is how I use this function: I want to solve a nonlinear system of ODEs and the above code generates the coefficient matrix that I use later in the following form:

odefun=@(t,u)(odefungPC(u,Dijk,PSI2Norm,Size));
options = odeset('MaxStep',(tn-t0)/(numt));

[T,Ug]=ode45(odefun,tspan,u0,options);

As I said it is a nonlinear system. Therefore, I don't know what u is before solving the system...

Thanks very much in advance!

[Doug]

Can you use a 3D array rather than a cell array for Dijk (Dijk(i,j,k) rather than Dijk{i}(j,k))? On my machine this change resulted in your example running 16 times as fast. There are other things you can do to get additional, though more modest, speedup. For example, you can compute 1/Psi2Norm before the loop, saving ~Size^3 divide operations, and you can compute u(j)/Psi2Norm(i) above the k loop, since it's not a function of k. There are other similar changes you can make.

[Mehdi]

I tried using a 3D array once before, it was slower in my machine! I checked it with a small problem. Even saving the cell array in the 3D array made the code slower...

Answer 1

[Roger Stafford]

In the inner two for-loops you have ordinary matrix multiplication. You should make use of it.

 function du = odefungPC(u,Dijk,Psi2Norm,Size)
 du=zeros(3*Size,1);
 for i=1:Size
   du(i) =           (u(Size+1:2*Size).'*Dijk{i}*u(2*Size+1:3*Size)) /Psi2Norm(i);
   du(i+Size) =      (u(1:Size).'       *Dijk{i}*u(2*Size+1:3*Size)) /Psi2Norm(i);
   du(2*Size+i) = -2*(u(1:Size).'       *Dijk{i}*u(Size+1:2*Size)  ) /Psi2Norm(i);
 end

[Mehdi]

Thanks Roger! I had written this code three years ago and at the time I didn't use it for large data and didn't see the difference much so didn't bother myself. This time when I used it for some large data just realize how slow it is. Unfortunately, I was a bit lazy to look at the details before. I just optimized the code as you said and the speed-up is incredible!

Thanks again!