On Aug 20, 12:11 pm, Jeremy Bailin <astroco...@gmail.com> wrote:
> On Aug 20, 8:22 am, "thomas.jagdhuber" <thomas.jagdhu...@gmail.com>
> wrote:
>
>> Dear experts,
>
>> I would like to matrix multiply two matrices with dimensions
>> [3,3,1500]. means: 1500 times a matrix multiplication of 2 matrices
>> with dimension [3,3]
>> I could do this with a for loop over the dimension [1500] but i
>> suppose this is not very elegant. Is there any other way to do this
>> time-efficient.
>
>> Best regards,
>
>> thomas
>
> My IDL-foo feels strong today... let's see:
>
> If A and B are 3x3, then C = A ## B is equivalent to:
>
> q = rebin(indgen(3),3,3)
> p = rebin(reform(indgen(3),1,3),3,3)
> c = total(rebin(reform(a[q,p],1,3,3),3,3,3) * reform(b[*,q],3,3,3), 2)
>
> Adding an extra dimension on the end of that gets tricky, but I think
> the following should work if A and B are each 3x3xNMATRIX and you want
> C to be a 3x3xNMATRIX array where each C[*,*,i] = A[*,*,i] ##
> B[*,*,i]:
>
> c = total(rebin(reform( (a[q,p,*])[0:2,3*lindgen(3),*],1,3,3,nmatrix),
> 3,3,3,nmatrix) * reform(b[*,q,*],
> 3,3,3,nmatrix), 2)
>
> Incidentally, does anyone have a better way of doing the (a[q,p,*])
> [0:2,3*lindgen(3),*] bit of it? The problem is that if A, Q and P are
> each 3x3 then A[Q,P] is 3x3, but if A is 3x3xN then A[Q,P,*] is 9x9xN.
>
> If I see another reform today, I'm going to scream. ;-)
>
> -Jeremy.
Nevermind I see what you're wanting... sorry for the confusion
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
