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

28 Aug 2001

      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" 
To: 
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@lists.sourceforge.net
http://lists.sourceforge.net/lists/listinfo/numpy-discussion
-- 
The Numeric Python EM Project

www.members.home.net/europax