Mailman 3 March 2002 - NumPy-Discussion

[Q] RandomArray module has a problem ????
by newton＠ns.kias.re.kr 05 Mar '02

05 Mar '02

Hi, RandomArray has a problem? I use python2.0, Numpy 20.3. simple source code is ... -------------------------------- #!/usr/bin/env python import math import Numeric import RandomArray import sys RandomArray.seed(1234,5678) i=0L while 1: i = i+1 a = RandomArray.randint(0,100) if a==100: print 'i=',i, 'a=',a --------------------------------- and result is --------------------------------- i= 70164640 a= 100 i= 152242967 a= 100 i= 159619195 a= 100 i= 173219763 a= 100 i= 200933959 a= 100 i= 233219191 a= 100 i= 276114822 a= 100 i= 313589319 a= 100 i= 340689813 a= 100 i= 402397265 a= 100 i= 456099215 a= 100 i= 506078935 a= 100 i= 547758957 a= 100 i= 559163554 a= 100 i= 570211180 a= 100 .......... --------------------------------- RandomArray.randint(0,100) has range 0<= RandomArray.randint(0,100) <100. But, result is not...somtime, a==100 arise. So, I upgrade to python 2.2 and Numpy 21b3. But, I met same problem. And, so, I change the os Mandrake 8.0 to Redhat 7.2. But, same problem... I don't know what is my mistake... Please help me ... Kee-Hyoung Joo ------------------------------------------------------------------ I love Jesus Christ who is my savior. He gives me meanning of life. In Christ, I have become shepherd and bible teacher. e-mail : newton(a)kias.re.kr home : http://newton.skku.ac.kr/~newton (My old home page) ------------------------------------------------------------------

1 0

[Q] RandomArray module has a problem ????
by newton＠ns.kias.re.kr 05 Mar '02

05 Mar '02

Hi, RandomArray has a problem? I use python2.0, Numpy 20.3. simple source code is ... -------------------------------- #!/usr/bin/env python import math import Numeric import RandomArray import sys RandomArray.seed(1234,5678) i=0L while 1: i = i+1 a = RandomArray.randint(0,100) if a==100: print 'i=',i, 'a=',a --------------------------------- and result is --------------------------------- i= 70164640 a= 100 i= 152242967 a= 100 i= 159619195 a= 100 i= 173219763 a= 100 i= 200933959 a= 100 i= 233219191 a= 100 i= 276114822 a= 100 i= 313589319 a= 100 i= 340689813 a= 100 i= 402397265 a= 100 i= 456099215 a= 100 i= 506078935 a= 100 i= 547758957 a= 100 i= 559163554 a= 100 i= 570211180 a= 100 .......... --------------------------------- RandomArray.randint(0,100) has range 0<= RandomArray.randint(0,100) <100. But, result is not...somtime, a==100 arise. So, I upgrade to python 2.2 and Numpy 21b3. But, I met same problem. And, so, I change the os Mandrake 8.0 to Redhat 7.2. But, same problem... I don't know what is my mistake... Please help me ... Kee-Hyoung Joo ------------------------------------------------------------------ I love Jesus Christ who is my savior. He gives me meanning of life. In Christ, I have become shepherd and bible teacher. e-mail : newton(a)kias.re.kr home : http://newton.skku.ac.kr/~newton (My old home page) ------------------------------------------------------------------

1 0

[Q] RandomArray has a problem ??
by 주기형 05 Mar '02

05 Mar '02

Hi, RandomArray has a problem? I use python2.0, Numpy 20.3. simple source code is ... -------------------------------- #!/usr/bin/env python import math import Numeric import RandomArray import sys RandomArray.seed(1234,5678) i=0L while 1: i = i+1 a = RandomArray.randint(0,100) if a==100: print 'i=',i, 'a=',a --------------------------------- and result is --------------------------------- i= 70164640 a= 100 i= 152242967 a= 100 i= 159619195 a= 100 i= 173219763 a= 100 i= 200933959 a= 100 i= 233219191 a= 100 i= 276114822 a= 100 i= 313589319 a= 100 i= 340689813 a= 100 i= 402397265 a= 100 i= 456099215 a= 100 i= 506078935 a= 100 i= 547758957 a= 100 i= 559163554 a= 100 i= 570211180 a= 100 .......... --------------------------------- RandomArray.randint(0,100) has range 0<= RandomArray.randint(0,100) <100. But, result is not...somtime, a==100 arise. So, I upgrade to python 2.2 and Numpy 21b3. But, I met same problem. And, so, I change the os Mandrake 8.0 to Redhat 7.2. But, same problem... I don't know what is my mistake... Please help me ... Kee-Hyoung Joo

1 0

RE: [Numpy-discussion] RE: Python 2.2 seriously crippled for numerical computation?
by Tim Peters 05 Mar '02

05 Mar '02

[Huaiyu Zhu] > ... > 1. The -lieee is indeed the most direct cure. On the specific platform tried. The libm errno rules changed between C89 and C99, and I'm afraid there's no errno behavior Python can rely on anymore. So I expect more changes will be needed in Python, regardless of how things turn out on this specific platform. > ... > 2. Is there a configure option to guarantee -lieee? If anyone can answer this question, please don't answer it here: it will just get lost. Attach it to Huaiyu's bug report instead: <http://sf.net/tracker/?func=detail&aid=525705&group_id=5470&atid=105470> Thanks. > ... > 3. errno 34 (or -lieee) may not be the sole reason. > > On a RedHat 6.1 upgraded to 7.1 (both gcc and glibc), errno 34 > is indeed raised in a C program linked without -lieee, and Python is > indeed compiled without -lieee, but Python does not raise > OverflowError. I expect you're missing something. Skip posted the Python code before, and if errno gets set, Python *does* raise OverflowError: errno = 0; /* Skip forgot to post this line, and it's important */ ... ix = pow(iv, iw); ... if (errno != 0) { /* XXX could it be another type of error? */ PyErr_SetFromErrno(PyExc_OverflowError); return NULL; If you can read C, believe your eyes <wink>. What you may be missing is what an utter mess C is here. How libm behaves may depend on compiler options, linker options, global system flags, and even options set for other libraries you link with. > ... > 4. Is there an easier way to debug such problems? The cause was obvious to the first person (Skip) who stepped into Python to see what the code did on a platform where it failed. It's not going to be obvious to someone who doesn't. > 5. How is 1e200**2 handled? It goes through exactly the same code. > Since both 1e-200**2 and 1e200**2 produce the same errno all the time, > but Python still raises OverflowError for 1e200**2 when linked with > -lieee, there must be a separate mechanism at work. You're speculating from a false base: if platform pow(x, y) sets errno to any non-zero value, Python x**y raises OverflowError. What differs is when platform pow(x, y) does not set errno. In that case, Python synthesizes errno=ERANGE if the pow() result equals +- platform HUGE_VAL. > What is that and how can I override it? Sorry, you can't override it. > I know this is by design, but I think the design is dumb (to put it > politely). I won't get into an argument here. I'll write > up my rationale against this when I have some time. I'm afraid a rationale won't do any good. I'm in favor of supplying full 754 compatibility myself, but: A) Getting from here to there requires volunteers to design, implement, document, and test the code. Given the atrocious state of C's 754 story, and the subtlety of 754's requirements, this needs volunteers who are both 754 experts and platform C experts. That combination is rare. B) Python's floating-point users will never agree it's a good thing, so such a change requires careful backward compatibility work too. This isn't likely to get done by someone who characterizes the other side as "dumb (to put it politely)" <0.9 wink>. Note that the only significant floating-point code ever contributed to the Python core was the fpectl module, and its purpose is to *break* 754 "non-stop" exception semantics in the places Python happens to let them sneak through. > I do remember there being several discussions in the past, but I don't > remember any convincing argument for the current decision. Any URL > would be greatly appreciated, beside the one pointed by Tim. Which "current decision" do you have in mind? There is no design doc for Python's numerics, if that's what you're looking for. As the text at the URL I gave you said, much of Python's fp behavior is accidental, inherited from platform C quirks.

1 0

Interface for numpy C-API <-> simple C++ matrix class
by Stefan Heinrichs 04 Mar '02

04 Mar '02

Hello, we would like to write routines processing numerical python arrays in C++, so that at least boundary checking can be enabled at runtime. While there are a lot of matrix libraries available for C++, I could not find the glue that interfaces such a library to the C-API of numerical python. Seamless access to a minimal C++ library would make the C++ part of programming much easier. Has anyone already written some wrapper/glue code? Thanks and best regards, Stefan -- ------------------------------------------------------------------- Email: Stefan.Heinrichs(a)uni-konstanz.de Address: Fakulaet fuer Physik, Universitaet Konstanz, Universitaetsstr.10, 78457 Konstanz, Germany Phone: +49 7531 88 3814

4 4

RE: Matrix.py change breaks LinearAlgebra
by Paul F Dubois 04 Mar '02

04 Mar '02

Travis fixed the error he accidentally made when improving MLab.py. This incident points out that our test suite does not include any tests for the non-core modules. Recently the discussion over the meaning of cov shows this too: I didn't even have a test that showed what the INTENDED answer is. We need test suites for FFT, LinearAlgebra, MLab, etc. If you have the subject competence to make a test file for us, please volunteer. I'd like them to be like the ones in Test/test.py, but as a separate file so that I can test the modules separately.

1 0

Matrix.py change breaks LinearAlgebra
by Paul F Dubois 04 Mar '02

04 Mar '02

Your change to Matrix.py has a fatal flaw: Python 2.2 (#1, Mar 1 2002, 11:11:28) [GCC 2.96 20000731 (Red Hat Linux 7.1 2.96-81)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import LinearAlgebra Traceback (most recent call last): File "<stdin>", line 1, in ? File "/home/dubois/linux/lib/python2.2/site-packages/Numeric/LinearAlgebra.py ", line 10, in ? import MLab File "/home/dubois/linux/lib/python2.2/site-packages/Numeric/MLab.py", line 10, in ? import Matrix File "/home/dubois/linux/lib/python2.2/site-packages/Numeric/Matrix.py", line 8, in ? from LinearAlgebra import inverse ImportError: cannot import name inverse >>>

1 0

How can CDOUBLE_to_CDOUBLE work correctly?
by Pearu Peterson 03 Mar '02

03 Mar '02

Hi! I am trying to copy a 2-d complex array to another 2-d complex array in an extension module. Both arrays may be noncontiguous. But using the same routine (namely Travis's copy_ND_array, you can find it at the end of this messsage) as for real arrays seems not work. After some playing around and reading docs about strides and stuff, I found that the reason might be in how Numeric (20.3) defines the CDOUBLE_to_CDOUBLE function: static void CDOUBLE_to_CDOUBLE(double *ip, int ipstep, double *op, int opstep, int n) {int i; for(i=0;i<2*n;i++,ip+=ipstep,op+=opstep) {*op = (double)*ip;}} It seems not to take into account that real and imaginary part are always contiguous in memory, even if an array itself is not. Actually, I don't understand how this code can work (unless some magic is done in places where this code is used). I would have expected that the code for CDOUBLE_to_CDOUBLE to be analoguous to relative functions but for the real data. For example, DOUBLE_to_DOUBLE is defined as static void DOUBLE_to_DOUBLE(double *ip, int ipstep, double *op, int opstep, int n) {int i; for(i=0;i<n;i++,ip+=ipstep,op+=opstep) {*op = (double)*ip;}} So, I would have expected CDOUBLE_to_CDOUBLE to be static void CDOUBLE_to_CDOUBLE(double *ip, int ipstep, double *op, int opstep, int n) { int i; for(i=0;i<n;i++,ip+=ipstep,op+=opstep) { *op = (double)*ip; /* copy real part */ *(op+1) = (double)*(ip+1); /* copy imaginary part that always follows the real part in memory */ } } Could someone explain how Numeric can work with the current CDOUBLE_to_CDOUBLE? Because I don't understand which one is broken, the Numeric's CDOUBLE_to_CDOUBLE (and relative) functions, or my code. The latter may be the case, but to fix it, I need some clarification on the issue. Can you help me? Thanks, Pearu /************************* copy_ND_array *******************************/ #define INCREMENT(ret_ind, nd, max_ind) \ { \ int k; \ k = (nd) - 1; \ if (k<0) (ret_ind)[0] = (max_ind)[0]; else \ if (++(ret_ind)[k] >= (max_ind)[k]) { \ while (k >= 0 && ((ret_ind)[k] >= (max_ind)[k]-1)) \ (ret_ind)[k--] = 0; \ if (k >= 0) (ret_ind)[k]++; \ else (ret_ind)[0] = (max_ind)[0]; \ } \ } #define CALCINDEX(indx, nd_index, strides, ndim) \ { \ int i; \ indx = 0; \ for (i=0; i < (ndim); i++) \ indx += nd_index[i]*strides[i]; \ } extern int copy_ND_array(const PyArrayObject *in, PyArrayObject *out) { /* This routine copies an N-D array in to an N-D array out where both can be discontiguous. An appropriate (raw) cast is made on the data. */ /* It works by using an N-1 length vector to hold the N-1 first indices into the array. This counter is looped through copying (and casting) the entire last dimension at a time. */ int *nd_index, indx1; int indx2, last_dim; int instep, outstep; if (0 == in->nd) { in->descr->cast[out->descr->type_num]((void *)in->data,1, (void*)out->data,1,1); return 0; } if (1 == in->nd) { in->descr->cast[out->descr->type_num]((void *)in->data,1, (void*)out->data,1,in->dimensions[0]); return 0; } nd_index = (int *)calloc(in->nd-1,sizeof(int)); last_dim = in->nd - 1; instep = in->strides[last_dim] / in->descr->elsize; outstep = out->strides[last_dim] / out->descr->elsize; if (NULL == nd_index ) { fprintf(stderr,"Could not allocate memory for index array.\n"); return -1; } while(nd_index[0] != in->dimensions[0]) { CALCINDEX(indx1,nd_index,in->strides,in->nd-1); CALCINDEX(indx2,nd_index,out->strides,out->nd-1); /* Copy (with an appropriate cast) the last dimension of the array */ (in->descr->cast[out->descr->type_num])((void*)(in->data+indx1),instep, (void*)(out->data+indx2),outstep,in->dimensions[last_dim]); INCREMENT(nd_index,in->nd-1,in->dimensions); } free(nd_index); return 0; } /* EOF copy_ND_array */

1 1

Re: [Numpy-discussion] RE: Python 2.2 seriously crippled for numerical computation?
by Andrew MacIntyre 02 Mar '02

02 Mar '02

On 2 Mar 2002, Konrad Hinsen wrote: > Tim Peters <tim.one(a)comcast.net> writes: > > > > # Python 2.2 > > > > > > >>> 1e-200**2 > > > Traceback (most recent call last): > > > File "<stdin>", line 1, in ? > > > OverflowError: (34, 'Numerical result out of range') > > > > That one is surprising and definitely not intended: it suggests your > > platform libm is setting errno to ERANGE for pow(1e-200, 2.0), or that your > > platform C headers define INFINITY but incorrectly, or that your platform C > > headers define HUGE_VAL but incorrectly, or that your platform C compiler > > generates bad code, or optimizes incorrectly, for negating and/or comparing > > I just tested and found the same behaviour, on RedHat Linux 7.1 > running on a Pentium machine. Python 2.1, compiled and running on the > same machine, returns 0. So does the Python 1.5.2 that comes with the > RedHat installation. Although there might certainly be something wrong > with the C compiler and/or header files, something has likely changed > in Python as well in going to 2.2, the only other explanation I see > would be a compiler optimization bug that didn't have an effect with > earlier Python releases. Other examples... FreeBSD 4.4: Python 2.1.1 (#1, Sep 13 2001, 18:12:15) [GCC 2.95.3 20010315 (release) [FreeBSD]] on freebsd4 Type "copyright", "credits" or "license" for more information. >>> 1e-200**2 0.0 >>> 1e200**2 Inf Python 2.3a0 (#1, Mar 1 2002, 00:00:52) [GCC 2.95.3 20010315 (release) [FreeBSD]] on freebsd4 Type "help", "copyright", "credits" or "license" for more information. >>> 1e-200**2 0.0 >>> 1e200**2 Traceback (most recent call last): File "<stdin>", line 1, in ? OverflowError: (34, 'Result too large') Both builds built with "./configure --with-fpectl", although the optimisation settings for the 2.3a0 (CVS of yesterday) build were tweaked (2.1.1: "-g -O3"; 2.3a0: "-s -m486 -Os"). My 2.2 OS/2 EMX build (which uses gcc 2.8.1 -O2) produces exactly the same result as 2.3a0 on FreeBSD. -- Andrew I MacIntyre "These thoughts are mine alone..." E-mail: andymac(a)bullseye.apana.org.au | Snail: PO Box 370 andymac(a)pcug.org.au | Belconnen ACT 2616 Web: http://www.andymac.org/ | Australia

1 0

numarray interface and performance issues (for dot product and transpose)
by a.schmolck＠gmx.net 01 Mar '02

01 Mar '02

Hi, Numeric is an impressively powerful and in many respects easy and comfortable to use package (e.g. it's sophisticated slicing operations, not to mention the power and elegance of the underlying python language) and one would hope that it can one day replace Matlab (which is both expensive and a nightmare as a programming language) as a standard platform for numerical calculations. Two important drawbacks of Numeric (python) compared to Matlab currently seem to be the plotting functionality and the availability of specialist libraries (say for pattern recognition problems etc.). However python seems to be quickly catching up in this area (e.g. scipy and other projects are bustling with activity and, as far as I know, there is currently a promising undergraduate project on the way at my university to provide powerful plotting facilities for python). There is however a problem that, for the use to which I want to put Numeric, runs deeper and provides me with quite a headache: Two essential matrix operations (matrix-multiplication and transposition (which is what I am mainly using) are both considerably a) less efficient and b) less notationally elegant under Numeric than under Matlab. Before I expound on that, I'd better mention two things: Firstly, I, have to confess largely ignorant about the gritty details of implementing efficient numerical computation (and indeed, I can't even claim to possess a good grasp of the underlying Linear Algebra). Secondly, and partly as a consequence, I realize that I'm likely to view this from a rather narrow perspective. Nonetheless, I need to solve these problems at least for myself, even if such a solution might be inappropriate for a more general audience, so I'd be very keen to hear what suggestions or experiences from other Numeric users are. I also want to share what my phd supervisor, Richard Everson, and I have come up with so far in terms of a), as I think it might be quite useful for other people, too (see [2], more about this later). Ok, so now let me proceed to state clearly what I mean by a) and b): The following Matlab fragment M * (C' * C) * V' * u currently has to be written like this in python: dot(dot(dot(M, dot(transpose(C), C)), transpose(v)), u) Or, even worse if one doesn't want to pollute the namespace: Numeric.dot(Numeric.dot(Numeric.dot(Numeric.M, Numeric.dot(Numeric.transpose(C), C)), Numeric.transpose(v)), u) Part a) (efficiency) -------------------- Apart from the syntactic inconveniences, the Matlab above will execute __considerably__ faster if C is a rather larger matrix (say 1000x1000). There are AFAIK two chief reasons for this: 1. Matlab uses a highly optimized ATLAS[3] BLAS routine to compute the dot-product, whereas Numeric uses a straight-forward home-brewn solution. 2. Numeric performs unnecessary transpose operations (prior to 20.3, I think, more about this later). The transpose operation is really damaging with big matrices, because it creates a complete copy, rather than trying to do something lazy (if your memory is already almost half filled up with (matrix) C, then creating a (in principle superfluous) transposed copy is not going to do you any good). The above C' * C actually creates, AFAIK, _3_ versions of C, 2 of them transposed (prior to 20.3; dot(a,b) translates into innerproduct(a, swapaxes(b, -1, -2)) In newer versions of Numeric, this is replaced by multiarray.matrixproduct(a, b) which has the considerable advantage that it doesn't create an unnecessary copy and the considerable disadvantage that it seems to be factor 3 or so slower than the old (already not blazingly fast) version for large Matrix x Matrix multiplication, (see timing results [1])). Now fortunately, we made some significant progress on the performance issues. My supervisor valiantly patched Numeric's multiarrayobject.c to use ATLAS for {scalar,vector,matrix} * {scalar,vector,matrix} multiplications. As a consequence the multiplication of pair of 1000x1000 matrices is computed more than _40_ times faster on my athlon machine (see [1], it seems to work fine but both the timing results as well as the code should be treated with a bit of caution). The BLAS-enabled dot takes almost exactly the same time as Matlab doing the same thing. In my case, this makes the difference between waiting several minutes and a day, and between deciding to use python or something faster but more painful. Although to benefit from this patch, one has to install ATLAS this is straight-forward (see [2] and [3] for pointers). In addition it is easy to install Lapack at the same time, as ATLAS provides significantly optimised versions of some of the Lapack routines too. Using Lapack rather than the lapack_lite that comes with Numeric also has the advantage that it is extensively tested and works very reliably (e.g. the lapack_lite routine that is called by Heigenvalues is pretty broken and searching the net suggests that this has been the case for some time now). This still leaves point 2), transposes, as a potential source for improvement. As far as I understand, the Lapack routines operate on flattened arrays, i.e. the only difference on calling some routine on a matrix and its transpose is that one has to pass a different additional parameters in the second case that specifies the structure of the matrix. Therefore if one were to program C' * C in fortran, only C itself needed to exist in memory. Would it therefore be possible that operations like transpose, slicing etc. be carried out lazily (i.e. so that only _modifying_ the result of the transpose or slice or passing it to some function for some reason needs a real copy resulted in an actual copy operation being performed)? As far as I understand the slicing in numarray is no longer simple aliasing (which is a good thing, since it make arrays more compatible with lists), are there any plans for implementing such an on-demand scheme? Part b) (syntax) ________________ As I said, dot(dot(dot(M, dot(transpose(C), C)), transpose(v)), u) is pretty obscure compared to M * (C' * C) * V' * u) Although linear algebra is not the be all and end all of numerical calculations, many numerical algorithms end up involving a good deal of linear algebra. The awkwardness of python/Numeric's notation for common linear algebra operation severely impairs, in my opinion, the utility of writing, using and debugging algorithms in python/Numeric. This is a particular shame given python's generally elegant and expresssive syntax. So far, I've thought of the following possible solutions: 0. Do nothing: Just live with the awkward syntax. 1. Comment: add a comment like the Matlab line above above complex expressions. 2. Wait for fix. Hope that either: 2.1. numarray will overload '*' as matrix multiplication, rather than elementwise multiplication (highly unlikely, I suspect). 2.2. the core python team accepts one of the proposed operator additions (or even, fond phantasy, allows users to add new operators) 3. Wrap: create a DotMatrix class that overloads '*' to be dot and maybe self.t to return the transpose -- this also means that all the numerical libraries I frequently use need to be wrapped. 4. Preprocess: Write some preprocessor that translates e.g. A ~* B in dot(A,B). 5. Hack: create own customized version of array that overloads e.g. '<<' to mean dot-product. The possible downsides: 0. Do nothing: Seems unacceptable (unless as a temporary solution), because the code becomes much less readable (and consequently harder to understand and more error-prone), thus forgoing one of the main reasons to choose python in the first place. 1. Comment: Unsatisfactory, because there is no way to gurantee that comment and code are in synch, which very likely will lead to difficult to find bugs. 2. Wait for a fix: Either would be nice (I could live with having to write multiply(a,b) -- I suppose however some other people can't, because I can't think of another reason why Numeric didn't overload * for matrixproduct in the first place, moreover it would mean a significant interface change (and code breakage), so I guess it is rather unlikely). From what I gather from searching the net and looking at peps there is also not much of a chance to get a new operator anytime soon. 3. Wrap: Promises to be a maintenance nightmare (BTW: will the new array class be subclassable?), but otherwise looks feasible. Has anyone done this? 4. Preprocess: Would have the advantage that I could always get back to "standard" python/Numeric code, but would involve quite a bit of work and might also break tools that parse python. 5. Hack array: Seems evil, but in principle feasible, because AFAIK '<<' isn't used in Numeric itself and hopefully it wouldn't involve too much work. However, given a free choice '<<' is hardly the operator I would choose to represent the dot product. I am not completely sure what the policy and rationale behind the current division of labor between functions, methods and attributes in numarray is but as far as the lack of a missing transposition postfix operator is concerned, one reasonable approach to make the transpose operation more readable (without the need to change anything in python itself) would seem to me to provide arrays with an attribute .t or .T so that: a.t == transpose(a) While this would admittedly be a bit of a syntactic hack, the operation is so commonplace and the gain in readability (in my eyes) is so significant that it would seem to me to merit serious consideration (maybe a.cc or so might also be an idea for complex conjugate notation, but I'm less convinced about that?). If the addition of a single operator for the exclusive benefit of Numeric users is rejected by the core python team, maybe it's worthwhile lobbying for some mechanism that allows users to define new operators (like e.g. in Haskell)... OK, that's all -- thanks for bearing with me through this long email. Suggestions and comments about the patch [2] and possible solutions to issues raised are greatly appreciated, regards, alexander schmolck Footnotes: [1] timing results for the patched version of Numeric[2] comparing new and old 'dot' performance: http://www.dcs.ex.ac.uk/~aschmolc/Numeric/TimingsForAtlasDot.txt [2] A patch for Numeric 21.1b that speeds up Numeric's dot function by several factors for large matrices can be found here: http://www.dcs.ex.ac.uk/~aschmolc/Numeric/ [3] ATLAS (http://math-atlas.sourceforge.net/) is a project to provide a highly optimized version of BLAS (Basic Linear Algebra Subroutines, a standard and thouroughly tested implementation of basic linear algebra operations) and some LAPACK routines. The charm of ATLAS is that it is platform-independent and yet highly optimized, which it achieves by essentially fine tuning a number of parameters until optimum performance for the _particular_ machine on which it is built is reached. As a consequence complete builds can take some time, but binary versions for ATLAS for common processors are available from http://www.netlib.org/atlas/archives (moreover, even if one decides to build ATLAS oneself, the search space can be considerably cut down if one accepts the suggested "experience" values offered during the make process). -- Alexander Schmolck Postgraduate Research Student Department of Computer Science University of Exeter A.Schmolck(a)gmx.net http://www.dcs.ex.ac.uk/people/aschmolc/

7 6