;;; Level3-fft-fftw.lisp, fast fourier by fftw. 

;;; 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.

;;; There is a sligh missuse of notation here, since the method 
;;; specialize on matrix-foreign-zge, rather than matrix-fftw-zge, which
;;; would have a more correct name.

(in-package :lisplab)

(defmethod fft1 ((x matrix-foreign-dge) &key)
  (fft1 (convert-to-matrix-zge x)))

(defmethod ifft1 ((x matrix-foreign-dge) &key)
  (ifft1 (convert-to-matrix-zge x)))

(defmethod fft2 ((x matrix-foreign-dge) &key)
  (fft2 (convert-to-matrix-zge x)))

(defmethod ifft2 ((x matrix-foreign-dge) &key)
  (ifft2 (convert-to-matrix-zge x)))

(defun use-fftw-p ()
  cl-user::*lisplab-libfftw-path*)

(defun fftw-use-threads-p ()
  (and cl-user::*lisplab-libfftw-threads-path*
       (> *lisplab-num-threads* 1)))

(defmethod init-threads :after (num-threads)
  (when (fftw-use-threads-p)
    (fftw-ffi:fftw-init-threads num-threads)))

(defmethod cleanup-threads :after ()
  (when (fftw-use-threads-p)
    (fftw-ffi:fftw-cleanup-threads)))

(defmethod fft1 ((x matrix-foreign-zge) &key)
   (if (not (use-fftw-p))
      (call-next-method)
      (let* ((rows (rows x))
	     (cols (cols x))
	     (store-x (vector-store x))
	     (y (mcreate x))
	     (store-y (vector-store y)))	 
	(dotimes (i cols)
	  (fftw-ffi:fftw-fft1 
	   rows
	   store-x
	   (* 2 i rows)
	   store-y
	   (* 2 i rows)
	   fftw-ffi:+FFTW-FORWARD+
	   fftw-ffi:+FFTW-ESTIMATE+))
	y)))

(defmethod ifft1 ((x matrix-foreign-zge) &key)
   (if (not (use-fftw-p))
      (call-next-method)
      (let* ((rows (rows x))
	     (cols (cols x))
	     (store-x (vector-store x))
	     (y (mcreate x))
	     (store-y (vector-store y)))	 
	(dotimes (i cols)
	  (fftw-ffi:fftw-fft1 
	   rows
	   store-x
	   (* 2 i rows)
	   store-y
	   (* 2 i rows)
	   fftw-ffi:+FFTW-BACKWARD+
	   fftw-ffi:+FFTW-ESTIMATE+))
	y)))

(defmethod fft2 ((x matrix-foreign-zge) &key)
  (if (not (use-fftw-p))
      (call-next-method)
      (let ((y (mcreate x)))
	(fftw-ffi:fftw-fft2 
	 (rows x) 
	 (cols x) 
	 (vector-store x) 
	 (vector-store y) 
	 fftw-ffi:+fftw-forward+ 
	 fftw-ffi:+FFTW-ESTIMATE+)
	y)))

(defmethod ifft2 ((x matrix-foreign-zge) &key)
  (if (not (use-fftw-p))
      (call-next-method)
      (let ((y (mcreate x)))
	(fftw-ffi:fftw-fft2 
	 (rows x) 
	 (cols x) 
	 (vector-store x) 
	 (vector-store y) 
	 fftw-ffi:+fftw-backward+ 
	 fftw-ffi:+FFTW-ESTIMATE+)
	y)))

;;; The below methods and functions with optimizations, such 
;;; as destructive transforms and real to complex transforms. Note
;;; that these are less tested than the other and more likly to have bugs.

;;; TODO: I should change the destructive methods so that they have explicite 
;;;       output also.

(defun fft1!-forward-or-backward (x direction)
  (let* ((rows (rows x))
	 (cols (cols x))
	 (store (vector-store x)))
    (dotimes (i cols)
      (fftw-ffi:fftw-fft1   
       rows
       store
       (* 2 i rows)
       store
       (* 2 i rows)
       direction
       fftw-ffi:+FFTW-ESTIMATE+)))
  x)

(defmethod fft1! ((x matrix-foreign-zge) &key)
  (if (not (use-fftw-p))
      (call-next-method)
      (fft1!-forward-or-backward x fftw-ffi:+fftw-forward+)))


(defmethod ifft1! ((x matrix-foreign-zge) &key)
  (if (not (use-fftw-p))
      (call-next-method)
      (fft1!-forward-or-backward x fftw-ffi:+fftw-backward+)))

(defmethod fft2! ((x matrix-foreign-zge) &key)
  (if (not (use-fftw-p))
      (call-next-method)
      (progn 
	(fftw-ffi:fftw-fft2 
	 (rows x) 
	 (cols x) 
	 (vector-store x) 
	 (vector-store x) 
	 fftw-ffi:+fftw-forward+ 
	 fftw-ffi:+FFTW-ESTIMATE+)
	x)))

(defmethod ifft2! ((x matrix-foreign-zge) &key)
  (if (not (use-fftw-p))
      (call-next-method)
      (progn
	(fftw-ffi:fftw-fft2 
	 (rows x) 
	 (cols x) 
	 (vector-store x) 
	 (vector-store x) 
	 fftw-ffi:+fftw-backward+ 
	 fftw-ffi:+FFTW-ESTIMATE+)
	x)))

(defmethod fftw1-r2c! ((x matrix-foreign-dge) (y matrix-foreign-zge))
  "Real to complex transform. y must have dimensions 1+(rows x)/2 X (cols x)."
  (let* ((rows-x (rows x))
	 (cols-x (cols x))
	 (rows-y (rows y))
	 (store-x (vector-store x))
	 (store-y (vector-store y)))
    (dotimes (i cols-x)
      (fftw-ffi:fftw-fft1-r2c   
       rows-x
       store-x
       (* i rows-x)
       store-y
       (* 2 i rows-y)
       fftw-ffi:+FFTW-ESTIMATE+)))
  y)

(defmethod fftw1-c2r! ((x matrix-foreign-zge) (y matrix-foreign-dge))
  "Complx to real inverse transform. x must have dimensions 1+(rows y)/2 X (cols y)."
  (let* ((rows-y (rows y))
	 (cols-y (cols y))
	 (rows-x (rows x))
	 (store-x (vector-store x))
	 (store-y (vector-store y)))
    (dotimes (i cols-y)
      (fftw-ffi:fftw-fft1-c2r   
       rows-y
       store-x
       (* 2 i rows-x)
       store-y
       (* i rows-y)
       fftw-ffi:+FFTW-ESTIMATE+)))
  y)

(defmethod fftw2-r2c! ((x matrix-foreign-dge) (y matrix-foreign-zge))
  "Real to complex transform. y must have dimensions 1+(rows x)/2 X (cols x)."
  (fftw-ffi:fftw-fft2-r2c  
   (rows x) 
   (cols x) 
   (vector-store x) 
   (vector-store y) 
   fftw-ffi:+FFTW-ESTIMATE+)
  y)

(defmethod fftw2-c2r! ((x matrix-foreign-zge) (y matrix-foreign-dge))
  "Complx to real inverse transform. x must have dimensions 1+(rows y)/2 X (cols y)."
  (fftw-ffi:fftw-fft2-c2r  
   (rows y) 
   (cols y) 
   (vector-store x) 
   (vector-store y) 
   fftw-ffi:+FFTW-ESTIMATE+)
  y)