Changes to improve performance on small matricies
We are porting some code from Matlab that does many thousands of operations on small matrices. We have found that some small changes to LinearAlgebra.py, Numeric.py, and multiarraymodule.c make our code run approximately twice as fast as with the unmodified Numeric-20.2 code. Our changes specifically replace LinearAlgebra._castCopyAndTranspose with a version we call _fastCopyAndTranspose that does most of the work in C instead of calling Numeric.transpose, which calls arange, which is quite expensive. We also optimized Numeric.dot to call a new function multiarray.matrixproduct which does the axis swap on the fly instead of calling swapaxes (which calls arange) and then calling innerproduct (which as the very first step copies the transposed matrix to make it contiguous). Formerly our code spent much of its time in arange (which we never explicitly call), dot, and _castCopyAndTranspose. The changes described above eliminate this overhead. I'm writing to ask if these changes might make a worthy patch to NumPy? We have tested them with on Windows2k (both Native and under Cygwin) and on Linux. Soon, we'll have a test on Mac OS X.1. If anyone is interested, I will figure out how to generate a patch file to submit. Thanks gb
Gary Bishop wrote:
I'm writing to ask if these changes might make a worthy patch to NumPy? We have tested them with on Windows2k (both Native and under Cygwin) and on Linux. Soon, we'll have a test on Mac OS X.1.
I'd say yes, but we do have to be satisfied with the robustness of the code, and I have no idea how to go about doing that. This brings up a good point. In Numeric.Py, the is the following comment: #The following functions are considered builtin, they all might be #in C some day I'd love to see that happen. As evidenced by this poster's example, Numeric could be made a lot faster by optimizing a number of build-in functions. What would it take to get a "Numeric Optimization project" going?? A few ideas: unit tests: If we had unit tests designed for all the built-in functions, we would have a good way to know if a newly coded C version is working right. How-To_Docs: I have tried to write some Numeric extensions myself, and it is very difficult. Writing one that does something unusual for my particular application is not so hard with the current docs. In that case it usually only has to work for one type, I know what shape the array should be etc. Writing a universal extension to Numeric that would work on any Numeric Array is a whole lot harder!! I am currently working on two of these, and the hard parts for me are: Doing something to all the data in a perhaps discontiguous array of arbitrary shape. I have soemthing working now, but it's pretty slow. Having a function work for all types of Arrays: dealing this this kin dof dynamicism on C is very difficult for me. The best place to look is clearly the source itself, but I have tried, and it's pretty darn dense. A How-To would be pretty darn helpful. I don't know who's going to write these things, but I think you could enlarge the pool of NumPy developers quite a bit. If anyone is interested, I have written "fastclip()" which does what clip() does, but it does it in place, and only works on Floats (see above) It is very fast for contiguous arrays, and disappointingly slow on discontiguous arrays (also see above). Utill I get it working on all types, there isn't much point in submitting it for inclusion in NumPy. The other two I have written I use for dealing with binary data coming in from files: byteswap() does a bytswap in place, whioch is faster an more memory efficient that the existing bytswapped() method. changetype() changes the type of the array, using the same data: changetype(m,typecode) = fromstring(m.tostring(),typecode) It is a lot faster and memory efficient (the memory was the issue for me, I'm dealing with up to 10MB arrays), but kind of an obscure need to, so this may not be a candidate for inclusion in NumPy either. I still have some questions about my code for that one: I'll post a separate message about that. -Chris -- Christopher Barker, Ph.D. ChrisHBarker@home.net --- --- --- http://members.home.net/barkerlohmann ---@@ -----@@ -----@@ ------@@@ ------@@@ ------@@@ Oil Spill Modeling ------ @ ------ @ ------ @ Water Resources Engineering ------- --------- -------- Coastal and Fluvial Hydrodynamics -------------------------------------- ------------------------------------------------------------------------
HI all, I've written a function (in C) that changes the type of an array, in place, using the same data block. The reason I need this is because I am reading in a large block of data from a file that is essentially a set of records, each one a C struct that has a couple of longs and a char in it. I read it in as a big array, slice iot to separate the fields, and then need to change the type. I had been using: A = fromstring(A.tostring(),Int32) This works fine, but needs a lot of memory, and it kind of slow (althought enough faster than the file read that it matters little) The array can be on order of 10MB however, so the memory issue is a problem. Anyway, I wrote a function that changes the type of data in an array, without making any copies of the data. It seems to work but I have a few questions (code to follow): 1) Is there an easier way to do this? I just notices that PyArray_As1D takes a type as an input... would it do what I want?? 2) Is there a way , at the C levbel to reshape an array. PyArray_Reshape takes a PyObject, which is kind of awkward, I'd love to pass in a few integers, or maybe a nd, and a pointer to a dimensions array. 3) What I did do is use PyArray_DescrFromType() to create a new desciption, and then assign the description pointer of the array to the new one. This seems to work, but I expect that the old description memory never gets freed. How do I free it (can you tell I'm a C newbie?). Note that I have called this function is a loop thousands of times, and not seen memory use go up, but it's a pretty small structure. 4) If I call PyArray_DescrFromType() with an invalid typcode, I get a crash. How can I check if the typecode is valid?? The code follows, I'd love any feedback anyone can offer. -Chris static PyObject * NumericExtras_changetype(PyObject *self, PyObject *args) { PyArrayObject *Array; //int *new_strides; //int *new_dimensions; PyArray_Descr *NewDescription; int TotalBytes; int NewElementSize; char typecode; if (!PyArg_ParseTuple(args, "O!c", &PyArray_Type, &Array, &typecode)) { return NULL; } if (!PyArray_ISCONTIGUOUS(Array)){ PyErr_SetString (PyExc_ValueError, "m must be contiguous"); return NULL; } TotalBytes = PyArray_NBYTES(Array); // This will crash if an invalid typecode is passed in. NewDescription = PyArray_DescrFromType(typecode); NewElementSize = NewDescription->elsize; if (TotalBytes % NewElementSize > 0){ PyErr_SetString (PyExc_ValueError, "The input array is the wrong size for the requested type"); return NULL; } Array->descr = NewDescription; // does the old descr need to be freed?? how?? Array->nd = 1; Array->strides[0] = NewElementSize; Array->dimensions[0] = TotalBytes / NewElementSize; return Py_BuildValue(""); } -- Christopher Barker, Ph.D. ChrisHBarker@home.net --- --- --- http://members.home.net/barkerlohmann ---@@ -----@@ -----@@ ------@@@ ------@@@ ------@@@ Oil Spill Modeling ------ @ ------ @ ------ @ Water Resources Engineering ------- --------- -------- Coastal and Fluvial Hydrodynamics -------------------------------------- ------------------------------------------------------------------------
participants (2)
-
Chris Barker -
Gary Bishop