fftw-fft.lisp

   1 (declaim (optimize (speed 2) (debug 3) (safety 3)))
   2 (in-package :vol)
   3
   4 (defconstant +forward+ 1)
   5 (defconstant +backward+ -1)
   6 (defconstant +estimate+ (ash 1 6))
   7
   8 (defmacro def-ffi (&optional (library "libfftw3.so") (prefix "fftw"))
   9   `(progn
  10      (load-shared-object ,(format nil "~a" library))
  11      (define-alien-routine (,(format nil "~a_execute" prefix)
  12                              ,(vol:format-symbol "~a-execute" prefix))
  13          void
  14        (plan (* int)))
  15      ,@(loop for rank in '(1 2 3) collect
  16             `(define-alien-routine
  17                  (,(format nil "~a_plan_dft_~ad" prefix rank)
  18                    ,(format-symbol "~a-plan-dft-~ad" prefix rank))
  19                  (* int)
  20                ,@(loop for i below rank collect
  21                       `(,(format-symbol "n~d" i) int))
  22                (in (* double-float))
  23                (out (* double-float))
  24                (sign int)
  25                (flags unsigned-int)))))
  26
  27 ;; define the foreign functions
  28 ;; fftw{f}-execute, fftw{f}-plan-dft-{1,2,3}d
  29 (def-ffi "libfftw3.so.3" "fftw")
  30 (def-ffi "libfftw3f.so.3" "fftwf")
  31
  32 #+nil
  33 (progn
  34   (load-shared-object "/usr/lib/libfftw3_threads.so")
  35
  36   (define-alien-routine ("fftw_init_threads" init-threads)
  37       int)
  38   (define-alien-routine ("fftw_plan_with_nthreads" plan-with-nthreads)
  39       void
  40     (nthreads int))
  41
  42   (defun init-ft ()
  43     (init-threads)
  44     (plan-with-nthreads 4)))
  45
  46 ;; to clean up completely call void fftw_cleanup_threads(void)
  47
  48
  49
  50 (def-generator (ft (rank type))
  51   (let* ((prefix (ecase type
  52                   (cdf 'fftw)
  53                   (csf 'fftwf)))
  54          (rlong-type (ecase type
  55                        (csf (get-long-type 'sf))
  56                        (cdf (get-long-type 'df))))
  57          (plan (format-symbol "~a-plan-dft-~ad" prefix rank))
  58          (execute (format-symbol "~a-execute" prefix)))
  59     `(defun ,name (in &key (forward t))
  60        (declare ((simple-array ,long-type ,rank) in)
  61                 (boolean forward)
  62                (values (simple-array ,long-type ,rank) &optional))
  63       (let* ((dims (array-dimensions in))
  64              (out (make-array dims :element-type ',long-type))
  65              (in-sap (sb-sys:vector-sap
  66                       (sb-ext:array-storage-vector in)))
  67              (out-sap (sb-sys:vector-sap
  68                        (sb-ext:array-storage-vector out)))
  69              (dir (if forward +forward+ +backward+)))
  70         (sb-sys:with-pinned-objects (in out)
  71           ,(ecase rank
  72                   (1 `(destructuring-bind (x) dims
  73                         (let ((p (,plan x in-sap out-sap
  74                                         dir +estimate+)))
  75                           (,execute p))))
  76                   (2 `(destructuring-bind (y x) dims
  77                         (let ((p (,plan y x in-sap out-sap
  78                                         dir +estimate+)))
  79                           (,execute p))))
  80                   (3 `(destructuring-bind (z y x) dims
  81                         (let ((p (,plan z y x in-sap out-sap
  82                                         dir +estimate+)))
  83                           (,execute p))))))
  84         (when forward ;; normalize if forward
  85           (s* (/ ,(coerce 1 rlong-type) (array-total-size out)) out))
  86         out))))
  87
  88 #+nil
  89 (def-ft-rank-type 1 csf)
  90
  91 (defmacro def-ft-functions (ranks types)
  92   (let* ((specifics nil)
  93          (cases nil)
  94          (name (format-symbol "ft")))
  95     (loop for rank in ranks do
  96          (loop for type in types do
  97               (let ((def-name (format-symbol "def-~a-rank-type" name))
  98                     (specific-name (format-symbol "~a-~a-~a" name rank type)))
  99                 (push `(,def-name ,rank ,type) specifics)
 100                 (push `((simple-array ,(get-long-type type) ,rank)
 101                         (,specific-name a :forward forward))
 102                       cases))))
 103     (store-new-function name)
 104     `(progn ,@specifics
 105             (defun ,name (a &key forward)
 106                (etypecase a
 107                  ,@cases
 108                  (t (error "The given type can't be handled with a generic ~a function." ',name)))))))
 109
 110 (def-ft-functions (1 2 3) (csf cdf))
 111
 112 (defmacro ift (in)
 113   `(ft ,in :forward nil))
 114
 115 #+nil
 116 (progn
 117   (time
 118    (let* ((nx 128)
 119           (ny nx)
 120           (nz ny)
 121           (a (convert-3-ub8/csf-mul
 122               (draw-sphere-ub8 20.0 nz ny nx))))
 123      (write-pgm "/home/martin/tmp/fftw.pgm"
 124                 (normalize-2-csf/ub8-abs
 125                  (cross-section-xz-csf (ft a)))))))