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.
|