[lisplab-cvs] r35 - in src: fft matrix
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Fri May 22 09:02:11 UTC 2009
Author: jivestgarden
Date: Fri May 22 05:02:08 2009
New Revision: 35
Log:
optimized and cleaned up
Modified:
src/fft/level3-fft-blas.lisp
src/fft/level3-fft-interface.lisp
src/matrix/level1-classes.lisp
src/matrix/level1-interface.lisp
src/matrix/level1-matrix.lisp
src/matrix/level1-util.lisp
src/matrix/level2-matrix-dge.lisp
src/matrix/level2-matrix-zge.lisp
Modified: src/fft/level3-fft-blas.lisp
==============================================================================
--- src/fft/level3-fft-blas.lisp (original)
+++ src/fft/level3-fft-blas.lisp Fri May 22 05:02:08 2009
@@ -92,19 +92,20 @@
(complex double-float))
(setf ref-blas-complex-store2)))
-(declaim (inline ref-blas-complex-store2 (setf ref-blas-complex-store)))
+(declaim (inline ref-blas-complex-store2 (setf ref-blas-complex-store2)))
(defun ref-blas-complex-store2 (store i start step)
"Accessor for the complex blas store"
(declare (type-blas-idx i start step)
(type-blas-store store))
- (let ((idx (truly-the type-blas-idx
- (* 2 (+
- (truly-the type-blas-idx (* step i))
- start)))))
- (declare (type-blas-idx idx))
- (complex (aref store idx)
- (aref store (1+ idx)))))
+ (let* ((idx (truly-the type-blas-idx
+ (* 2 (+ (truly-the type-blas-idx (* step i))
+ start))))
+ (val (complex (aref store idx)
+ (aref store (1+ idx)))))
+ (declare (type type-blas-idx idx)
+ (type (complex double-float) val))
+ val))
(defun (setf ref-blas-complex-store2) (value store i start step)
(declare (type-blas-idx i start step)
@@ -136,10 +137,10 @@
(let* ((dual (expt 2 bit))
(W #C(1.0 0.0))
(tmp (- (exp (/ (* sign %i pi) dual)) 1.0 )))
- (declare (type-blas-idx dual)
- ((integer 0 30) bit)
- ((complex double-float) W tmp))
- (loop for b from 0 to (- n 1) by (* 2 dual) do
+ (declare (type type-blas-idx dual)
+ (type (integer 0 30) bit)
+ (type (complex double-float) W tmp))
+ (loop for b from 0 below n by (* 2 dual) do
(let* ((wd (ref-blas-complex-store2 ftx (truly-the type-blas-idx (+ b dual)) start step)))
(declare (type-blas-idx b)
((complex double-float) Wd))
@@ -162,8 +163,6 @@
wd))))))
ftx))
-
-
(defun bit-reverse-blas-complex-store! (vec n start step)
;; This is the Goldrader bit-reversal algorithm
(let ((j 0))
@@ -177,8 +176,29 @@
(setf (ref-blas-complex-store2 vec i start step) (ref-blas-complex-store2 vec j start step)
(ref-blas-complex-store2 vec j start step) tmp)))
(do () ((> k j))
- (setf j (the fixnum (- j k))
+ (setf j (the type-blas-idx (- j k))
k (floor k 2)))
(incf j k))))
vec)
+(defmethod fft-shift ((k matrix-base))
+ "Only for 2D. TODO 1d."
+ (let ((out (copy k))
+ (r/2 (/ (rows k) 2))
+ (c/2 (/ (cols k) 2)))
+ (dotimes (i (rows k))
+ (dotimes (j (cols k))
+ (setf (mref out i j)
+ (cond ((and (< i r/2) (< j c/2))
+ (mref k (+ i r/2) (+ j c/2)))
+ ((and (< i r/2) (>= j c/2))
+ (mref k (+ i r/2) (- j c/2)))
+ ((and (>= i r/2) (< j c/2))
+ (mref k (- i r/2) (+ j c/2)))
+ (t
+ (mref k (- i r/2) (- j c/2)))))))
+ out))
+
+(defmethod ifft-shift ((k matrix-base))
+ "Currently the same as fft-shift since only grids with power 2 sized grids are allowed."
+ (fft-shift k))
\ No newline at end of file
Modified: src/fft/level3-fft-interface.lisp
==============================================================================
--- src/fft/level3-fft-interface.lisp (original)
+++ src/fft/level3-fft-interface.lisp Fri May 22 05:02:08 2009
@@ -49,7 +49,7 @@
(:documentation "Inverse fast fourier transform on all rows and columns. Destructive"))
(defgeneric fft-shift (x)
- (:documentation "Resucturing of Brillouin zones"))
+ (:documentation "Restructuring of Brillouin zones"))
(defgeneric ifft-shift (x)
- (:documentation "Resucturing of Brillouin zones"))
+ (:documentation "Inverse restructuring of Brillouin zones"))
Modified: src/matrix/level1-classes.lisp
==============================================================================
--- src/matrix/level1-classes.lisp (original)
+++ src/matrix/level1-classes.lisp Fri May 22 05:02:08 2009
@@ -30,6 +30,8 @@
;; Do we need the others?
))
+(declaim (inline matrix-store))
+
(defclass matrix-base () ())
;;; The matrix element tells the element type of the matrix
@@ -69,14 +71,12 @@
:initarg :rows
:initform 0
:reader rows
- :type type-blas-idx
- :documentation "Number of rows in the matrix")
+ :type type-blas-idx)
(cols
:initarg :cols
:initform 0
:reader cols
- :type type-blas-idx
- :documentation "Number of columns in the matrix")
+ :type type-blas-idx)
(size
:reader size
:type type-blas-idx)))
Modified: src/matrix/level1-interface.lisp
==============================================================================
--- src/matrix/level1-interface.lisp (original)
+++ src/matrix/level1-interface.lisp Fri May 22 05:02:08 2009
@@ -68,7 +68,7 @@
(defgeneric (setf rank) (value matrix))
(defgeneric rows (matrix)
- (:documentation "The number of columns, ie (dim 0)."))
+ (:documentation "The number of rows, ie (dim 0)."))
(defgeneric (setf rows) (value matrix))
Modified: src/matrix/level1-matrix.lisp
==============================================================================
--- src/matrix/level1-matrix.lisp (original)
+++ src/matrix/level1-matrix.lisp Fri May 22 05:02:08 2009
@@ -66,47 +66,47 @@
;;; Spcialized for blas-dge
(defmethod mref ((matrix matrix-base-dge) row col)
- (aref (the type-blas-store (matrix-store matrix))
- (truly-the type-blas-idx (column-major-idx (truly-the type-blas-idx row)
- (truly-the type-blas-idx col)
- (rows matrix)))))
-
-(defmethod (setf mref) (value (matrix matrix-base-dge) row col)
- (setf (aref (the type-blas-store (matrix-store matrix))
- (column-major-idx (truly-the type-blas-idx row)
- (truly-the type-blas-idx col)
- (rows matrix)))
- (truly-the double-float (coerce value 'double-float))))
+ (ref-blas-real-store (matrix-store matrix) row col (rows matrix)))
+
+(defmethod (setf mref) (value (matrix matrix-base-dge) row col)
+ (let ((val2 (coerce value 'double-float)))
+ (declare (type double-float val2))
+ (setf (ref-blas-real-store (matrix-store matrix) row col (rows matrix))
+ val2)
+ val2))
(defmethod vref ((matrix matrix-base-dge) idx)
(aref (the type-blas-store (matrix-store matrix)) idx))
(defmethod (setf vref) (value (matrix matrix-base-dge) idx)
- (setf (aref (the type-blas-store (matrix-store matrix)) idx)
- (the double-float (coerce value 'double-float))))
-
+ (let ((val2 (coerce value 'double-float)))
+ (declare (type double-float val2))
+ (setf (aref (the type-blas-store (matrix-store matrix)) idx)
+ val2)
+ val2))
;;; Spcialized for blas-zge
(defmethod mref ((matrix matrix-base-zge) row col)
- (ref-blas-complex-store (matrix-store matrix)
- (column-major-idx row col (rows matrix))
- 0 1))
+ (ref-blas-complex-store (matrix-store matrix)
+ row col (rows matrix)))
(defmethod (setf mref) (value (matrix matrix-base-zge) row col)
- (setf (ref-blas-complex-store (matrix-store matrix)
- (column-major-idx row col (rows matrix))
- 0 1)
- (coerce value '(complex double-float)))
- value)
+ (let ((val2 (coerce value '(complex double-float))))
+ (declare (type (complex double-float) val2))
+ (setf (ref-blas-complex-store (matrix-store matrix) row col (rows matrix))
+ val2)
+ val2))
(defmethod vref ((matrix matrix-base-zge) i)
(ref-blas-complex-store (matrix-store matrix) i 0 1))
(defmethod (setf vref) (value (matrix matrix-base-zge) i)
- (setf (ref-blas-complex-store (matrix-store matrix) i 0 1)
- (coerce value '(complex double-float)))
- value)
+ (let ((val2 (coerce value '(complex double-float))))
+ (declare (type (complex double-float) val2))
+ (setf (ref-blas-complex-store (matrix-store matrix) i 0 1)
+ val2)
+ val2))
;;; Function matrix
Modified: src/matrix/level1-util.lisp
==============================================================================
--- src/matrix/level1-util.lisp (original)
+++ src/matrix/level1-util.lisp Fri May 22 05:02:08 2009
@@ -38,6 +38,10 @@
(deftype type-blas-idx ()
'(MOD 536870911))
+(declaim (inline column-major-idx))
+(declaim (inline ref-blas-real-store (setf ref-blas-real-store)))
+(declaim (inline ref-blas-complex-store (setf ref-blas-complex-store)))
+
(declaim (ftype (function
(type-blas-idx
type-blas-idx
@@ -45,10 +49,6 @@
type-blas-idx)
column-major-idx))
-(declaim (inline column-major-idx))
-(defun column-major-idx (i j rows)
- (truly-the type-blas-idx (+ i (truly-the type-blas-idx (* j rows)))))
-
(declaim (ftype (function
(type-blas-store
type-blas-idx
@@ -66,8 +66,26 @@
)
double-float)
(setf ref-blas-real-store)))
+(declaim (ftype (function
+ (type-blas-store
+ type-blas-idx
+ type-blas-idx
+ type-blas-idx)
+ (complex double-float))
+ ref-blas-complex-store))
-(declaim (inline ref-blas-real-store (setf ref-blas-real-store)))
+(declaim (ftype (function
+ ((complex double-float)
+ type-blas-store
+ type-blas-idx
+ type-blas-idx
+ type-blas-idx
+ )
+ (complex double-float))
+ (setf ref-blas-complex-store)))
+
+(defun column-major-idx (i j rows)
+ (truly-the type-blas-idx (+ i (truly-the type-blas-idx (* j rows)))))
(defun ref-blas-real-store (store row col rows)
"Accessor for the real blas store"
@@ -90,27 +108,6 @@
:initial-element
(coerce initial-element 'double-float)))
-
-(declaim (ftype (function
- (type-blas-store
- type-blas-idx
- type-blas-idx
- type-blas-idx)
- (complex double-float))
- ref-blas-complex-store))
-
-(declaim (ftype (function
- ((complex double-float)
- type-blas-store
- type-blas-idx
- type-blas-idx
- type-blas-idx
- )
- (complex double-float))
- (setf ref-blas-complex-store)))
-
-(declaim (inline ref-blas-complex-store (setf ref-blas-complex-store)))
-
(defun ref-blas-complex-store (store row col rows)
"Accessor for the complet blas store"
(let ((idx (truly-the type-blas-idx
@@ -136,6 +133,7 @@
(rv (coerce (realpart value) 'double-float))
(iv (coerce (imagpart value) 'double-float))
(store (allocate-real-store 2size iv)))
+ (declare (type type-blas-idx 2size))
(loop for i from 0 below 2size by 2 do
(setf (aref store i) rv))
store))
Modified: src/matrix/level2-matrix-dge.lisp
==============================================================================
--- src/matrix/level2-matrix-dge.lisp (original)
+++ src/matrix/level2-matrix-dge.lisp Fri May 22 05:02:08 2009
@@ -22,11 +22,12 @@
(defmethod fill! ((a matrix-lisp-dge) value)
(let ((x (coerce value 'double-float))
(store (matrix-store a)))
+ (declare (type type-blas-store store))
(fill store x)))
(defmethod copy ((matrix matrix-lisp-dge))
(make-instance (class-name (class-of matrix))
- :store (copy-seq (matrix-store matrix))
+ :store (copy-seq (the type-blas-store (matrix-store matrix)))
:rows (rows matrix)
:cols (cols matrix)))
Modified: src/matrix/level2-matrix-zge.lisp
==============================================================================
--- src/matrix/level2-matrix-zge.lisp (original)
+++ src/matrix/level2-matrix-zge.lisp Fri May 22 05:02:08 2009
@@ -26,10 +26,9 @@
(setf (aref store i) rx
(aref store (1+ i)) cx))))
-
-(defmethod copy ((matrix matrix-base-zge))
+(defmethod copy ((matrix matrix-lisp-zge))
(make-instance (class-name (class-of matrix))
- :store (copy-seq (matrix-store matrix))
+ :store (copy-seq (the type-blas-store (matrix-store matrix)))
:rows (rows matrix)
:cols (cols matrix)))
@@ -51,10 +50,6 @@
(copy-contents a b #'abs)
b))
-
-
-
-
(defmacro def-binary-op-blas-complex (new old)
;;; TODO speed up for real numbers
(let ((a (gensym "a"))
More information about the lisplab-cvs
mailing list