From: Rob <eu...@ho...> - 2001-08-26 22:19:11
|
The good news is that your routine reduces execution time by 30%. The bad news is that the wrong numbers are coming out of it. I'm trying some otner permutations of "take" and "add.reduce" in the function to see if I can get it. One method that worked was splitting up "c" into c1,c2,c3, such that: c1=Numeric.take(a.NodeCord[:,0],TrngleNode) etc and then using Numeric.add.reduce(s*c1,1) etc This gives the right results, but is slower than plain Python. I'll keep at it. Thanks again. "Paul F. Dubois" wrote: > > def cgqp(QuadPoint, TrngleNode, a): > s = a.Qpnt[Quadpoint,:] > c = Numeric.take(a.NodeCord, TrngleNode) > return Numeric.add.reduce(s * c, axis=1) > > This may or may not be right. The data structures I would have to set up > to test it are too much for Sunday morning. > > -----Original Message----- > From: num...@li... > [mailto:num...@li...] On Behalf Of Rob > Sent: Sunday, August 26, 2001 8:39 AM > To: num...@li.... > Subject: [Numpy-discussion] Can this function by Numpy-ized? > > I finally got my FEM EM code working. I profiled it and this function > uses up a big hunk of time. It performs gaussian integration over a > triangle. I am trying to figure out how to slice the arrays so as to > push it down into the C level. Does anyone have any ideas? Thanks, > Rob. > > ps. it looks to be intractible to me. Maybe I need to look at writing a > C extension. I've never done that before. > > ##********************************************************************** > ***** > ##Prototype: void ComputeGaussQuadPoint(int QuadPoint, int > *a.TrngleNode, > ## double *SrcPointCol) > ##Description: To compute the coordinates of 7-point Gauss nodes of > ## a triangular patch. > ##Input value: > ## int QuadPoint --- node index, it can be from 0 to 6. > ## int *a.TrngleNode --- the three nodes of a tringular patch. > ## double *SrcPointCol --- where to store the results > ##Return value: none > ##Global value used: a.NodeCord, a.Qpnt > ##Global value modified: none > ##Subroutines called: none > ##Note: Not very sure. > ************************************************************************ > **** > def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a): > > SrcPointCol=zeros((3),Float) > > SrcPointCol[0] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\ > + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\ > + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0] > > SrcPointCol[1] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\ > + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\ > + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1] > > SrcPointCol[2] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\ > + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\ > + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2] > > > return SrcPointCol > > > -- > The Numeric Python EM Project > > www.members.home.net/europax > > _______________________________________________ > Numpy-discussion mailing list > Num...@li... > http://lists.sourceforge.net/lists/listinfo/numpy-discussion -- The Numeric Python EM Project www.members.home.net/europax |