[lisplab-cvs] r176 - in trunk/src: draft fft
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Mon May 24 14:02:41 UTC 2010
Author: jivestgarden
Date: Mon May 24 10:02:40 2010
New Revision: 176
Log:
bugfix in fft + new destructive optimized methods for real to complex transforms
Modified:
trunk/src/draft/level0-expression.lisp
trunk/src/fft/fftw-ffi-package.lisp
trunk/src/fft/fftw-ffi.lisp
trunk/src/fft/level3-fft-fftw.lisp
Modified: trunk/src/draft/level0-expression.lisp
==============================================================================
--- trunk/src/draft/level0-expression.lisp (original)
+++ trunk/src/draft/level0-expression.lisp Mon May 24 10:02:40 2010
@@ -19,6 +19,229 @@
(in-package :lisplab)
+(defclass expression-base ()
+ ()
+ (:documentation "Represents symbolic expressions."))
+
+(defclass whatever (expression-base)
+ ()
+ (:documentation "Matches anything."))
+
+(defclass symbol-expression (expression-base)
+ ((symbol :initarg :symbol :accessor expression-symbol :initform nil))
+ (:documentation "Matches the same symbol. "))
+
+
+
+(defclass list-expression (expression-base)
+ ((list :accessor expression-list :initform nil :initarg :list))
+ (:documentation "A collection of expressions."))
+
+(defclass rule-base ()
+ ((name :accessor rule-name :initarg :name :initform nil)))
+
+(defclass function-rule (rule-base)
+ ((arg :accessor rule-arg :initarg :arg :initform nil)
+ (val :accessor rule-val :initarg :val :initform nil)))
+
+
+(defclass relation-rule (rule-base)
+ ((a :accessor rule-a :initarg :a :initform nil)
+ (b :accessor rule-b :initarg :b :initform nil)
+ (val :accessor rule-val :initarg :val :initform nil)))
+
+(defvar *rules* nil)
+
+(defun add-rule (rule)
+ (push rule *rules*))
+
+(defgeneric applicable-p (rule expr))
+
+(defgeneric apply-rule (rule expr))
+
+(defgeneric match-p (expr pat))
+
+(defmethod match-p (expr pat)
+ nil)
+
+(defmethod match-p (expr (pat whatever))
+ t)
+
+(defmethod match-p ((expr symbol-expression) (pat symbol-expression))
+ (eql (expression-symbol expr)
+ (expression-symbol pat)))
+
+(defmethod match-p ((expr list-expression) (pat list-expression))
+ (every #'match-p
+ (expression-list expr)
+ (expression-list pat)))
+
+
+
+(defmethod applicable-p (rule expr) nil)
+
+(defmethod applicable-p ((rule function-rule) (expr list-expression))
+ (let ((elms (expression-list expr)))
+ (if (< (length elms) 2)
+ nil
+ (and (match-p (rule-name rule) (car elms))
+ (match-p (rule-arg rule) (cadr elms))))))
+
+
+;;; Just simple integration with lisplab
+
+(defmethod print-object ((ex symbol-expression) stream)
+ (prin1 (expression-symbol ex) stream))
+
+(defmethod print-object ((ex list-expression) stream)
+ (prin1 (expression-list ex) stream))
+
+
+;;; Bellow is just cut and paste code to simplify debugging. Must do it properly later
+
+(defmethod .= ((a symbol) (b symbol) &optional acc)
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.=)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .= ((a symbol) (b number) &optional acc)
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.=)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .= ((a number) (b symbol) &optional acc)
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.=)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .= ((a expression-base) (b symbol) &optional acc)
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.=)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .= ((a symbol) (b expression-base) &optional acc)
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.=)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .= ((a expression-base) (b number) &optional acc)
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.=)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .= ((a number) (b expression-base) &optional acc)
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.=)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .add ((a symbol) (b symbol))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.+)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .add ((a symbol) (b number))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.+)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .add ((a number) (b symbol))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.+)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .add ((a expression-base) (b symbol))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.+)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .add ((a symbol) (b expression-base))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.+)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .add ((a expression-base) (b number))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.+)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .add ((a number) (b expression-base))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.+)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+
+
+
+(defmethod .mul ((a symbol) (b symbol))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.*)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .mul ((a symbol) (b number))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.*)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .mul ((a number) (b symbol))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.*)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .mul ((a expression-base) (b symbol))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.*)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .mul ((a symbol) (b expression-base))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.*)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .mul ((a expression-base) (b number))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.*)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+(defmethod .mul ((a number) (b expression-base))
+ (make-instance 'list-expression
+ :list (list (make-instance 'symbol-expression :symbol '.*)
+ (make-instance 'symbol-expression :symbol a)
+ (make-instance 'symbol-expression :symbol b))))
+
+
+
+
+
+
+
+
+
+
+
+
+#|
+
+
(defclass expression ()
((list :accessor expression-list :initform nil :initarg :list)))
@@ -274,4 +497,6 @@
(case b
(0 0)
(1 a)
- (t (call-next-method))))
\ No newline at end of file
+ (t (call-next-method))))
+
+|#
\ No newline at end of file
Modified: trunk/src/fft/fftw-ffi-package.lisp
==============================================================================
--- trunk/src/fft/fftw-ffi-package.lisp (original)
+++ trunk/src/fft/fftw-ffi-package.lisp Mon May 24 10:02:40 2010
@@ -22,6 +22,10 @@
"+FFTW-BACKWARD+"
"FFTW-FFT1"
"FFTW-FFT2"
+ "FFTW-FFT1-R2C"
+ "FFTW-FFT1-C2R"
+ "FFTW-FFT2-R2C"
+ "FFTW-FFT2-C2R"
"FFTW-INIT-THREADS"
"FFTW-CLEANUP-THREADS")
(:documentation "Simple ffi for fftw."))
Modified: trunk/src/fft/fftw-ffi.lisp
==============================================================================
--- trunk/src/fft/fftw-ffi.lisp (original)
+++ trunk/src/fft/fftw-ffi.lisp Mon May 24 10:02:40 2010
@@ -16,8 +16,6 @@
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-;;; TODO: the calls should be wrapped in unwind protect
-;;; to avoid memory leaks
(in-package :fftw-ffi)
@@ -29,38 +27,10 @@
(defconstant +FFTW-BACKWARD+ 1)
-(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)))
+;;; The interface functions
(defun fftw-fft1 (n a astart b bstart direction flag)
"One dimensional fft by forign call to fftw."
- ;; TODO we should handle conditions to avoid mem-leaks
(let ((astart (* astart +double-float-bytes+))
(bstart (* bstart +double-float-bytes+)))
(with-pinned-objects (a b)
@@ -77,7 +47,6 @@
(defun fftw-fft2 (m n in out direction flag)
"Two dimensional fft by forign call to fftw."
- ;; TODO we should handle conditions to avoid mem-leaks
(with-pinned-objects (in out)
(let ((plan (|fftw_plan_dft_2d|
n ; swap n and m due to row major order
@@ -91,6 +60,100 @@
(|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
@@ -114,3 +177,38 @@
(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))
\ No newline at end of file
Modified: trunk/src/fft/level3-fft-fftw.lisp
==============================================================================
--- trunk/src/fft/level3-fft-fftw.lisp (original)
+++ trunk/src/fft/level3-fft-fftw.lisp Mon May 24 10:02:40 2010
@@ -58,13 +58,12 @@
(y (mcreate x))
(store-y (matrix-store y)))
(dotimes (i cols)
- ;; Could be made in parallel
(fftw-ffi:fftw-fft1
rows
store-x
- (* i cols)
+ (* 2 i rows)
store-y
- (* i cols)
+ (* 2 i rows)
fftw-ffi:+FFTW-FORWARD+
fftw-ffi:+FFTW-ESTIMATE+))
y)))
@@ -78,13 +77,12 @@
(y (mcreate x))
(store-y (matrix-store y)))
(dotimes (i cols)
- ;; Could be made in parallel
(fftw-ffi:fftw-fft1
rows
store-x
- (* i cols)
+ (* 2 i rows)
store-y
- (* i cols)
+ (* 2 i rows)
fftw-ffi:+FFTW-BACKWARD+
fftw-ffi:+FFTW-ESTIMATE+))
y)))
@@ -115,7 +113,12 @@
fftw-ffi:+FFTW-ESTIMATE+)
y)))
-;;; TODO: remove the destructive mothods below. They only mess things up
+;;; 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))
@@ -125,25 +128,27 @@
(fftw-ffi:fftw-fft1
rows
store
- (* i cols)
+ (* 2 i rows)
store
- (* i cols)
+ (* 2 i rows)
direction
fftw-ffi:+FFTW-ESTIMATE+)))
x)
(defmethod fft1! ((x matrix-blas-zge) &key)
- (if cl-user::*lisplab-libfftw-path*
- (fft1!-forward-or-backward x fftw-ffi:+fftw-forward+)
- (call-next-method)))
+ (if (not (use-fftw-p))
+ (call-next-method)
+ (fft1!-forward-or-backward x fftw-ffi:+fftw-forward+)))
+
(defmethod ifft1! ((x matrix-blas-zge) &key)
- (if cl-user::*lisplab-libfftw-path*
- (fft1!-forward-or-backward x fftw-ffi:+fftw-backward+)
- (call-next-method)))
+ (if (not (use-fftw-p))
+ (call-next-method)
+ (fft1!-forward-or-backward x fftw-ffi:+fftw-backward+)))
(defmethod fft2! ((x matrix-blas-zge) &key)
- (if cl-user::*lisplab-libfftw-path*
+ (if (not (use-fftw-p))
+ (call-next-method)
(progn
(fftw-ffi:fftw-fft2
(rows x)
@@ -152,11 +157,11 @@
(matrix-store x)
fftw-ffi:+fftw-forward+
fftw-ffi:+FFTW-ESTIMATE+)
- x)
- (call-next-method)))
+ x)))
(defmethod ifft2! ((x matrix-blas-zge) &key)
- (if cl-user::*lisplab-libfftw-path*
+ (if (not (use-fftw-p))
+ (call-next-method)
(progn
(fftw-ffi:fftw-fft2
(rows x)
@@ -165,5 +170,58 @@
(matrix-store x)
fftw-ffi:+fftw-backward+
fftw-ffi:+FFTW-ESTIMATE+)
- x)
- (call-next-method)))
+ x)))
+
+(defmethod fftw1-r2c! ((x matrix-blas-dge) (y matrix-blas-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 (matrix-store x))
+ (store-y (matrix-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-blas-zge) (y matrix-blas-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 (matrix-store x))
+ (store-y (matrix-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-blas-dge) (y matrix-blas-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)
+ (matrix-store x)
+ (matrix-store y)
+ fftw-ffi:+FFTW-ESTIMATE+)
+ y)
+
+(defmethod fftw2-c2r! ((x matrix-blas-zge) (y matrix-blas-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)
+ (matrix-store x)
+ (matrix-store y)
+ fftw-ffi:+FFTW-ESTIMATE+)
+ y)
\ No newline at end of file
More information about the lisplab-cvs
mailing list