;;; Foreign function interfaces for FFTW version 3.

;;; Copyright (C) 2009 Joern Inge Vestgaarden
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License along
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.


(in-package :fftw-ffi)

(defconstant +double-float-bytes+ (truncate (sb-alien:ALIEN-SIZE sb-alien:double-float) 8))

(defconstant +FFTW-ESTIMATE+ 64)

(defconstant +FFTW-FORWARD+ -1)

(defconstant +FFTW-BACKWARD+ 1)

;;; The interface functions

(defun fftw-fft1 (n a astart b bstart direction flag)
  "One dimensional fft by forign call to fftw."
  (let ((astart (* astart +double-float-bytes+))
	(bstart (* bstart +double-float-bytes+)))
    (with-pinned-objects (a b)  
      (let ((plan (|fftw_plan_dft_1d| 
		   n  
		   (sap+ (vector-sap a) astart) 
		   (sap+ (vector-sap b) bstart)
		   direction
		   flag)))
	(unwind-protect 
	     (|fftw_execute| plan)
	  (|fftw_destroy_plan| plan))))
    b))

(defun fftw-fft2 (m n in out direction flag)
  "Two dimensional fft by forign call to fftw."
  (with-pinned-objects (in out)  
    (let ((plan (|fftw_plan_dft_2d| 
		 n ; swap n and m due to row major order
		 m
		 (vector-sap in)
		 (vector-sap out)		 
		 direction
		 flag)))
      (unwind-protect 
	   (|fftw_execute| plan)
	(|fftw_destroy_plan| plan))))
  out)

(defun fftw-fft1-r2c (n a astart b bstart flag)
  "One dimensional real fft by forign call to fftw. 
The length of b must at least be n/2+1."
  (let ((astart (* astart +double-float-bytes+))
	(bstart (* bstart +double-float-bytes+)))
    (with-pinned-objects (a b)  
      (let ((plan (|fftw_plan_dft_r2c_1d| 
		   n  
		   (sap+ (vector-sap a) astart) 
		   (sap+ (vector-sap b) bstart)
		   flag)))
	(unwind-protect 
	     (|fftw_execute| plan)
	  (|fftw_destroy_plan| plan))))
    b))

(defun fftw-fft1-c2r (n a astart b bstart flag)
  "One dimensional reverse fft transform by forign call to fftw."
  (let ((astart (* astart +double-float-bytes+))
	(bstart (* bstart +double-float-bytes+)))
    (with-pinned-objects (a b)
      (let ((plan (|fftw_plan_dft_c2r_1d| 
		   n  
		   (sap+ (vector-sap a) astart) 
		   (sap+ (vector-sap b) bstart)
		   flag)))
	(unwind-protect 
	     (|fftw_execute| plan)
	  (|fftw_destroy_plan| plan))))
    b))

(defun fftw-fft2-r2c (m n in out flag)
  "Two dimensional fft by forign call to fftw. Real to complex. Length of complex must
at least be n/2+1."
  (with-pinned-objects (in out)  
    (let ((plan (|fftw_plan_dft_r2c_2d| 
		 n ; swap n and m due to row major order
		 m
		 (vector-sap in)
		 (vector-sap out)		 
		 flag)))
      (unwind-protect 
	   (|fftw_execute| plan)
	(|fftw_destroy_plan| plan))))
  out)

(defun fftw-fft2-c2r (m n in out flag)
  "Two dimensional inverse fft by forign call to fftw. Real to complex. 
Length of complex must at least be n/2+1."
  (with-pinned-objects (in out)  
    (let ((plan (|fftw_plan_dft_c2r_2d| 
		 n ; swap n and m due to row major order
		 m
		 (vector-sap in)
		 (vector-sap out)		 
		 flag)))
      (unwind-protect 
	   (|fftw_execute| plan)
	(|fftw_destroy_plan| plan))))
  out)




;;; The FFI definitions

(declaim (inline |fftw_plan_dft_1d|))
(define-alien-routine |fftw_plan_dft_1d|
    (* t)
  (n int)
  (in (* double-float))
  (out (* double-float))
  (sign int)
  (flags int))
 
(declaim (inline |fftw_plan_dft_2d|))
(define-alien-routine |fftw_plan_dft_2d|
    (* t)
  (n0 int)
  (n1 int)
  (in (* double-float))
  (out (* double-float))
  (sign int)
  (flags int))

(declaim (inline |fftw_execute|))
(define-alien-routine |fftw_execute| 
    void
  (plan (* t)))

(declaim (inline |fftw_destroy_plan|))
(define-alien-routine |fftw_destroy_plan| 
    void 
  (plan (* t)))

;;;; Now multi-thread code

(declaim (inline |fftw_init_threads|))
(define-alien-routine |fftw_init_threads|
    int)
  
(declaim (inline |fftw_plan_with_nthreads|))
(define-alien-routine |fftw_plan_with_nthreads|
    void
  (nthreads int))

(declaim (inline |fftw_cleanup_threads|))
(define-alien-routine |fftw_cleanup_threads|
    void)

(defun fftw-init-threads (num-threads)
  ;; Note: assumes that |fftw_init_threads| has been called!
  (|fftw_plan_with_nthreads| num-threads))

(defun fftw-cleanup-threads ()
  (|fftw_cleanup_threads|))

;;; Now real to complex transforms

(declaim (inline |fftw_plan_dft_r2c_1d|))
(define-alien-routine |fftw_plan_dft_r2c_1d|
    (* t)
  (n int)
  (in (* double-float))
  (out (* double-float))
  (flags int))

(declaim (inline |fftw_plan_dft_c2r_1d|))
(define-alien-routine |fftw_plan_dft_c2r_1d|
    (* t)
  (n int)
  (in (* double-float))
  (out (* double-float))
  (flags int))

(declaim (inline |fftw_plan_dft_r2c_2d|))
(define-alien-routine |fftw_plan_dft_r2c_2d|
    (* t)
  (n0 int)
  (n1 int)
  (in (* double-float))
  (out (* double-float))
  (flags int))

(declaim (inline |fftw_plan_dft_c2r_2d|))
(define-alien-routine |fftw_plan_dft_c2r_2d|
    (* t)
  (n0 int)
  (n1 int)
  (in (* double-float))
  (out (* double-float))
  (flags int))