Re: vectorization quiz

[Jaroslav Hajek]

On Fri, Aug 28, 2009 at 10:44 AM, Carlo de Falco<address@hidden> wrote:
> Hi,
> I've just come back to work after holidays but my brain seems not to
> have noticed,
> I'm sure there is an easy way to vectorize the following code but I
> just can't see it at the moment...
>
> Here is the code:
>
>  for iel = 1:nel
>     for iqn = 1:nqn
>       shg(:,:,iqn,iel) = J(:, :, iqn, iel) * shg(:,:,iqn,iel);
>     endfor
>   endfor
>
> where nel and nqn are possibly very large and size(J(:, :, iqn, iel))
> = size(shg(:,:,iqn,iel)) = [2 2]
>
> Can anyone suggest a way to speed this up?
> Thanks in advance,
> c.

Since the leading dims are so tiny, you could even write the matrix
multiply loop explicitly
for i=1:2; for j = 1:2
  shg1(i,j,:,:) = J(i,1,:,:) .* shg(1,j,:,:) + J(i,2,:,:) .* shg(2,j,:,:)
endfor;endfor
Since each slice on RHS is extracted twice, this could also be sped up
by caching them first in a cell array.

For a more universal solution to multiplying arrays of matrices, see
http://octave.sourceforge.net/doc/f/arraymm.html

regards

-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz