[Numpy-discussion] Trick for fast

Søren Gammelmark gammelmark@gmail....
Fri Feb 3 08:43:25 CST 2012


What about this?

A = einsum("i,ij->", mass, x ** 2)
B = einsum("i,ij,ik->jk", mass, x, x)
I = A * eye(3) - B

/Søren

On 3 February 2012 15:10, <josef.pktd@gmail.com> wrote:

> On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac <alan.isaac@gmail.com> wrote:
> > On 2/3/2012 5:16 AM, santhu kumar wrote:
> >> x = nX3 vector.
> >> mass = nX1 vector
> >> inert = zeros((3,3))
> >> for i in range(n):
> >>        ri = x[i,:].reshape(1,3)
> >>        inert = inert + mass[i,]*(sum(ri*ri)*eye(3) - dot(ri.T,ri))
> >>
> >
> >
> > This should buy you a bit.
> >
> > xdot = (x*x).sum(axis=1)
> > for (massi,xi,xdoti) in zip(mass.flat,x,xdot):
> >       temp = -np.outer(xi,xi)
> >       temp.flat[slice(0,None,4)] += xdoti
> >       inert += massi*temp
> >
> > Alan Isaac
>
> maybe something like this, (self contained example and name spaces to
> make running it easier)
>
> import numpy as np
> n = 15
> x = np.arange(n*3.).reshape(-1,3) #nX3 vector.
> mass = np.linspace(1,2,n)[:,None] #nX1 vector
> inert = np.zeros((3,3))
> for i in range(n):
>      ri = x[i,:].reshape(1,3)
>       inert = inert + mass[i,]*(sum(ri*ri)*np.eye(3) - np.dot(ri.T,ri))
> print inert
>
> print  np.diag((mass * x**2).sum(0))  - np.dot(x.T, mass*x)
>
> [[     0.  -16755.  -17287.5]
>  [-16755.       0.  -17865. ]
>  [-17287.5 -17865.       0. ]]
> [[     0.  -16755.  -17287.5]
>  [-16755.       0.  -17865. ]
>  [-17287.5 -17865.       0. ]]
>
> Josef
>
>
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > NumPy-Discussion@scipy.org
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20120203/d1faa546/attachment.html 


More information about the NumPy-Discussion mailing list