;
; Copyright 2011,2012 Lajos Foldy
;
; IDL interface to FFTW3 (call_external)
;
; Result = FL_FFTW3( Array [, Direction] [, DIMENSION=], [, /DOUBLE] [, /INVERSE] [, /OVERWRITE] )
;
; License: GPLv2 (http://www.gnu.org/licenses/gpl-2.0.html)
;
; 20111031  Original version, tested in Linux only
; 20120808  Multidimensional FFT fixed
; 20121005  DIMENSION added
;
; FFTW3: download, configure (with --with-combined-threads), compile and install.
; Change the library locations at line 51.
;

function fl_fftw3, data, direction, dimension=dimension, double=double, inverse=inverse, overwrite=overwrite

dir=-1l
if n_elements(direction) ne 0 then $
   if direction gt 0 then dir=1l
if keyword_set(inverse) then dir=1l

dble=-1
if n_elements(double) ne 0 then dble=keyword_set(double)

size=size(data)
ndim=size[0]
type=size[ndim+1]

dimen=0
if n_elements(dimension) ne 0 then begin
   if n_elements(dimension) gt 1 then message, 'scalar required: DIMENSION'
   dimen=long(dimension[0])
   if dimen lt 0 || dimen gt ndim then message, 'invalid value: DIMENSION'
   if ndim eq 1 then dimen=0
endif

over=keyword_set(overwrite)

if type eq 0 or type eq 7 or type eq 8 or type eq 10 or type eq 11 then message, 'numeric data required'
if ndim eq 0 then message, 'array required'
dims=long(size[1:ndim])

if dble eq -1 then begin
   rtype=6
   if type eq 5 or type eq 9 then rtype=9
endif else rtype=6+3*dble

library= rtype eq 6 ? '/tmp/fftw3/lib/libfftw3f.so' : '/tmp/fftw3/lib/libfftw3.so'
prefix= rtype eq 6 ? 'fftwf_' : 'fftw_'

x=call_external(library, prefix+'init_threads', /auto_glue)
x=call_external(library, prefix+'plan_with_nthreads', long(!cpu.tpool_nthreads), value=[1], /auto_glue)

if dimen eq 0 then begin
   if over then begin
      if rtype ne type then data=fix(data, type=rtype)
      plan=call_external(library, prefix+'plan_dft', ndim, reverse(dims), $
                         data, data, dir, 64l, value=[1,0,0,0,1,1], /ul64_value, /auto_glue)
   endif else begin
      out=rtype eq type ? data : fix(data, type=rtype)
      plan=call_external(library, prefix+'plan_dft', ndim, reverse(dims), out, out, $
                         dir, 64l, value=[1,0,0,0,1,1], /ul64_value, /auto_glue)
   endelse
   x=call_external(library, prefix+'execute', plan, value=[1], /auto_glue)
   x=call_external(library, prefix+'destroy_plan', plan, value=[1], /auto_glue)
   scale= dir eq -1l ? 1.0/size[ndim+2] : 0.0
endif else begin
   n1=dimen gt 1 ? long(product(dims[0:dimen-2], /int)) : 1l
   n2=long(dims[dimen-1])
   n3=dimen lt ndim ? long(product(dims[dimen:*], /int)) : 1l
   n=[n2]
   out=rtype eq type ? data : fix(data, type=rtype)
   if n1 eq 1 then begin
      plan=call_external(library, prefix+'plan_many_dft', 1l, n, n3, out, n, 1l, n2, out, n, 1l, n2, $
                         dir, 64l, value=[1,0,1,0,0,1,1,0,0,1,1,1,1], /ul64_value, /auto_glue)
      x=call_external(library, prefix+'execute', plan, value=[1], /auto_glue)
      x=call_external(library, prefix+'destroy_plan', plan, value=[1], /auto_glue)
   endif else begin
      out=reform(out, n1,n2,n3, /over)
      for j=0l,n3-1 do begin
          tmp=out[*,*,j]
          plan=call_external(library, prefix+'plan_many_dft', 1l, n, n1, tmp, n, n1, 1l, tmp, n, n1, 1l, $
                             dir, 64l, value=[1,0,1,0,0,1,1,0,0,1,1,1,1], /ul64_value, /auto_glue)
          x=call_external(library, prefix+'execute', plan, value=[1], /auto_glue)
          x=call_external(library, prefix+'destroy_plan', plan, value=[1], /auto_glue)
          out[*,*,j]=tmp
      endfor
      out=reform(out, dims, /over)
   endelse
   scale= dir eq -1l ? 1.0/n2 : 0.0
   if over then data=out
endelse

if scale ne 0.0 then $
   if over then data=data*scale else out=out*scale

return, over ? data : out
end