Yes, I believe that results from fftw3 agree with fftw2. The tests below pass with the following exceptions:
- those involving A34. I fancy we're getting different random values. Certainly I get different random values on w32 than I do on w64 (also they're different machines; I don't know which difference is important). Therefore FFTA34 doesn't match. - The sentence f?.(6?.10)$100 causes J to exit unexpectedly on w32. (But this test does pass on w64). - I get the crash with FFTW2 as well, so I have at least not made this worse. - I don't know whether this is a J bug or a bug in FFTW or a bug in the script joining them together - The interface has changed enough in FFTW3 that the final set of tests are not correct - the in/out parameters are now passed to createplan, and fftwndone is replaced by fftwexecute David "bill lam" <[email protected]> wrote in message news:[email protected]... > Does results from fftw3 agree with that from fftw2? Can you also > verify fftw pass these tests? > > > NB. test fftw system > NB. > NB. load this file with 0!:2 > NB. > NB. 0!:2 <'\jx\addon\fftw\fftw\testfftw.ijs' > > load '~system\packages\fftw\fftw.ijs' > > clean=: 1e_10&$: : (4 : 0) > if. (3!:0 y) e. 16 16384 do. > j./"1 y* x <: | y=. +.y > else. > y * x <: |y > end. > ) > > dfft=: 3 : '+/ y * ^ (#y) %~ (- o. 0j2 ) * */~ i.#y' > matchclean=: 0: *./ . = clean @ , @: - > round=: [ * [: <. 0.5"_ + %~ > > A34=: j./?.1+i.2 3 4 > > FFTA34=: ". ;._2 (0 : 0) > 39j99 _10j_10 10.3661j18.8301 9.3205j_6.6795 > _8j6 _33.0788j_17.2417 _13.9282j3.0526 8.634j10.1698 > _1j13 _25.3205j_41.3205 24.0789j_39.7584 _0.0718j_35.0526 > ) > > NB. ========================================================= > NB. test dat = ifftw fftw dat > > f=: matchclean if...@fftw > f i.8 > f i.4096 > f i.2 3 4 > f ?.(6?.10)$100 > > NB. ========================================================= > NB. test against direct calculation > > f=: dfft matchclean fftw > f i.8 > f i.132 > > NB. ========================================================= > NB. known examples: > a=: 2 0 1 0 0 0 1 0 > b=: 4 2 0 2 4 2 0 2 > > a matchclean ifftw b > b matchclean fftw a > FFTA34 matchclean 1e_4 round fftw A34 > > NB. ========================================================= > NB. createplan > NB. existimate: > P=: createplan_jfftw_ 3 4;_1;0 > FFTA34 matchclean 1e_4 round fftwndone_jfftw_ P;_1;A34 > FFTA34 matchclean 1e_4 round -: fftwndone_jfftw_ P;_1;A34*2 > destroyplan_jfftw_ P > > NB. performance measure: > P=: createplan_jfftw_ 3 4;_1;1 > FFTA34 matchclean 1e_4 round fftwndone_jfftw_ P;_1;A34 > FFTA34 matchclean 1e_4 round -: fftwndone_jfftw_ P;_1;A34*2 > destroyplan_jfftw_ P > > > dim, 21 Mar 2010, David Hotham skribis: >> "Therefore I suspect the problem originated from the window binary" >> >> It would be nice to put the blame elsewhere but, alas, this is not so. >> >> The issue seems to be that the definition of FFTW_ESTIMATE has changed >> between FFTW2 and FFTW3. What was zero is now (1U << 6). >> >> Here's a revised version of the script that: >> - incorporates your changes >> - removes some crud (or, at any rate, it seems to me to be crud; and >> removing it makes no difference on Windows) >> - fixes FFTW_ESTIMATE >> >> ... and (on w32 and w64) gives sensible results first time. >> >> David >> >> >> NB. built from project: ~Addons/math/fftw/fftw >> NB. z definitions: >> >> script_z_ '~system\main\dll.ijs' >> >> coclass 'jfftw' >> >> >> fftw_z_=: (_1 & fftwnd_jfftw_) :. (1 & fftwnd_jfftw_) >> ifftw_z_=: (1 & fftwnd_jfftw_) :. (_1 & fftwnd_jfftw_) >> >> NB. fftw utils >> NB. >> NB. cd 15!:0 >> NB. clean clean numbers near 0 >> NB. info cover for wdinfo >> NB. matchclean if clean x-y is all 0 >> >> zzero=: 1j1-1j1 >> >> 3 : 0'' >> if. IFUNIX do. >> DLL=: 'libfftw3.so.3' >> else. >> DLL=: '"',~'"',jpath '~addons\math\fftw\libfftw3-3.dll' >> end. >> ) >> >> FFTW_FORWARD=: _1 >> FFTW_BACKWARD=: 1 >> FFTW_ESTIMATE=: 6 (33 b.) 1 >> FFTW_MEASURE=: 0 >> >> FFTW_VERSION=: 3.2 >> >> cd=: 15!:0 >> >> info=: wdinfo @ ('FFTW'&;) >> >> matchclean=: 0: *./ . = clean @ , @: - >> >> NB. ========================================================= >> NB.*clean v clean numbers to tolerance (default 1e_10) >> NB. sets values less than tolerance to 0 >> clean=: (1e_10&$:) : (4 : 0) >> if. L. y do. >> x clean each y >> else. >> if. 16 ~: 3!:0 y do. >> y * x <: |y >> else. >> j./"1 y* x <: | y=. +.y >> end. >> end. >> ) >> >> >> NB. fftw >> >> NB. ========================================================= >> NB.*createplan v create a plan >> NB. y = shape ; in ; out ; direction; flag >> NB. >> NB. direction = FFTW_FORWARD | FFTW_BACKWARD >> NB. flag = FFTW_ESTIMATE | FFTW_MEASURE >> createplan=: 3 : 0 >> 'shape in out dir flag'=. y >> assert dir e. FFTW_FORWARD,FFTW_BACKWARD >> assert flag e. FFTW_ESTIMATE, FFTW_MEASURE >> shape=. ,shape >> dir=. dir >> flag=. flag >> cmd=. DLL,' fftw_plan_dft + x i *i *j *j i i' >> 0 pick cmd cd (#shape);shape;in;out;dir;flag >> ) >> >> NB. ========================================================= >> NB.*destroyplan v destroy a plan >> destroyplan=: 3 : 0 >> cmd=. DLL,' fftw_destroy_plan + n x' >> 1 [ cmd cd y >> ) >> >> NB. ========================================================= >> NB.*fftwnd d n-dimensional FFT >> NB. x = _1 forward >> NB. 1 backward >> NB. y = data >> fftwnd=: 4 : 0 >> shp=. $y >> in=. zzero + , |: y >> out=. in * 0 >> assert x e. _1 1 >> dir=. x >> plan=. createplan shp;in;out;dir;FFTW_ESTIMATE >> fftwexecute plan >> destroyplan plan >> res=. |: (|.shp) $ out >> if. dir=1 do. res % */shp end. >> ) >> >> NB. ========================================================= >> NB.*fftwexecute d one call to n-dimensional FFT >> NB. y = plan >> fftwexecute=: 3 : 0 >> cmd=. DLL,' fftw_execute + n x' >> 1 [ cmd cd y >> ) >> >> >> "bill lam" <[email protected]> wrote in message >> news:[email protected]... >> [---=| TOFU protection by t-prot: 136 lines snipped |=---] > > -- > regards, > ==================================================== > GPG key 1024D/4434BAB3 2008-08-24 > gpg --keyserver subkeys.pgp.net --recv-keys 4434BAB3 > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
