Hi all,
I just wrote ctypes bindings to fftw3 (see
http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html
for the post to scipy).
Now I have a couple of numpy related questions:
In order to be able to use simd instructions I
create an ndarray subclass, which uses fftw_malloc to allocate the
memory and fftw_free to free the memory when the array is deleted. This
works fine for inplace operations however if someone does something like
this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an
even bigger problem when executing fftw plans), however
type(a) still gives me
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
.
This might help some: http://www.scipy.org/Subclasses Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma
On Mon, 2009-01-26 at 19:25 -0600, Ryan May wrote:
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
. This might help some:
http://www.scipy.org/Subclasses
Ryan
Thanks, I had read about __array_finalize__, but not about __array_wrap__. I'm not quite sure if I understand how I can use this to get around my problem, can you explain? Cheers Jochen
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
.
I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff... cheers, David
On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote:
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag:
http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff...
cheers,
David
Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan. You can then either execute the plan object, which executes the fftw_plan, or you can call the fftw guru method fftw_execute_dft(inarray, outarray), where I don't do any checking on the arrays, to keep the performance high. I agree that such an approach is not quite as intuitive as x=fft(y), but for my application (propagation equations) I usually perform a large number of FFTs on two arrays, so I just opted for staying close to the fftw way of doing things. So my problem is not really so much with the alignment, it's more with executing the plans, even if someone has created a plan from two unaligned arrays, if he does something like x=x+1 somewhere in the calculations and later executes a plan which was created using x, x will not change because it's memory is not pointing to the original location. I was looking for some way to prevent this. Cheers Jochen
______________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Jochen wrote:
On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote:
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag:
http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff...
cheers,
David
Hi David,
I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan.
I am not sure I follow you when you say the "fftw way": you can avoid having to store the arrays altogether while still using fftw plans, there is nothing "unfftw" about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated. David
On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote:
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag:
http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff...
cheers,
David
Hi David,
I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan.
I am not sure I follow you when you say the "fftw way": you can avoid having to store the arrays altogether while still using fftw plans, there is nothing "unfftw" about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated.
David
Sorry maybe I wasn't quite clear, what I mean by the "fftw way" is creating a plan and executing the plan, instead of doing x=fft(y). As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change. BTW memmap arrays have the same problem if I create a memmap array and later do something like a=a+1 all later changes will not be written to the file. BTW I just answered to you in the other thread on scipy-users as well, I'll try to just keep the rest of my posts here so people don't have to look at two lists if they want to follow things. Cheers Jochen
______ _________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Jochen wrote:
On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote:
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag:
http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff...
cheers,
David
Hi David,
I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan.
I am not sure I follow you when you say the "fftw way": you can avoid having to store the arrays altogether while still using fftw plans, there is nothing "unfftw" about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated.
David
Sorry maybe I wasn't quite clear, what I mean by the "fftw way" is creating a plan and executing the plan, instead of doing x=fft(y).
I guess I was not very clear, because my suggestion has nothing to do with the above.
As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change.
I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO. Another solution would be to work on numpy itself, so that it use aligned buffers, but that's obviously much longer term - that's just one of those things on my TODO list that I never took the time to do, David
On Tue, 2009-01-27 at 13:54 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote:
Jochen wrote:
Hi all,
I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions:
In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this:
a = fftw3.AlignedArray(1024,complex)
a = a+1
a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me
. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag:
http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff...
cheers,
David
Hi David,
I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan.
I am not sure I follow you when you say the "fftw way": you can avoid having to store the arrays altogether while still using fftw plans, there is nothing "unfftw" about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated.
David
Sorry maybe I wasn't quite clear, what I mean by the "fftw way" is creating a plan and executing the plan, instead of doing x=fft(y).
I guess I was not very clear, because my suggestion has nothing to do with the above.
As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change.
I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO.
Ah ok, I think I understand you now :). I agree that using fftw_execute is not the most pythonic way of doing things. However every other way I could think of, you would need to do a number of checks and possibly create new plans, as well as keeping old plans around when doing a fft. (I believe that's what you did in the old fftw bindings in scipy?). So I decided to create fftw bindings first for maximum performance (I know, just doing the whole calculation in c is probably more appropriate if you want maximum performance ;), it also fits my needs. However I've been thinking about ways to make the whole thing more intuitive. At least now I can compare how much overhead I'm adding if I start to do checks in order to make the handling easier :)
Another solution would be to work on numpy itself, so that it use aligned buffers, but that's obviously much longer term - that's just one of those things on my TODO list that I never took the time to do,
I agree that would be nice, I'd think a couple of operations would probably benefit from that.
David _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Jochen wrote:
On Tue, 2009-01-27 at 13:54 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote:
Jochen wrote:
> Hi all, > > I just wrote ctypes bindings to fftw3 (see > http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html > for the post to scipy). > Now I have a couple of numpy related questions: > > In order to be able to use simd instructions I > create an ndarray subclass, which uses fftw_malloc to allocate the > memory and fftw_free to free the memory when the array is deleted. This > works fine for inplace operations however if someone does something like > this: > > a = fftw3.AlignedArray(1024,complex) > > a = a+1 > > a.ctypes.data points to a different memory location (this is actually an > even bigger problem when executing fftw plans), however > type(a) still gives me
. > > > > I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff...
cheers,
David
Hi David,
I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan.
I am not sure I follow you when you say the "fftw way": you can avoid having to store the arrays altogether while still using fftw plans, there is nothing "unfftw" about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated.
David
Sorry maybe I wasn't quite clear, what I mean by the "fftw way" is creating a plan and executing the plan, instead of doing x=fft(y).
I guess I was not very clear, because my suggestion has nothing to do with the above.
As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change.
I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO.
Ah ok, I think I understand you now :). I agree that using fftw_execute is not the most pythonic way of doing things. However every other way I could think of, you would need to do a number of checks and possibly create new plans, as well as keeping old plans around when doing a fft.
Yes, but I don't see the problem with keeping old plans/creating new ones. It was a pain in C++, because I had to handle many backends, not just fftw, but in your case, doing it in python is not very difficult I think.
(I believe that's what you did in the old fftw bindings in scipy?).
that's what was done in scipy.fftpack originally, but not by me.
So I decided to create fftw bindings first for maximum performance (I know, just doing the whole calculation in c is probably more appropriate if you want maximum performance ;), it also fits my needs.
Caching plans will not affect performances; on the contrary, you will be able to use more aggressively optimized plans if you add the ability to save/load plans to the disk (using wisdow mechanism). Checking for aligned arrays is one line in python, and should not affect much performances, specially if you do fft on big arrays, David
On Tue, 2009-01-27 at 14:46 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 13:54 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote:
Jochen wrote:
On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote:
> Jochen wrote: > > > >> Hi all, >> >> I just wrote ctypes bindings to fftw3 (see >> http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html >> for the post to scipy). >> Now I have a couple of numpy related questions: >> >> In order to be able to use simd instructions I >> create an ndarray subclass, which uses fftw_malloc to allocate the >> memory and fftw_free to free the memory when the array is deleted. This >> works fine for inplace operations however if someone does something like >> this: >> >> a = fftw3.AlignedArray(1024,complex) >> >> a = a+1 >> >> a.ctypes.data points to a different memory location (this is actually an >> even bigger problem when executing fftw plans), however >> type(a) still gives me
. >> >> >> >> > I can't comment about subclassing ndarrays, but I can give you a hint > about aligned allocator problem: you could maintain two list of cached > plans, automatically detect whether your arrays are aligned or not, and > use the appropriate list of plans; one list is for aligned arrays, one > for unaligned. Before removing support for fftw, I played with some C++ > code to do exactly that. You can tell fftw to create plans for unaligned > arrays by using FFTW_UNALIGNED flag: > > http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/ff... > > cheers, > > David > > > Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan.
I am not sure I follow you when you say the "fftw way": you can avoid having to store the arrays altogether while still using fftw plans, there is nothing "unfftw" about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated.
David
Sorry maybe I wasn't quite clear, what I mean by the "fftw way" is creating a plan and executing the plan, instead of doing x=fft(y).
I guess I was not very clear, because my suggestion has nothing to do with the above.
As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change.
I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO.
Ah ok, I think I understand you now :). I agree that using fftw_execute is not the most pythonic way of doing things. However every other way I could think of, you would need to do a number of checks and possibly create new plans, as well as keeping old plans around when doing a fft.
Yes, but I don't see the problem with keeping old plans/creating new ones. It was a pain in C++, because I had to handle many backends, not just fftw, but in your case, doing it in python is not very difficult I think.
(I believe that's what you did in the old fftw bindings in scipy?).
that's what was done in scipy.fftpack originally, but not by me.
So I decided to create fftw bindings first for maximum performance (I know, just doing the whole calculation in c is probably more appropriate if you want maximum performance ;), it also fits my needs.
Caching plans will not affect performances; on the contrary, you will be able to use more aggressively optimized plans if you add the ability to save/load plans to the disk (using wisdow mechanism). Checking for aligned arrays is one line in python, and should not affect much performances, specially if you do fft on big arrays,
I agree checking for aligned arrays does not affect performance much, I guess the lookup mechanism for the plans is what would take up the most time IMO. I'll definitely look into this though :) Cheers Jochen
On 1/27/2009 6:03 AM, Jochen wrote:
BTW memmap arrays have the same problem if I create a memmap array and later do something like a=a+1 all later changes will not be written to the file.
= is Python's rebinding operator. a = a + 1 rebinds a to a different object. As for ndarray's, I'd like to point out the difference between a[:] = a + 1 and a = a + 1 S.M.
On 1/27/2009 1:26 AM, Jochen wrote:
a = fftw3.AlignedArray(1024,complex)
a = a+1
= used this way is not assignment, it is name binding. It is easy to use function's like fftw_malloc with NumPy: import ctypes import numpy fftw_malloc = ctypes.cdll.fftw.fftw_malloc fftw_malloc.argtypes = [ctypes.c_ulong,] fftw_malloc.restype = ctypes.c_ulong def aligned_array(N, dtype): d = dtype() address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong if (address = 0): raise MemoryError, 'fftw_malloc returned NULL' class Dummy(object): pass d = Dummy() d.__array_interface__ = { 'data' = (address, False), 'typestr' : dtype.str, 'descr' : dtype.descr, 'shape' : shape, 'strides' : strides, 'version' : 3 } return numpy.asarray(d) If you have to check for a particular alignment before calling fftw, that is trivial as well: def is_aligned(array, boundary): address = array.__array_interface__['data'][0] return not(address % boundary)
there a way that I get a different object type? Or even better is there a way to prevent operations like a=a+1 or make them automatically in-place operations?
a = a + 1 # rebinds the name 'a' to another array a[:] = a + 1 # fills a with the result of a + 1 This has to do with Python syntax, not NumPy per se. You cannot overload the behaviour of Python's name binding operator (=). Sturla Molden
On 1/27/2009 12:37 PM, Sturla Molden wrote:
It is easy to use function's like fftw_malloc with NumPy:
Besides this, if I were to write a wrapper for FFTW in Python, I would consider wrapping FFTW's Fortran interface with f2py. It is probably safer, as well as faster, than using ctypes. It would also allow the FFTW library to be linked statically to the Python extension, avoiding DLL hell. S.M.
On 1/27/2009 12:37 PM, Sturla Molden wrote:
import ctypes import numpy
fftw_malloc = ctypes.cdll.fftw.fftw_malloc fftw_malloc.argtypes = [ctypes.c_ulong,] fftw_malloc.restype = ctypes.c_ulong
def aligned_array(N, dtype): d = dtype() address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong if (address = 0): raise MemoryError, 'fftw_malloc returned NULL' class Dummy(object): pass d = Dummy() d.__array_interface__ = { 'data' = (address, False), 'typestr' : dtype.str, 'descr' : dtype.descr, 'shape' : shape, 'strides' : strides, 'version' : 3 } return numpy.asarray(d)
Or if you just want to use numpy, aligning to a 16 byte boundary can be done like this: def aligned_array(N, dtype): d = dtype() tmp = numpy.array(N * d.nbytes + 16, dtype=numpy.uint8) address = tmp.__array_interface__['data'][0] offset = (16 - address % 16) % 16 return tmp[offset:offset+N].view(dtype=dtype) S.M.
On Tue, 2009-01-27 at 14:16 +0100, Sturla Molden wrote:
On 1/27/2009 12:37 PM, Sturla Molden wrote:
import ctypes import numpy
fftw_malloc = ctypes.cdll.fftw.fftw_malloc fftw_malloc.argtypes = [ctypes.c_ulong,] fftw_malloc.restype = ctypes.c_ulong
def aligned_array(N, dtype): d = dtype() address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong if (address = 0): raise MemoryError, 'fftw_malloc returned NULL' class Dummy(object): pass d = Dummy() d.__array_interface__ = { 'data' = (address, False), 'typestr' : dtype.str, 'descr' : dtype.descr, 'shape' : shape, 'strides' : strides, 'version' : 3 } return numpy.asarray(d)
Or if you just want to use numpy, aligning to a 16 byte boundary can be done like this:
def aligned_array(N, dtype): d = dtype() tmp = numpy.array(N * d.nbytes + 16, dtype=numpy.uint8) address = tmp.__array_interface__['data'][0] offset = (16 - address % 16) % 16 return tmp[offset:offset+N].view(dtype=dtype)
S.M.
Ah, I didn't think about doing it in python, cool thanks.
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Tue, 2009-01-27 at 14:16 +0100, Sturla Molden wrote:
def aligned_array(N, dtype): d = dtype() tmp = numpy.array(N * d.nbytes + 16, dtype=numpy.uint8) address = tmp.__array_interface__['data'][0] offset = (16 - address % 16) % 16 return tmp[offset:offset+N].view(dtype=dtype)
Ah, I didn't think about doing it in python, cool thanks.
Doing it from Python means you don't have to worry about manually deallocating the array afterwards. It seems the issue of 16-byte alignment has to do with efficient data alignment for SIMD instructions (SSE, MMX, etc). So this is not just an FFTW issue. I would just put a check for 16-byte alignment in the wrapper, and raise an exception (e.g. MemoryError) if the array is not aligned properly. Raising an exception will inform the user of the problem. I would not attempt to make a local copy if the array is errorneously aligned. That is my personal preference. Sturla Molden
On Tue, 2009-01-27 at 12:37 +0100, Sturla Molden wrote:
On 1/27/2009 1:26 AM, Jochen wrote:
a = fftw3.AlignedArray(1024,complex)
a = a+1
= used this way is not assignment, it is name binding.
It is easy to use function's like fftw_malloc with NumPy:
import ctypes import numpy
fftw_malloc = ctypes.cdll.fftw.fftw_malloc fftw_malloc.argtypes = [ctypes.c_ulong,] fftw_malloc.restype = ctypes.c_ulong
def aligned_array(N, dtype): d = dtype() address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong if (address = 0): raise MemoryError, 'fftw_malloc returned NULL' class Dummy(object): pass d = Dummy() d.__array_interface__ = { 'data' = (address, False), 'typestr' : dtype.str, 'descr' : dtype.descr, 'shape' : shape, 'strides' : strides, 'version' : 3 } return numpy.asarray(d)
I actually do it slightly different, because I want to free the memory using fftw_free. So I'm subclassing ndarray and in the __new__ function of the subclass I create a buffer object from fftw_malloc using PyBuffer_FromReadWriteMemory and pass that to ndarray.__new__. In __del__ I check if the .base is a buffer object and do an fftw_free.
If you have to check for a particular alignment before calling fftw, that is trivial as well:
def is_aligned(array, boundary): address = array.__array_interface__['data'][0] return not(address % boundary)
Yes I knew that, btw do you know of any systems which need alignment to boundaries other than 16byte for simd operations?
there a way that I get a different object type? Or even better is there a way to prevent operations like a=a+1 or make them automatically in-place operations?
a = a + 1 # rebinds the name 'a' to another array
a[:] = a + 1 # fills a with the result of a + 1
I knew about this one, this is what I have been doing.
This has to do with Python syntax, not NumPy per se. You cannot overload the behaviour of Python's name binding operator (=).
Yes I actually thought about it some more yesterday when going home and realised that this would really not be possible. I guess I just have to document it clearly so that people don't run into it. Cheers Jochen
Sturla Molden
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
participants (4)
-
David Cournapeau -
Jochen -
Ryan May -
Sturla Molden