numpy-discussion — Discussion list for all users of Numerical Python

[Numpy-discussion] Happy New Year Nummies -- and some discussion about the future

[Paul F. D. <pau...@ho...>]

A millenium-end report from the Head Nummie (this name is a joke; see the
DEVELOPERS file):

There have been a steady set of messages on the subject of I should do this
or that to make it easier to make RPMs.  It is impossible for me to act on
these: I don't know much about RPMs, and if I did, I don't know if making
the change suggested is good or bad for someone doing something else, like
making Windows installers. Therefore my policy is to rely on the Distutils
people to work this out. Those who wish to make it easier to make a binary
installer for platform xyz should figure out what would be required by the
Distutils bdist family of commands.

That is not to say that I don't appreciate people trying to help. I'm
grateful for all the support I get from the community. I think that relying
on division of labor in this case is the right thing to do, so that we take
advantage of the Distutils effort. If I'm wrong, I'll listen.

There are a number of bug reports on the sourceforge site. I would be
grateful for patches. In particular there are two reports dealing with FFTs.
I lack the expertise to deal with these.

The masked array package MA has been getting more of a workout as I put it
into production at my site. I believe that it fills not only the immediate
need for dealing with missing values, but can serve as a model for how to
make a "Numeric-workalike" with extra properties, since we can't inherit
from Numeric's array. Since MA is improved fairly often, I recommend keeping
up via cvs if you are a user.

I have new manuals but have had difficulty with the transport up to
SourceForge. Anybody else having such a problem? I used scp from a Linux box
and it sat there and didn't do anything.

The rest of this is for developers.

Actually, once you get into it, it isn't all that clear that inheritance
would help very much. For example, suppose you want an array class F that
has masked values but also has a special behavior f() controlled by a
parameter set at creation, beta. Suppose therefore you decide to inherit
from class MA. Thus the constructor of your new class must take the same
arguments as an MA array but add a beta=somedefault. OK, we do that. Now we
can construct a couple of F's:
f1 = F([1.,2.,3.], beta=1.)
f2 = F([4.,2.,1.], beta=2.)
Great. Now we can do f1.f(), f2.f(). Maybe we redefine __str__ so we can
print beta when we print f1.

Now try to do something. Anything. Say,
f3 = f1 + f2

Oops. f3 is an MA, not an F. We might have written __add__ in MA so the it
used self.__class__ to construct the answer. But since the constructor now
needs a new parameter, there is no way MA.__add__ can make an F. It doesn't
know how. Doesn't know how to call F(), doesn't know what value to use for
beta anyway.

So now we redefine all the methods in MA. Besides muttering that maybe
inheriting didn't buy me a lot, I am still nowhere, for the next thing I
realize is that every function f(a) that takes an MA as an argument and
returns an MA, still returns an MA. If any of these make sense for an
instance of F, I have to replace them, just as MA replaced sin, cos, sqrt,
take, etc. from Numeric.

I have therefore come to the conclusion that we have been barking up the
wrong tree. There might be a few cases where inheritance would buy you
something, but essentially Numeric and MA are useless as parents. Instead,
what would be good to have is a python class Numeric built upon a suite of C
routines that did all the real work, such as adding, transposing, iterating
over elements, applying a function to each element, etc. Since it is now
possible to build a Python class with C methods, which it was not when
Numeric was created, we ought to think about it. Such an API could be used
to make other classes with good performance. We could lose the artificial
layer that is in there now that makes it so tedious to add a function. (I
counted something like five or six places I had to modify when I added
"put".)

Best regards to all,

Paul Dubois

[Numpy-discussion] CXX versions

[<pf...@mm...>]

Hi,

I'm playing around with CXX because I'm going to be writing an extension in
the next few days.  I'm on WinNT running Microsoft Visual C++ 6.0 with
service pack 3 installed.  Also, Python 2.0.  I've gotten all of CXX-5.0b
to compile, and I compiled the demo.  This required a slight edit of
cxx_extensions.cxx which was documented in the SourceForge buglist for CXX.
I was pleasantly surprised that VC++ was able to compile this code.
However, when I import CXX.example python crashes.  Is this a known
problem, is there a workaround, or should I try CXX-4.2 instead?  For some
reason, I am really adverse to writing another extension using C.

Tom Fenn

[Numpy-discussion] Numpy doesn't build properly

[Tony S. <ant...@ie...>]

Because the Numpy sub-packages use '#include "Numeric/arratobject.h", the
Numeric package will not build from scratch.  This may not matter if one
does a 'python setup_all.py install' as long as the base Numeric package is
installed before the sub-packages are built.  I don't think that this is
what setup_all.py does because the core setup.py is called with no
arguments.  In any case, I'm trying to build RPMs and the build and install
processes are seperated.

In order to build the packages seprately I am now putting a symbolic link in
each package's Include directory called Numeric that points to the Include
directory for the core package.  I would like to suggest that the core
header files be put in a directory 'Include/Numeric' rather than directly in
'Include'.  That way the header file locations more acurately mirror the way
that they are installed.

Tony

[Numpy-discussion] 17.2.0

[Paul F. D. <pau...@ho...>]

I uploaded Numeric 17.2.0 but then got errors from Sourceforge's web pages
for actually releasing it. I'm going to wait a day or two and try again. I
may be forced to renumber around the one that is already there. I don't
think users can actually see this file yet.

CVS is up to date with what will become 17.2.? when I finally succeed.

Re: [Numpy-discussion] Tutorial in Numerical Python? (fwd)

[Herb S. <Her...@gr...>]

>---------- Forwarded message ----------
>Date: 13 Dec 2000 11:27:08 +0000
>From: José Luis Gómez Dans <j.l...@sh...>
>Newsgroups: comp.lang.python
>Subject: Tutorial in Numerical Python?
>
>
>Hi,
>         I am wondering whether such a thing exists. Is there a nice
>tutorial for NumPy that guides you through the basic steps? Say, a
>little hands-on introduction to Linear Algebra, Transforms and
>possibly the inclusion of some of the graphic modules?

Hi,

  Take a look at: <http://www.wag.caltech.edu/home/rpm/python_course/>
-- 
---------------------------------
Herb Schilling
NASA Glenn Research Center
Cleveland, OH 44135
Her...@gr...

If all our misfortunes were laid in one common heap whence everyone must take
an equal portion, most people would be contented to take their own and
depart. -Socrates (469?-399 B.C.)

[Numpy-discussion] Tutorial in Numerical Python? (fwd)

[Fredrik S. <su9...@na...>]

Hi,

	I took the liberty to forward your email to
num...@li....

The only document i know of is located at the Python
homepage http://numpy.sourceforge.net.

/fredrik

---------- Forwarded message ----------
Date: 13 Dec 2000 11:27:08 +0000
From: José Luis Gómez Dans <j.l...@sh...>
Newsgroups: comp.lang.python
Subject: Tutorial in Numerical Python?


Hi,
        I am wondering whether such a thing exists. Is there a nice
tutorial for NumPy that guides you through the basic steps? Say, a
little hands-on introduction to Linear Algebra, Transforms and
possibly the inclusion of some of the graphic modules?

        There are some willing students to be converted around these
parts, and a hands-on tutorial would be great. Unfortunately, I am a
beginner in these waters myself, but can see the opportunity for this.

        Thanks in advance,
José
-- 
José L Gómez Dans			PhD student
Tel: +44 114 222 5582			Radar & Communications Group
FAX; +44 870 132 2990			Department of Electronic Engineering
					University of Sheffield UK

Re: [Numpy-discussion] matrix file format

[Konrad H. <hi...@cn...>]

> I'm trying to find out what format a matrix input file should have in
> numpy.

Any format you like, but you have to supply the I/O code yourself.
Which however is really simple. For example, the following
function from ScientificPython:

  from Scientific.IO.TextFile import TextFile
  import string, Numeric

  def readArray(filename):
      "Return an array containing the data from file |filename|."
      data = []
      for line in TextFile(filename):
	  if line[0] != '#':
	      data.append(map(eval, string.split(line)))
      a = Numeric.array(data)
      if a.shape[0] == 1 or a.shape[1] == 1:
	  a = Numeric.ravel(a)
      return a

Konrad.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: hi...@cn...
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------

[Numpy-discussion] matrix file format

[Christopher L. <cl...@gn...>]

You have many options:  

I frequently use Konrad Hinsen's "Scientific" package
(http://dirac.cnrs-orleans.fr/programs/scientific.html).  It make importing
files from user input or other packages like matlab easy and supports
several file formats.  The simpliest is a basic white-space-delimitted text file

if you have a text file "arrayfile.txt" that looks like this:

# arrayfile.txt (lines beginning  with '#' are ignored)
1 2 3
4 5 6
6 8 9

You can fire up python and read in your matrix like this:

   >>> from Scientific.IO.ArrayIO import readArray
   >>> A = readArray('arrayfile.txt')
   >>> print A
   [[1 2 3]
    [4 5 6]
    [6 8 9]]

Re: [Numpy-discussion] matrix file format

[Janko H. <jh...@if...>]

In principle there is no such format in the sense of other numerical
environments. But it is possible to simple put the data into a python
module and import that, like

# Datafile data.py
import Numeric
data=Numeric.array([[1,2,3],[4,5,6],[7,8,9]])

# in python (other module or interactive session)
from data import data
print data

HTH
__Janko


Peter I. Hansen writes:
 > Hi,
 > 
 > I'm trying to find out what format a matrix input file should have in
 > numpy.
 > example :
 > 
 > I have the matrix :
 > [[1 2 3
 >    4 5 6
 >    7 8 9]]
 > 
 > and want to read it to nympy from a text-filen that is easy for a user
 > to generate in an editor.
 > 
 > Can you tell me how to do that. I read the manual and it didnøt given
 > much af on idea.
 > 
 > I thank you in advance
 > 
 > Peter I. Hansen
 > 
 > _______________________________________________
 > Numpy-discussion mailing list
 > Num...@li...
 > http://lists.sourceforge.net/mailman/listinfo/numpy-discussion

[Numpy-discussion] matrix file format

[Peter I. H. <pi...@oe...>]

Hi,

I'm trying to find out what format a matrix input file should have in
numpy.
example :

I have the matrix :
[[1 2 3
   4 5 6
   7 8 9]]

and want to read it to nympy from a text-filen that is easy for a user
to generate in an editor.

Can you tell me how to do that. I read the manual and it didnøt given
much af on idea.

I thank you in advance

Peter I. Hansen

[Numpy-discussion] CVS updated but cannot release 17.2.0

[Paul F. D. <pau...@ho...>]

I have CVS updated with fixes to MA, and the optimization for ones, etc.

I can't release the files because they lost the upload server that does that
and until they get the replacement online there are no new file releases.

Developers, there are bug reports for FFT and savespace that need attention
from the corresponding experts.

The patch list is clear.

We have completed the refereed papers program for IPC9 and it will be
released shortly after notification of authors. The program is quite good,
plus there are tutorials and two other tracks. I encourage all Nummies to
attend (March, Long Beach CA; see python9.org). Registration has just
opened.

[Numpy-discussion] Bug in MA

[Paul F. D. <du...@ll...>]

It has come to my attention that the Masked Array package, MA, has a rather
severe bug. The bug occurs when doing a mixed scalar/array operation in which
the array has its first element masked. The bug causes the value computed by
this operation to be wrong. I should have a fix for this problem shortly but
since this bug is so severe I thought I should put out a warning immediately.

This bug was introduced a few versions ago. 

Re: [Numpy-developers] Re: [Numpy-discussion] RE: Numpy

[Travis O. <Oli...@ma...>]

> I agree -- I also had the thought of changing to an if statement after
> I sent my mail.  Much better than the switch since it can't break.

I took the liberty to change this in the CVS tree.  (I haven't been very
active lately and wanted to help in some way.)

I haven't released a new version with the change, though.

-Travis

Re: [Numpy-discussion] RE: Numpy

[<rl...@st...>]

>From ta...@rh... Wed Dec  6 13:23 EST 2000
>
>I think a more portable solution is to just make it all an if statement.
>I'm not comfortable with presuming that floats are always 4 bytes.

I agree -- I also had the thought of changing to an if statement after
I sent my mail.  Much better than the switch since it can't break.

>I've thrown in the int and float cases to make sure we cover every
>machine that may be out there.

I suspect that if float is 8 bytes then int will be too.  Just char,
short, float, and double should suffice.  If the 4-byte version is
missing it will just fall back on the default (slower but otherwise
OK), and if the C compiler doesn't support shorter data items then
probably the user won't be using them anyway...
					Rick White

Re: [Numpy-discussion] RE: Numpy

[Lee T. <ta...@rh...>]

I think a more portable solution is to just make it all an if statement.
I'm not comfortable with presuming that floats are always 4 bytes.


        if (src_nd == 0 && dest_nd == 0) {
            if (elsize == sizeof(char)) {
                memset(dest, *src, copies);
            } else if (elsize == sizeof(short)) {
                for(j=copies; j; --j, dest += sizeof(short))
                    *(short*)dest = *(short*)src;
            } else if (elsize == sizeof(int)) {
                for(j=copies; j; --j, dest += sizeof(int))
                    *(int*)dest = *(int*)src;
            } else if (elsize == sizeof(long)) {
                for(j=copies; j; --j, dest += sizeof(long))
                    *(long*)dest = *(long*)src;
            } else if (elsize == sizeof(float)) {
                for(j=copies; j; --j, dest += sizeof(float))
                    *(float*)dest = *(float*)src;
            } else if (elsize == sizeof(double)) {
                for(j=copies; j; --j, dest += sizeof(double))
                    *(double*)dest = *(double*)src;
            } else {
                for(j=copies; j; --j, dest += elsize)
                    memcpy(dest, src, elsize);
            }


I've thrown in the int and float cases to make sure we cover every
machine that may be out there.

Lee Taylor

On Wed, 6 Dec 2000 rl...@st... wrote:

> A suggestion for an easy way to fix this problem: use float instead of
> long.  Presumably floats are always 4 bytes.
> 				Rick
> 
> Paul Dubois writes:
> 
> >This was a patch that I put in before making this release. If you get 17.1.1
> >it won't be in there.
> >
> > [...]
> >
> >cc: Error: Src/arrayobject.c, line 79: The switch statement containing this
> >case label already has a case label for "8". (dupcase)
> >            case sizeof(double):
> >------------^
> >error: command 'cc' failed with exit status 1
> >tc05:/usr/tmp/chase/Numeric-17.1.2[408]
> >
> >
> >int do_sliced_copy(char *dest, int *dest_strides, int *dest_dimensions,
> >                   int dest_nd, char *src, int *src_strides,
> >                   int *src_dimensions, int src_nd, int elsize,
> >                   int copies) {
> >        int i, j;
> >
> >        if (src_nd == 0 && dest_nd == 0) {
> >            switch(elsize) {
> >            case sizeof(char):
> >                memset(dest, *src, copies);
> >                break;
> >            case sizeof(short):
> >                for(j=copies; j; --j, dest += sizeof(short))
> >                    *(short*)dest = *(short*)src;
> >                break;
> >            case sizeof(long):
> >                for(j=copies; j; --j, dest += sizeof(int))
> >                    *(int*)dest = *(int*)src;
> >                break;
> >            case sizeof(double):
> >                for(j=copies; j; --j, dest += sizeof(double))
> >                    *(double*)dest = *(double*)src;
> >                break;
> >            default:
> >                for(j=copies; j; --j, dest += elsize)
> >                    memcpy(dest, src, elsize);
> >            }
> >            return 0;
> >        }
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> http://lists.sourceforge.net/mailman/listinfo/numpy-discussion
> 

Re: [Numpy-discussion] RE: Numpy

[<rl...@st...>]

A suggestion for an easy way to fix this problem: use float instead of
long.  Presumably floats are always 4 bytes.
				Rick

Paul Dubois writes:

>This was a patch that I put in before making this release. If you get 17.1.1
>it won't be in there.
>
> [...]
>
>cc: Error: Src/arrayobject.c, line 79: The switch statement containing this
>case label already has a case label for "8". (dupcase)
>            case sizeof(double):
>------------^
>error: command 'cc' failed with exit status 1
>tc05:/usr/tmp/chase/Numeric-17.1.2[408]
>
>
>int do_sliced_copy(char *dest, int *dest_strides, int *dest_dimensions,
>                   int dest_nd, char *src, int *src_strides,
>                   int *src_dimensions, int src_nd, int elsize,
>                   int copies) {
>        int i, j;
>
>        if (src_nd == 0 && dest_nd == 0) {
>            switch(elsize) {
>            case sizeof(char):
>                memset(dest, *src, copies);
>                break;
>            case sizeof(short):
>                for(j=copies; j; --j, dest += sizeof(short))
>                    *(short*)dest = *(short*)src;
>                break;
>            case sizeof(long):
>                for(j=copies; j; --j, dest += sizeof(int))
>                    *(int*)dest = *(int*)src;
>                break;
>            case sizeof(double):
>                for(j=copies; j; --j, dest += sizeof(double))
>                    *(double*)dest = *(double*)src;
>                break;
>            default:
>                for(j=copies; j; --j, dest += elsize)
>                    memcpy(dest, src, elsize);
>            }
>            return 0;
>        }

[Numpy-discussion] RE: Numpy

[Paul F. D. <du...@us...>]

This was a patch that I put in before making this release. If you get 17.1.1
it won't be in there.

Developers -- I'm busy right now with the Python conference; can someone fix
this for me and make a new distribution?

-----Original Message-----
From: Lila Chase [mailto:ch...@tc...]
Sent: Tuesday, December 05, 2000 3:45 PM
To: du...@ll...
Subject: Numpy


Paul,  First attempt at installing Numeric Python on Teracluster; fails.
Lila

tc05:/usr/tmp/chase/Numeric-17.1.2[407] !393
python setup_all.py install
running install
running build
running build_py
not copying Lib/ArrayPrinter.py (output up-to-date)
not copying Lib/Numeric.py (output up-to-date)
not copying Lib/Precision.py (output up-to-date)
not copying Lib/UserArray.py (output up-to-date)
not copying Lib/numeric_version.py (output up-to-date)
running build_ext
building '_numpy' extension
skipping Src/_numpymodule.c
(build/temp.osf1-V4.0-alpha-2.0/Src/_numpymodule.o up-to-date)
cc -O -Olimit 1500 -IInclude -I/g/g16/chase/include/python2.0 -c
Src/arrayobject.c -o build/temp.osf1-V4.0-alpha-2.0/Src/arrayobject.o
cc: Error: Src/arrayobject.c, line 79: The switch statement containing this
case label already has a case label for "8". (dupcase)
            case sizeof(double):
------------^
error: command 'cc' failed with exit status 1
tc05:/usr/tmp/chase/Numeric-17.1.2[408]


int do_sliced_copy(char *dest, int *dest_strides, int *dest_dimensions,
                   int dest_nd, char *src, int *src_strides,
                   int *src_dimensions, int src_nd, int elsize,
                   int copies) {
        int i, j;

        if (src_nd == 0 && dest_nd == 0) {
            switch(elsize) {
            case sizeof(char):
                memset(dest, *src, copies);
                break;
            case sizeof(short):
                for(j=copies; j; --j, dest += sizeof(short))
                    *(short*)dest = *(short*)src;
                break;
            case sizeof(long):
                for(j=copies; j; --j, dest += sizeof(int))
                    *(int*)dest = *(int*)src;
                break;
            case sizeof(double):
                for(j=copies; j; --j, dest += sizeof(double))
                    *(double*)dest = *(double*)src;
                break;
            default:
                for(j=copies; j; --j, dest += elsize)
                    memcpy(dest, src, elsize);
            }
            return 0;
        }

[Numpy-discussion] numpy on alpha/64bit

[Matthias K. <do...@cs...>]

Doh, found a problem and am still trying to figure out what the best way
to fix it is.  In Src/arrayobject.c, starting at around line 70:

int do_sliced_copy(char *dest, int *dest_strides, int *dest_dimensions,
                   int dest_nd, char *src, int *src_strides,
                   int *src_dimensions, int src_nd, int elsize,
                   int copies) {
        int i, j;               

        if (src_nd == 0 && dest_nd == 0) {
            switch(elsize) {
            case sizeof(char):
                memset(dest, *src, copies);
                break;                
            case sizeof(short):
                for(j=copies; j; --j, dest += sizeof(short))
                    *(short*)dest = *(short*)src;               
                break;
            case sizeof(long):
                for(j=copies; j; --j, dest += sizeof(int))
                    *(int*)dest = *(int*)src;
                break;
            case sizeof(double):
                for(j=copies; j; --j, dest += sizeof(double))
                    *(double*)dest = *(double*)src;
                break;  
            default: 
                for(j=copies; j; --j, dest += elsize)
                    memcpy(dest, src, elsize);
            }
            return 0;
        }

The problem is, that sizeof(long) = sizeof(double) on Alpha (and probably
on the other 64-bit archs).  So, what do you think the best method of
handling this is?

I don't know the package well anymore or else I'd fix it :-P

C

[Numpy-discussion] Re: [Python-Help] Numeric bug, or my misunderstanding?

[Daniel Y. <dy...@hk...>]

On Fri, 1 Dec 2000, chad scherrer wrote:

> >>> import Numeric
> >>> def mod2(x): return x%2
> >>> apply(mod2, [-3])
> 1
> >>> apply(mod2, Numeric.array([-3]))
> 1
> >>> apply(mod2, Numeric.array([[-3]]))
> array([-1])
> 
> Should that not be array([1]) ?

Hmmm.  Good point; This looks weird to me too.  You should probably talk
to the Numeric Python people at the Numpy-discussion newsgroup:

    http://lists.sourceforge.net/mailman/listinfo/numpy-discussion

I'll try to forward this message to that group, so perhaps they'll be able
to figure out if it's a true bug.

I hope that your answer will be answered quickly.  Good luck!

RE: [Numpy-discussion] matrixmultiply implementation

[Paul F. D. <pau...@ho...>]

Not to put a damper on anything, but it would be a big change to REQUIRE
dgemm in the Numeric core.
The issues we shunted aside to the subpackage would return.

-----Original Message-----
From: num...@li...
[mailto:num...@li...]On Behalf Of
Christopher Lee
Sent: Friday, December 01, 2000 4:43 PM
To: num...@li...
Subject: [Numpy-discussion] matrixmultiply implementation



I'm curious to see if anyone has looked at the implementation of
dot/matrixmultiply recently.  I was doing a benchmark lapack-based matrix
inversion and was surprised to find that the accuracy check was taking
longer than the inversion.  The bottleneck was the matrix multiply.  It
appears that a blas dgemm implementation might give a large improvement.
Here is a sample benchmark using PyLapack's dblas.dgemm for comparison.
Results for a 1000 x 1000 matrix inversion:

  time elapsed for inverse (sec) 10.050768
  time elapsed for matrixmultiply (sec) 24.142714
  time elapsed for dgemm-matrixmultiply (sec) 5.997188
  max error is: 3.30398486348e-10

Given the "dot" code doesn't look too bad, so I might be able to add dgemm
support myself, at least for some cases.

-chris

p.s. I'll include the test code below
#####################################
import time
from Numeric import *
from LinearAlgebra import *
import dblas

for N in [5,100, 1000]:
    a = reshape(arange(float(N*N)), (N,N)) + identity(N)
    # print "typecode of array is: ", a.typecode()
    start = time.time()
    inv_a = inverse(a)
    stop = time.time()

    # print inv_a

    startmult = time.time()
    b =  matrixmultiply(a, inv_a)
    stopmult = time.time()
    residual = b - identity(N)

    #DGEMM(TRANSA, TRANSB, M, N, K, ALPHA,A, LDA, B, LDB, BETA,C,LDC)
    C = zeros(a.shape, a.typecode())

    alpha = array([1.0])
    beta  = array([1.0])
    startdgemm = time.time()
    dblas.dgemm('N','N', N,N,N, alpha,
                 a,N,       # A, LDA
                 inv_a, N,  # B, LDB
                 beta,  # BETA
                 C, N)
    stopdgemm = time.time()
    # print C

    print "%d x %d matrix inversion" % (N,N)
    print "max error is:",
    print maximum.reduce(fabs(ravel(residual)))
    print "time elapsed for inverse (sec) %f" % (stop-start)
    print "time elapsed for matrixmultiply (sec) %f" % (stopmult-startmult)
    print "time elapsed for dgemm-matrixmultiply (sec) %f" %
(stopdgemm-startdgemm)
    print
_______________________________________________
Numpy-discussion mailing list
Num...@li...
http://lists.sourceforge.net/mailman/listinfo/numpy-discussion

[Numpy-discussion] matrixmultiply implementation

[Christopher L. <cl...@gn...>]

I'm curious to see if anyone has looked at the implementation of
dot/matrixmultiply recently.  I was doing a benchmark lapack-based matrix
inversion and was surprised to find that the accuracy check was taking
longer than the inversion.  The bottleneck was the matrix multiply.  It
appears that a blas dgemm implementation might give a large improvement.
Here is a sample benchmark using PyLapack's dblas.dgemm for comparison.
Results for a 1000 x 1000 matrix inversion:

  time elapsed for inverse (sec) 10.050768
  time elapsed for matrixmultiply (sec) 24.142714
  time elapsed for dgemm-matrixmultiply (sec) 5.997188
  max error is: 3.30398486348e-10

Given the "dot" code doesn't look too bad, so I might be able to add dgemm
support myself, at least for some cases.

-chris

p.s. I'll include the test code below
#####################################
import time
from Numeric import *
from LinearAlgebra import *
import dblas

for N in [5,100, 1000]:
    a = reshape(arange(float(N*N)), (N,N)) + identity(N)
    # print "typecode of array is: ", a.typecode()
    start = time.time()
    inv_a = inverse(a)
    stop = time.time()

    # print inv_a

    startmult = time.time()
    b =  matrixmultiply(a, inv_a)
    stopmult = time.time()
    residual = b - identity(N)

    #DGEMM(TRANSA, TRANSB, M, N, K, ALPHA,A, LDA, B, LDB, BETA,C,LDC) 
    C = zeros(a.shape, a.typecode())

    alpha = array([1.0])
    beta  = array([1.0])
    startdgemm = time.time()
    dblas.dgemm('N','N', N,N,N, alpha, 
                 a,N,       # A, LDA
                 inv_a, N,  # B, LDB
                 beta,  # BETA
                 C, N)
    stopdgemm = time.time()
    # print C

    print "%d x %d matrix inversion" % (N,N)
    print "max error is:",
    print maximum.reduce(fabs(ravel(residual)))
    print "time elapsed for inverse (sec) %f" % (stop-start)
    print "time elapsed for matrixmultiply (sec) %f" % (stopmult-startmult)
    print "time elapsed for dgemm-matrixmultiply (sec) %f" % (stopdgemm-startdgemm)
    print

[Numpy-discussion] RE: [numpy - Help] Numeric-17.1.2.tar.gz missing

[Paul F. D. <pau...@ho...>]

I succeeded in getting the file at 2 p.m. Friday Dec 1 Pacific Standard
Time.

I had received a similar report earlier in the day so I don't think you're
crazy, just that something at SF was down.

-----Original Message-----
From: no...@so... [mailto:no...@so...]
Sent: Thursday, November 30, 2000 9:23 AM
To: no...@so...
Subject: [numpy - Help] Numeric-17.1.2.tar.gz missing


Read and respond to this message at:
http://sourceforge.net/forum/message.php?msg_id=80631
By: vanandel

There are no Numeric 17 tar files on http://download.sourceforge.net/numpy/

Could someone replace them?

Thanks!

______________________________________________________________________
You are receiving this email because you elected to monitor this forum.
To stop monitoring this forum, login to SourceForge and visit:
http://sourceforge.net/forum/monitor.php?forum_id=3848

Re: [Numpy-discussion] (no subject)

[<ro...@ho...>]

>>>>> "s" == szport  <sz...@ru...> writes:


 s> Is the following definition faster or not?

 s> def ones(shape, typecode=\'l\', savespace=0):
 s>    a = zeros( (product(shape),), typecode, savespace)
 s>    a[:] = 1
 s>    a.shape = shape
 s>    return a

This is yet another way to write a[...]=1 manually via a small
detour. For small arrays (size 10...1000) it is about twice slower
than writing [...] (0.1ms slower per call on my workstation).

Rob

-- 
=====   ro...@ho...          http://www.hooft.net/people/rob/  =====
=====   R&D, Nonius BV, Delft  http://www.nonius.nl/             =====
===== PGPid 0xFA19277D ========================== Use Linux! =========

[Numpy-discussion] (no subject)

[<sz...@ru...>]

> I had another look at the definition of \"ones\" and of
another routine
>           > I frequently use: arange. It appears that even
without rewriting them
>           > in C, some speedup can be achieved:
>           > 
>           >  - in ones(), the + 1 should be done \"in
place\", saving about 15%, more
>           >    if you run out of processor cache:
>
>           I\'d also try assignment in place:
>
>           def ones(shape, typecode=\'l\', savespace=0):
>              a = zeros(shape, typecode, savespace)
>              a[len(shape)*[slice(0, None)]] = 1
>              return a
>
>           Konrad.

Is the following definition faster or not?

def ones(shape, typecode=\'l\', savespace=0):
    a = zeros( (product(shape),), typecode, savespace)
    a[:] = 1
    a.shape = shape
    return a

Zaur

Re: [Numpy-discussion] Initialization of array?

[<ro...@ho...>]

 KH> I had another look at the definition of "ones" and of another routine
 KH> I frequently use: arange. It appears that even without rewriting them
 KH> in C, some speedup can be achieved:
 KH> 
 KH> - in ones(), the + 1 should be done "in place", saving about 15%, more
 KH> if you run out of processor cache:

 KH> I'd also try assignment in place:

 KH> def ones(shape, typecode='l', savespace=0):
 KH>    a = zeros(shape, typecode, savespace)
 KH>    a[len(shape)*[slice(0, None)]] = 1
 KH>    return a

This is even faster, but it is better to write "a[...] = 1", because
your manual calculation of "..." gives a large overhead for small arrays.

On another machine this time:
    Numeric.ones             10 ->    0.254ms
    Numeric.ones            100 ->    0.268ms
    Numeric.ones           1000 ->    0.340ms
    Numeric.ones          10000 ->    1.960ms
    Numeric.ones         100000 ->   29.300ms
    Numeric.zeros            10 ->    0.055ms
    Numeric.zeros           100 ->    0.059ms
    Numeric.zeros          1000 ->    0.068ms
    Numeric.zeros         10000 ->    0.430ms
    Numeric.zeros        100000 ->    9.800ms
    Add inplace              10 ->    0.246ms
    Add inplace             100 ->    0.255ms
    Add inplace            1000 ->    0.312ms
    Add inplace           10000 ->    1.270ms
    Add inplace          100000 ->   18.100ms
    Assign inplace           10 ->    0.192ms
    Assign inplace          100 ->    0.201ms
    Assign inplace         1000 ->    0.242ms
    Assign inplace        10000 ->    1.010ms
    Assign inplace       100000 ->   16.300ms
    Reshape 1                10 ->    0.842ms
    Reshape 1               100 ->    1.175ms
    Reshape 1              1000 ->    4.100ms
    Reshape 1             10000 ->   35.100ms
    Reshape 1            100000 ->  368.600ms

Rob

-- 
=====   ro...@ho...          http://www.hooft.net/people/rob/  =====
=====   R&D, Nonius BV, Delft  http://www.nonius.nl/             =====
===== PGPid 0xFA19277D ========================== Use Linux! =========