[Numpy-discussion] Gaussian Quadrature routine Numpy-ization :)

[Rob]

Tim, I used your routine, computing all of the triangles at once and
then feeding the resulting array back into the needy functions.  I
shaved 30% off the execution time!!  There is another routine which
computes correction terms which can use the "take()" method.  I'll work
on that tomorrow.  Thanks alot!

Rob.

ps. I'll have to study the Numpy manual to figure out how your routine
works.


Tim Hochberg wrote:
> 
> Hi Rob,
> 
> It looks like you are trying to use the function below to integrate over a
> single triangle. I'm almost certain that you will _not_ be able to get this
> to run fast; there's just too much overhead per triangle from all the
> function calls and what not. Instead you need to structure things so that
> you do more work per call. One way is to pass a list of triangles and get
> back a list of answers. This way you spread your area over many triangles.
> Another possibility, assuming that you need to evaluate the integral for all
> seven possible values of QPoint each time, is to compute the answer for all
> values of QPoint at once so that you reduce the overhead per triangle by
> seven. For example (minimally tested):
> 
> Qpnt=reshape(arange(7*3), (7,3))
> # Other code from other messages elided....
> 
> def Quads(TrngleNode, Qpnt, NodeCord):
>      c = take(NodeCord, TrngleNode)
>      SrcPointCols= add.reduce(Qpnt[...,NewAxis] *  c[NewAxis,],1)
>      return SrcPointCols
> 
> # Initialize stuff to minimize effects of timing....
> from time import clock
> q1 = [None]*7
> q2 = [None]*7
> rng = range(7)
> 
> # Time the three methods....
> 
> t0 = clock()
> for i in rng:
>     q1[i] = ComputeGaussQuadPoint(i,TrngleNode,Qpnt,NodeCord)
> t1 = clock()
> for i in rng:
>     q2[i] = Quad(i,TrngleNode,Qpnt,NodeCord)
> t2 = clock()
> q3 = Quads(TrngleNode, Qpnt, NodeCord)
> t3 = clock()
> 
> # And the results...
> 
> print "Times (us):", (t1-t0)*1e6, (t2-t1)*1e6, (t3-t2)*1e6
> for i in range(7):
>     print "The Answers:", q1[i], q2[i], q3[i]
> 
> If this code is still a bottleneck, you could do both (compute the values
> for all nodes and all values of QPoint at once, but this may be enough to
> move the bottleneck elsewhere; by my timing it's over 7 times as fast.
> 
> -tim
> 
> [snip]
> 
> > Evidently transpose is not very fast even for
> > smal matrices.
> 
> Actually, transpose should be fast, and should look progressively faster for
> larger matrices. Typically all that happens is that the strides are changed.
> 
> ----- Original Message -----
> From: "Rob" <europax at home.com>
> To: <numpy-discussion at lists.sourceforge.net>
> Sent: Monday, August 27, 2001 7:49 AM
> Subject: [Numpy-discussion] Gaussian Quadrature routine Numpy-ization :)
> 
> > I finally got it to work, but the Numpy-ized version runs slower than
> > the plain Python one.  I think that I can transpose the NodeCord matrix
> > once in the program and feed that in, rather than the scratch matrix
> > that is generated here.  Evidently transpose is not very fast even for
> > smal matrices. Here is my test program:
> >
> >
> > from Numeric import *
> >
> >
> >
> > Qpnt=array(([20,21,22],[23,24,25],[26,27,28]))
> > NodeCord=array(([1,2,3],[4,5,6],[7,8,9]))
> > TrngleNode=array((1,2,0))
> >
> >
> > #the original routine
> > def ComputeGaussQuadPoint(QuadPoint,TrngleNode,Qpnt,NodeCord):
> >
> >     SrcPointCol=zeros((3))
> >
> >
> >     SrcPointCol[0] =   Qpnt[QuadPoint,0]*NodeCord[TrngleNode[0],0]\
> >                      + Qpnt[QuadPoint,1]*NodeCord[TrngleNode[1],0]\
> >                      + Qpnt[QuadPoint,2]*NodeCord[TrngleNode[2],0]
> >
> >     SrcPointCol[1] =   Qpnt[QuadPoint,0]*NodeCord[TrngleNode[0],1]\
> >                      + Qpnt[QuadPoint,1]*NodeCord[TrngleNode[1],1]\
> >                      + Qpnt[QuadPoint,2]*NodeCord[TrngleNode[2],1]
> >
> >     SrcPointCol[2] =   Qpnt[QuadPoint,0]*NodeCord[TrngleNode[0],2]\
> >                      + Qpnt[QuadPoint,1]*NodeCord[TrngleNode[1],2]\
> >                      + Qpnt[QuadPoint,2]*NodeCord[TrngleNode[2],2]
> >
> >     return SrcPointCol
> >
> >
> > #the yet-to-be-faster routine
> >
> > def Quad(QuadPoint, TrngleNode, Qpnt,NodeCord):
> >     s = Qpnt[QuadPoint,:]
> >     c= take(NodeCord, TrngleNode)
> >     SrcPointCol= add.reduce(s *
> > transpose(c),1)
> >
> >     return SrcPointCol
> >
> > QuadPoint=1
> >
> > print "The Correct:"
> > print ComputeGaussQuadPoint(QuadPoint,TrngleNode,Qpnt,NodeCord)
> >
> > print "The New"
> > print  Quad(QuadPoint,TrngleNode,Qpnt,NodeCord)
> >
> >
> >
> >
> >
> > --
> > The Numeric Python EM Project
> >
> > www.members.home.net/europax
> >
> > _______________________________________________
> > Numpy-discussion mailing list
> > Numpy-discussion at lists.sourceforge.net
> > http://lists.sourceforge.net/lists/listinfo/numpy-discussion

-- 
The Numeric Python EM Project

www.members.home.net/europax