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

Tim Hochberg tim.hochberg at ieee.org
Mon Aug 27 12:10:01 EDT 2001

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.



> 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

More information about the NumPy-Discussion mailing list