I have the following MATLAB code to calculate the strain energy. I want to optimize this code, because this is where most of my time in an optimization problem is spent.

 stiffness = [6 x 6 matrix] % For a MWE, try magic(6) magic(6)';

 for idr = 1:length(r) % Index of r-location
     for idt = 1:length(theta) % Index of theta location
         e = [exx(idr,idt);eyy(idr,idt);ezz(idr,idt);2*exz(idr,idt);2*eyz(idr,idt);2*exy(idr,idt)]; % strain vector
         E(idr,idt) = e'*stiffness*e; %$ Strain energy at this location
     end
 end

There are two problems here:

• Since length(r)==length(theta)==100, this piece of code gets called 10k times.
• About a third of the time spent in the loop is on assembling e without actually doing any calculation.

I think there should be some ways to speed this up significantly by properly using vectorization, but I'm not entirely sure how to proceed in this particular case.

Note: there are some similar questions, e.g. this, so perhaps there are some ideas that can be translated to MATLAB.

CodePudding user response：

Another option, faster by about a factor ~x2 than @bla's answer,

eall=([exx(:),eyy(:),ezz(:),2*exz(:),2*eyz(:),2*exy(:)]);
E = reshape(sum(stiffness*eall.*eall),[100 100])

taking advantage of the fact that e is a vector so we can do the sum of element-wise multiplication instead of the vector product.

CodePudding user response：

Here is one way of vectorizing your code, doing the vector'*matrix*vector operation on all elements at once by extending the info to the 3rd dimension:

eall=([exx(:),eyy(:),ezz(:),2*exz(:),2*eyz(:),2*exy(:)]);
E=(reshape(pagemtimes(permute(eall,[3 2 1]),pagemtimes(stiffness,permute(eall,[2 3 1]))),[100 100]));

I'm not sure how significant will be the increase of performance, I get ~x5 on the 100x100 example you had...