I am moving this over to numpy-discussion maillist... I don't have a firm answer for you, but I did notice one issue in your code. You call arange(len(dx) - 1) for your loops, but you probably really need arange(1, len(dx) - 1) because you are accessing elements both after *and* before the current index. An index of -1 is actually valid because that means the last element of the array, and may not be what you intended.
Ben Root On Fri, Jul 2, 2010 at 1:15 PM, Nicolas Bigaouette <nbigaoue...@gmail.com>wrote: > Hi all, > > I don't really know where to ask, so here it is. > > I was able to vectorize the normalization calculation in quantum mechanics: > <phi|phi>. Basically it's a volume integral of a scalar field. Using: > >> norm = 0.0 >> for i in numpy.arange(len(dx)-1): >> for j in numpy.arange(len(dy)-1): >> for k in numpy.arange(len(dz)-1): >> norm += psi[k,j,i]**2 * dx[i] * dy[j] * dz[k] >> > if dead slow. I replaced that with: > >> norm = (psi**2 * >> dx*dy[:,numpy.newaxis]*dz[:,numpy.newaxis,numpy.newaxis]).sum() >> > which is almost instantanious. > > I want to do the same for the calculation of the kinetic energy: > <phi|p^2|phi>/2m. There is a laplacian in the volume integral which > complicates things: > >> K = 0.0 >> for i in numpy.arange(len(dx)-1): >> for j in numpy.arange(len(dy)-1): >> for k in numpy.arange(len(dz)-1): >> K += -0.5 * m * phi[k,j,i] * ( >> (phi[k,j,i-1] - 2.0*phi[k,j,i] + phi[k,j,i+1]) / >> dx[i]**2 >> + (phi[k,j-1,i] - 2.0*phi[k,j,i] + phi[k,j+1,i]) / >> dy[j]**2 >> + (phi[k-1,j,i] - 2.0*phi[k,j,i] + phi[k+1,j,i]) / >> dz[k]**2 >> ) >> > > My question is, how would I vectorize such loops? I don't know how I would > manage the "numpy.newaxis" code-foo with neighbours dependency... Any idea? > > Thanx! > > > ------------------------------------------------------------------------------ > This SF.net email is sponsored by Sprint > What will you do first with EVO, the first 4G phone? > Visit sprint.com/first -- http://p.sf.net/sfu/sprint-com-first > _______________________________________________ > Matplotlib-users mailing list > matplotlib-us...@lists.sourceforge.net > https://lists.sourceforge.net/lists/listinfo/matplotlib-users > >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
