Re: vectorization quiz

[Carlo de Falco]

2009/8/29 Jaroslav Hajek <address@hidden>:
> 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

yes that's actually what I came up with myself but I don't like it very much...

> 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

This looks much better! I was not aware of this function, thanks for
the pointer!

> regards

c.