[lisplab-cvs] r45 - src/fft src/matlisp system
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Wed May 27 17:56:15 UTC 2009
Author: jivestgarden
Date: Wed May 27 13:56:14 2009
New Revision: 45
Log:
added ffiw lib. Not yet tested
Added:
src/fft/fftw-ffi-package.lisp
src/fft/fftw-ffi.lisp
src/fft/level3-fft-fftw.lisp
src/matlisp/f77-package.lisp
Modified:
src/fft/level3-fft-blas.lisp
src/matlisp/geev.lisp
system/lisplab.asd
system/package.lisp
Added: src/fft/fftw-ffi-package.lisp
==============================================================================
--- (empty file)
+++ src/fft/fftw-ffi-package.lisp Wed May 27 13:56:14 2009
@@ -0,0 +1,25 @@
+;;; 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.
+
+(defpackage "FFTW-FFI"
+ (:use "COMMON-LISP" "SB-ALIEN" "SB-SYS")
+ (:export "+DOUBLE_FLOAT-BYTES+"
+ "+FFTW-ESTIMATE+"
+ "+FFTW-FORWARD+"
+ "+FFTW-BACKWARD+"
+ "FFTW-FFT1"
+ "FFTW-FFT2")
+ (:documentation "Simple ffi for fftw."))
Added: src/fft/fftw-ffi.lisp
==============================================================================
--- (empty file)
+++ src/fft/fftw-ffi.lisp Wed May 27 13:56:14 2009
@@ -0,0 +1,89 @@
+;;; Foreign function interfaces for 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.
+
+(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)
+
+(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)))
+
+(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+)))
+ (without-gcing
+ (let ((plan (|fftw_plan_dft_1d|
+ n
+ (sap+ (vector-sap a) astart)
+ (sap+ (vector-sap b) bstart)
+ direction
+ flag)))
+ (|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."
+ ;; TODO we should handle conditions to avoid mem-leaks
+ (without-gcing
+ (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)))
+ (|fftw_execute| plan)
+ (|fftw_destroy_plan| plan)))
+ out)
+
+
Modified: src/fft/level3-fft-blas.lisp
==============================================================================
--- src/fft/level3-fft-blas.lisp (original)
+++ src/fft/level3-fft-blas.lisp Wed May 27 13:56:14 2009
@@ -52,6 +52,8 @@
(defmethod fft2 ((x matrix-lisp-zge))
(fft2! (copy x)))
+;;;; The implementing methods
+
(defmethod fft1! ((x matrix-lisp-zge))
(dotimes (i (cols x))
(fft-radix-2-blas-complex-store! :f (matrix-store x) (rows x) (* (rows x) i) 1))
Added: src/fft/level3-fft-fftw.lisp
==============================================================================
--- (empty file)
+++ src/fft/level3-fft-fftw.lisp Wed May 27 13:56:14 2009
@@ -0,0 +1,59 @@
+;;; 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.
+
+(in-package :lisplab)
+
+(defun fft1!-forward-or-backward (x direction)
+ (let* ((rows (rows x))
+ (cols (cols x))
+ (store (matrix-store x)))
+ (print "allois")
+ (dotimes (i cols)
+ (fftw-ffi:fftw-fft1
+ rows
+ store
+ (* i cols)
+ store
+ (* i cols)
+ direction
+ fftw-ffi:+FFTW-ESTIMATE+)))
+ x)
+
+(defmethod fft1! ((x matrix-blas-zge))
+ (fft1!-forward-or-backward x fftw-ffi:+fftw-forward+))
+
+(defmethod ifft1! ((x matrix-blas-zge))
+ (fft1!-forward-or-backward x fftw-ffi:+fftw-backward+))
+
+(defmethod fft2! ((x matrix-blas-zge))
+ (fftw-ffi:fftw-fft2
+ (rows x)
+ (cols x)
+ (matrix-store x)
+ (matrix-store x)
+ fftw-ffi:+fftw-forward+
+ fftw-ffi:+FFTW-ESTIMATE+))
+
+(defmethod fft2! ((x matrix-blas-zge))
+ (fftw-ffi:fftw-fft2
+ (rows x)
+ (cols x)
+ (matrix-store x)
+ (matrix-store x)
+ fftw-ffi:+fftw-backward+
+ fftw-ffi:+FFTW-ESTIMATE+))
Added: src/matlisp/f77-package.lisp
==============================================================================
--- (empty file)
+++ src/matlisp/f77-package.lisp Wed May 27 13:56:14 2009
@@ -0,0 +1,5 @@
+
+(defpackage "F77"
+ (:use "COMMON-LISP" "SB-EXT" "FORTRAN-FFI-ACCESSORS")
+ ;; I'm too lazy to export. Just force way to get the functions.
+ (:documentation "Wrappers for Fortran functions"))
\ No newline at end of file
Modified: src/matlisp/geev.lisp
==============================================================================
--- src/matlisp/geev.lisp (original)
+++ src/matlisp/geev.lisp Wed May 27 13:56:14 2009
@@ -49,7 +49,7 @@
(defmethod rearrange-eigenvector-matrix ((evals matrix-blas-zge) (p matrix-blas-dge))
(let* ((n (size evals))
- (evec (cnew 0 n n)))
+ (evec (znew 0 n n))) ; TODO aggregate input type
(do ((col 0 (incf col)))
((>= col n))
(if (zerop (imagpart (vref evals col )))
@@ -165,7 +165,7 @@
(let* ((n (rows a))
(2n (* 2 n))
(xxx (allocate-real-store 2))
- (w (cnew 0 n 1))
+ (w (znew 0 n 1)) ; TODO aggregate type
(vl (if vl-mat (matrix-store vl-mat) xxx))
(vr (if vr-mat (matrix-store vr-mat) xxx))
(lwork (zgeev-workspace-size n (if vl-mat t nil) (if vr-mat t nil)))
Modified: system/lisplab.asd
==============================================================================
--- system/lisplab.asd (original)
+++ system/lisplab.asd Wed May 27 13:56:14 2009
@@ -121,7 +121,7 @@
:pathname "../src/matlisp/"
:serial t
:components
- (
+ ((:file "f77-package")
(:file "f77-mangling")
(:file "ffi-sbcl")
(:file "blas")
@@ -138,6 +138,21 @@
(dolist (lib asdf::*lisplab-external-libraries*)
(format t "Loads alien object <~A>" lib))))
+
+ ;;
+ ;; Blas and Lapack implmentations (Level 3)
+ ;;
+ (:module :fftw
+ :depends-on (:matrix :fft)
+ :pathname "../src/fft/"
+ :serial t
+ :components
+ ((:file "fftw-ffi-package")
+ (:file "fftw-ffi")
+ (:file "level3-fft-fftw"))
+ :perform (asdf:load-op :after (op c)
+ (sb-alien:load-shared-object "/usr/lib/libfftw3.so.3")))
+
;;
;; Euler and Runge-Kutt solvers (Level 3)
;;
Modified: system/package.lisp
==============================================================================
--- system/package.lisp (original)
+++ system/package.lisp Wed May 27 13:56:14 2009
@@ -17,6 +17,8 @@
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+;;; TODO split this out
+
(defpackage "LISPLAB"
(:nicknames "LL")
(:use "COMMON-LISP" "SB-EXT")
@@ -29,10 +31,6 @@
"INCF-SAP"
"WITH-VECTOR-DATA-ADDRESSES")
(:documentation "Fortran foreign function interface"))
-
-(defpackage "F77"
- (:use "COMMON-LISP" "SB-EXT" "FORTRAN-FFI-ACCESSORS")
- (:documentation "Wrappers for Fortran functions"))
(defpackage "LISPLAB-BLAS"
(:nicknames "LL-BLAS")
More information about the lisplab-cvs
mailing list