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

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