Pau, can you confirm that this is the same as the routine you're interested in?
def dot2(a,b): '''Returns dot product of last two dimensions of two 3-D arrays, threaded over first dimension.''' try: assert a.shape[1] == b.shape[2] assert a.shape[0] == b.shape[0] except AssertionError: print "incorrect input shapes" res = zeros( (a.shape[0], a.shape[1], a.shape[1]), dtype=float ) for i in range(a.shape[0]): res[i,...] = dot( a[i,...], b[i,...] ) return res
yes, that is what I would like. I would like it even with more dimensions and with all the broadcasting rules ;-) These can probably be achieved by building actual 'arrays of matrices' (an array of matrix objects) and then using the ufunc machinery. But I think that a simple dot2 function (with an other name of course) will still very useful. pau