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

Rob europax at home.com
Mon Aug 27 19:40:30 CDT 2001

```Hi Tim,  Yes this is part of a FEM-MOM hybrid EM simulator.  It was
originally written in C, so I've had to deal with the C way of doing
things.  I'd like to consolidate many of these types of operations so
that the program becomes streamlined and easier to understand- "The
Python Way" TM.  :)

I will try your method to see how it works.  By contrast, I have a FDTD
simulator which is a true speed demon with Numpy.  But it was organized
in a way that I would easily use slicing to time march the 6 scalar
Maxwell's equations.

In any case I'm having fun.  I do this at home.  At work I have Ansoft

Rob.

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
> 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....
>
>      c = take(NodeCord, TrngleNode)
>      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:
> t1 = clock()
> for i in rng:
> 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
> >
> >     SrcPointCol=zeros((3))
> >
> >
> >
> >
> >
> >     return SrcPointCol
> >
> >
> > #the yet-to-be-faster routine
> >
> >     c= take(NodeCord, TrngleNode)
> > transpose(c),1)
> >
> >     return SrcPointCol
> >
> >
> > print "The Correct:"
> >
> > print "The New"
> >
> >
> >
> >
> >
> > --
> > 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
>
> _______________________________________________
> 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

```