On Fri, Nov 18, 2011 at 5:18 AM, Dag Sverre Seljebotn < [email protected]> wrote:
> On 11/18/2011 12:58 PM, David Cournapeau wrote: > > On Fri, Nov 18, 2011 at 11:42 AM, Robert Kern<[email protected]> > wrote: > >> On Fri, Nov 18, 2011 at 11:19, Dag Sverre Seljebotn > >> <[email protected]> wrote: > >>> I've been in touch with Martin Reinecke, author of the libpsht code for > >>> spherical harmonic transforms, about licensing issues. > >>> > >>> libpsht itself will remain under the GPL, but he is likely to release > >>> his C port of FFTPACK under BSD in the near future, as it is based on > >>> the public domain FFTPACK. > >>> > >>> I'm grateful for this change for my own purposes (allows releasing my > >>> own competing SHT library under the BSD) -- but it could perhaps be > >>> useful for NumPy or SciPy as well, depending on how complete the port > >>> is? E.g., perhaps make numpy.fft more complete (is the > >>> numpy.fft/scipy.fftpack split simply because of the Fortran > dependency?). > >> > >> It used to be the case that scipy.fftpack allowed one to build against > >> multiple different, usually faster, FFT libraries like FFTW. I think > >> we have backed away from that since the cost of maintaining the build > >> configuration for all of those different backends was so high. It's > >> worth noting that numpy.fft is already using a C translation of > >> FFTPACK. I'm not sure what the differences are between this > >> translation and Martin's. > > Here's some more info forwarded from Martin: > > """ > - only FFTs are supported (no DCTs/DSTs) > - only double precision is supported (extension to single precision might > not be much work, though) > - both complex and real FFTs are supported > - real FFTs allow various storage schemes for the (half)complex frequency > domain data (classic FFTPACK scheme, FFTW or halfcomplex scheme, > uncompressed complex storage) > - precision of transforms involving large prime factors should be slightly > better than with original FFTPACK > - Bluestein's algorithm is automatically selected if considered profitable > - small accuracy self-testing code is provided. > > Fairly complete interface documentation is available in Doxygen format. > I'll prepare a source package later in the afternoon and send it around. > > Best regards, > Martin > """ > > > > > Having a Bluestein transformation alone would be worthwhile, as it > > would avoid the N^2 penalty for prime sizes. > > > > I am wondering about precision issues, though (when I tried > > implementing bluestein transforms on top of fftpack, it gave very bad > > results numerically-wise). A comparison with fftw would be good here. > > Well, there's an indirect comparison: My SHT code currently uses FFTW3, > and it manages to agree with Reinecke's SHT code to a precision of > better than ~1e-12 for full SHTs. That includes several other sources of > errors though. > > (That's an average over several different-sized FFTs, of which half has > n=8192 and the other half all have different size, starting from 4 and > increasing up to 8192 in steps of 4 -- meaning prime factors on the > order of 1000). > > I agree, a more direct comparison with FFTW would be good. > > In more detail from the README: > > I replaced the iterative sine and cosine calculations in radfg() and > radbg() > 17 by an exact calculation, which slightly improves the transform > accuracy for > 18 real FFTs with lengths containing large prime factors. > > The current numpy FFT code could use a review in any case, IIRC, there were some possible integer related issues in the arguments. Having clean documented C code would probably be an improvement and also help with maintenance. I suspect there were iterative computations of sin/cos that were spoiling the precision of the single precision versions (scipy), so fixing those might also be useful. Chuck
_______________________________________________ NumPy-Discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
