numpy-discussion

[Numpy-discussion] Can this function by Numpy-ized?

[Rob <eu...@ho...>]

I finally got my FEM EM code working.  I profiled it and this function
uses up a big hunk of time.  It performs gaussian integration over a
triangle.  I am trying to figure out how to slice the arrays so as to
push it down into the C level.  Does anyone have any ideas?  Thanks, 
Rob.

ps. it looks to be intractible to me.  Maybe I need to look at writing a
C extension.  I've never done that before.


##***************************************************************************
##Prototype:  void  ComputeGaussQuadPoint(int QuadPoint, int
*a.TrngleNode,
##                                            double *SrcPointCol) 
##Description: To compute the coordinates of 7-point Gauss nodes of 
##             a triangular patch. 
##Input value: 
##    int QuadPoint  --- node index, it can be from 0 to 6. 
##    int *a.TrngleNode  ---  the three nodes of a tringular patch. 
##    double *SrcPointCol --- where to store the results 
##Return value: none 
##Global value used: a.NodeCord, a.Qpnt 
##Global value modified:    none 
##Subroutines called:    none 
##Note: Not very sure.           
****************************************************************************
def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a):
    
    SrcPointCol=zeros((3),Float)
    
    SrcPointCol[0] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\
                     + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\
                     + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0]

    SrcPointCol[1] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\
                     + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\
                     + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1]

    SrcPointCol[2] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\
                     + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\
                     + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2]
                                     

    return SrcPointCol
  





-- 
The Numeric Python EM Project

www.members.home.net/europax

Re: [Numpy-discussion] Can this function by Numpy-ized?

[John J. L. <jj...@po...>]

On Sun, 26 Aug 2001, Rob wrote:

> I finally got my FEM EM code working.  I profiled it and this function
> uses up a big hunk of time.  It performs gaussian integration over a
> triangle.  I am trying to figure out how to slice the arrays so as to
> push it down into the C level.  Does anyone have any ideas?  Thanks,
> Rob.
>
> ps. it looks to be intractible to me.  Maybe I need to look at writing a
> C extension.  I've never done that before.
[...]

I'm not a good enough ufunc / array function hacker to come up with
anything, but I doubt it, and I wonder how much time it would save anyway,
given the size of the arrays involved.  From the look of your comments
before the function, it looks like a) you're a C programmer, not quite
comfortable with Python yet, or b) you wrote it in C first, then moved it
into Python.  If the latter, you might want to try wrapping that C
function with SWIG, though to be honest I'm not sure the overhead of
learning to use SWIG is any less than for writing a C extension manually
(but less error-prone I expect, and less work if you have a lot of C to
wrap).  I think the Scipy project has some typemaps for passing Numeric
arrays; if not, I've seen some collections of SWIG typemaps for Numeric
somewhere...

BTW, the best place to put those comments is in a docstring.  Here is a
slightly more Pythonically-formatted version (I must be bored today):

def ComputeGaussQuadPoint(QuadPoint, a):
    """Return coordinates of 7-point Gauss nodes of a triangular patch.

    QuadPoint: node index, from 0 to 6
    a: triangular patch?
    """
    SrcPointCol=zeros((3),Float)
    tn = a.TrngleNode  # the three nodes of the triangular patch
    SrcPointCol[0] = (a.Qpnt[QuadPoint,0]*a.NodeCord[tn[0],0] +
                      a.Qpnt[QuadPoint,1]*a.NodeCord[tn[1],0] +
                      a.Qpnt[QuadPoint,2]*a.NodeCord[tn[2],0])

    SrcPointCol[1] = (a.Qpnt[QuadPoint,0]*a.NodeCord[tn[0],1] +
                      a.Qpnt[QuadPoint,1]*a.NodeCord[tn[1],1] +
                      a.Qpnt[QuadPoint,2]*a.NodeCord[tn[2],1])

    SrcPointCol[2] = (a.Qpnt[QuadPoint,0]*a.NodeCord[tn[0],2] +
                      a.Qpnt[QuadPoint,1]*a.NodeCord[tn[1],2] +
                      a.Qpnt[QuadPoint,2]*a.NodeCord[tn[2],2])

    return SrcPointCol


John

RE: [Numpy-discussion] Can this function by Numpy-ized?

[Paul F. D. <pa...@pf...>]

def cgqp(QuadPoint, TrngleNode, a):
    s = a.Qpnt[Quadpoint,:]
    c = Numeric.take(a.NodeCord, TrngleNode)
    return Numeric.add.reduce(s * c, axis=1)

This may or may not be right. The data structures I would have to set up
to test it are too much for Sunday morning.

-----Original Message-----
From: num...@li...
[mailto:num...@li...] On Behalf Of Rob
Sent: Sunday, August 26, 2001 8:39 AM
To: num...@li....
Subject: [Numpy-discussion] Can this function by Numpy-ized?

I finally got my FEM EM code working.  I profiled it and this function
uses up a big hunk of time.  It performs gaussian integration over a
triangle.  I am trying to figure out how to slice the arrays so as to
push it down into the C level.  Does anyone have any ideas?  Thanks, 
Rob.

ps. it looks to be intractible to me.  Maybe I need to look at writing a
C extension.  I've never done that before.


##**********************************************************************
*****
##Prototype:  void  ComputeGaussQuadPoint(int QuadPoint, int
*a.TrngleNode,
##                                            double *SrcPointCol) 
##Description: To compute the coordinates of 7-point Gauss nodes of 
##             a triangular patch. 
##Input value: 
##    int QuadPoint  --- node index, it can be from 0 to 6. 
##    int *a.TrngleNode  ---  the three nodes of a tringular patch. 
##    double *SrcPointCol --- where to store the results 
##Return value: none 
##Global value used: a.NodeCord, a.Qpnt 
##Global value modified:    none 
##Subroutines called:    none 
##Note: Not very sure.           
************************************************************************
****
def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a):
    
    SrcPointCol=zeros((3),Float)
    
    SrcPointCol[0] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\
                     + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\
                     + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0]

    SrcPointCol[1] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\
                     + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\
                     + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1]

    SrcPointCol[2] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\
                     + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\
                     + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2]
                                     

    return SrcPointCol
  





-- 
The Numeric Python EM Project

www.members.home.net/europax

_______________________________________________
Numpy-discussion mailing list
Num...@li...
http://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] Can this function by Numpy-ized?

[Rob <eu...@ho...>]

Thanks Paul,  I will try it out.  That at least gives me some
direction.  I've been pouring over the Numpy manual, but there are so
many different functions!  I would like to try to write an extension
that includes some of these FEM routines, but thats for later.  

BTW, at work I have to use Windows NT, but I've grown to love the
calldll/DynWin extension.  You can even interact with the windows kernel
with that one.   Rob.


"Paul F. Dubois" wrote:
> 
> def cgqp(QuadPoint, TrngleNode, a):
>     s = a.Qpnt[Quadpoint,:]
>     c = Numeric.take(a.NodeCord, TrngleNode)
>     return Numeric.add.reduce(s * c, axis=1)
> 
> This may or may not be right. The data structures I would have to set up
> to test it are too much for Sunday morning.
> 
> -----Original Message-----
> From: num...@li...
> [mailto:num...@li...] On Behalf Of Rob
> Sent: Sunday, August 26, 2001 8:39 AM
> To: num...@li....
> Subject: [Numpy-discussion] Can this function by Numpy-ized?
> 
> I finally got my FEM EM code working.  I profiled it and this function
> uses up a big hunk of time.  It performs gaussian integration over a
> triangle.  I am trying to figure out how to slice the arrays so as to
> push it down into the C level.  Does anyone have any ideas?  Thanks,
> Rob.
> 
> ps. it looks to be intractible to me.  Maybe I need to look at writing a
> C extension.  I've never done that before.
> 
> ##**********************************************************************
> *****
> ##Prototype:  void  ComputeGaussQuadPoint(int QuadPoint, int
> *a.TrngleNode,
> ##                                            double *SrcPointCol)
> ##Description: To compute the coordinates of 7-point Gauss nodes of
> ##             a triangular patch.
> ##Input value:
> ##    int QuadPoint  --- node index, it can be from 0 to 6.
> ##    int *a.TrngleNode  ---  the three nodes of a tringular patch.
> ##    double *SrcPointCol --- where to store the results
> ##Return value: none
> ##Global value used: a.NodeCord, a.Qpnt
> ##Global value modified:    none
> ##Subroutines called:    none
> ##Note: Not very sure.
> ************************************************************************
> ****
> def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a):
> 
>     SrcPointCol=zeros((3),Float)
> 
>     SrcPointCol[0] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\
>                      + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\
>                      + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0]
> 
>     SrcPointCol[1] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\
>                      + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\
>                      + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1]
> 
>     SrcPointCol[2] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\
>                      + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\
>                      + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2]
> 
> 
>     return SrcPointCol
> 
> 
> --
> The Numeric Python EM Project
> 
> www.members.home.net/europax
> 
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> http://lists.sourceforge.net/lists/listinfo/numpy-discussion
> 
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> http://lists.sourceforge.net/lists/listinfo/numpy-discussion

-- 
The Numeric Python EM Project

www.members.home.net/europax

Re: [Numpy-discussion] Can this function by Numpy-ized?

[Rob <eu...@ho...>]

The good news is that your routine reduces execution time by 30%.  The
bad news is that the wrong numbers are coming out of it.  I'm trying
some otner permutations of "take" and "add.reduce" in the function to
see if I can get it.  One method that worked was splitting up "c" into
c1,c2,c3, such that:

c1=Numeric.take(a.NodeCord[:,0],TrngleNode)
etc

and then using  Numeric.add.reduce(s*c1,1)
etc

This gives the right results, but is slower than plain Python.  I'll
keep at it.  Thanks again.


"Paul F. Dubois" wrote:
> 
> def cgqp(QuadPoint, TrngleNode, a):
>     s = a.Qpnt[Quadpoint,:]
>     c = Numeric.take(a.NodeCord, TrngleNode)
>     return Numeric.add.reduce(s * c, axis=1)
> 
> This may or may not be right. The data structures I would have to set up
> to test it are too much for Sunday morning.
> 
> -----Original Message-----
> From: num...@li...
> [mailto:num...@li...] On Behalf Of Rob
> Sent: Sunday, August 26, 2001 8:39 AM
> To: num...@li....
> Subject: [Numpy-discussion] Can this function by Numpy-ized?
> 
> I finally got my FEM EM code working.  I profiled it and this function
> uses up a big hunk of time.  It performs gaussian integration over a
> triangle.  I am trying to figure out how to slice the arrays so as to
> push it down into the C level.  Does anyone have any ideas?  Thanks,
> Rob.
> 
> ps. it looks to be intractible to me.  Maybe I need to look at writing a
> C extension.  I've never done that before.
> 
> ##**********************************************************************
> *****
> ##Prototype:  void  ComputeGaussQuadPoint(int QuadPoint, int
> *a.TrngleNode,
> ##                                            double *SrcPointCol)
> ##Description: To compute the coordinates of 7-point Gauss nodes of
> ##             a triangular patch.
> ##Input value:
> ##    int QuadPoint  --- node index, it can be from 0 to 6.
> ##    int *a.TrngleNode  ---  the three nodes of a tringular patch.
> ##    double *SrcPointCol --- where to store the results
> ##Return value: none
> ##Global value used: a.NodeCord, a.Qpnt
> ##Global value modified:    none
> ##Subroutines called:    none
> ##Note: Not very sure.
> ************************************************************************
> ****
> def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a):
> 
>     SrcPointCol=zeros((3),Float)
> 
>     SrcPointCol[0] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\
>                      + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\
>                      + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0]
> 
>     SrcPointCol[1] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\
>                      + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\
>                      + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1]
> 
>     SrcPointCol[2] =   a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\
>                      + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\
>                      + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2]
> 
> 
>     return SrcPointCol
> 
> 
> --
> The Numeric Python EM Project
> 
> www.members.home.net/europax
> 
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> http://lists.sourceforge.net/lists/listinfo/numpy-discussion

-- 
The Numeric Python EM Project

www.members.home.net/europax