[Numpy-discussion] Need to eliminate a nested loop

eat e.antero.tammi at gmail.com
Wed May 11 08:00:07 EDT 2011


Hi,

On Wed, May 11, 2011 at 8:28 AM, Elfnor <elfnor at gmail.com> wrote:

>
> Hi
>
> The following code produces the desired result but has a slow triple loop
> iterating over the matrix multiplication.
>
> I'm sure it can be eliminated with a neat indexing trick but I can't figure
> out how.
>
> Any suggestions please?
> -----------------------------
> import numpy
> #define domain of function
> x = numpy.linspace(-5,5,64)
> y = numpy.linspace(-5,5,64)
> z = numpy.linspace(-5,5,64)
>
> #calculate f at each point in domain
> a = 5.0
> b = 3.0
> c = 2.0
> #ellipsoid
> E = numpy.array([[1/a**2,   0   ,   0  ,  0  ],
>                [   0   ,1/b**2 ,   0  ,  0  ],
>                [   0   ,   0   ,1/c**2,  0  ],
>                [   0   ,   0   ,   0  , -1  ]])
>
> f = numpy.zeros((x.size,y.size,z.size))
>
> for i,xi in enumerate(x):
>    for j,yj in enumerate(y):
>        for k,zk in enumerate(z):
>            X = numpy.array([xi,yj,zk,1])
>            f[i,j,k] = numpy.dot(numpy.dot(X.T,E),X)
> -----------------------------------
>
Something like this:
n= 64
u= np.linspace(-5, 5, n)[None, :]
u0= u.repeat(n** 2)[None, :]

u1= u.repeat(n)[None, :].repeat(n, axis= 0).reshape(1, -1)

u2= u.repeat(n** 2, axis= 0).reshape(1, -1)

U= np.r_[u0, u1, u2, np.ones((1, n** 3))]

f= (np.dot(E, U)* U).sum(0).reshape(n, n, n)


Regards,eat

>
> Thanks Eleanor
> --
> View this message in context:
> http://old.nabble.com/Need-to-eliminate-a-nested-loop-tp31591457p31591457.html
> Sent from the Numpy-discussion mailing list archive at Nabble.com.
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20110511/6b602523/attachment.html>


More information about the NumPy-Discussion mailing list