[lisplab-cvs] r121 - in src: fft matrix

[Jørn Inge Vestgården]

Author: jivestgarden
Date: Sun Dec 13 09:19:15 2009
New Revision: 121

Log:
optimized constuctores and fftshift

Modified:
   src/fft/level3-fft-zge.lisp
   src/matrix/level2-constructors.lisp

Modified: src/fft/level3-fft-zge.lisp
==============================================================================
--- src/fft/level3-fft-zge.lisp	(original)
+++ src/fft/level3-fft-zge.lisp	Sun Dec 13 09:19:15 2009
@@ -19,8 +19,53 @@
 
 ;;; TODO should use the normal ref-blas-complex-store ?
 
+;;; TODO make shift routines for complex matrices also.
+
 (in-package :lisplab)
 
+;; Optimized shift methods for real matrix. 
+(defmethod fft-shift ((m matrix-base-dge))
+  (let* ((rows (rows m))
+         (fr (floor rows 2))
+         (cr (ceiling rows 2))
+         (cols (cols m))
+         (fc (floor cols 2))
+         (cc (ceiling cols 2))
+         (store-m (matrix-store m))
+         (m2 (mcreate m))
+         (store-m2 (matrix-store m2)))
+         
+    (declare (type type-blas-store store-m store-m2)
+             (type type-blas-idx fr cr fc cc rows cols))
+    (dotimes (i rows)
+      (dotimes (j cols)
+        (setf (aref store-m2 (column-major-idx i j rows))
+              (aref store-m (column-major-idx  (if (< i fr) (+ i cr) (- i fr))
+                                               (if (< j fc) (+ j cc) (- j fc))
+                                               rows)))))
+    m2))
+
+(defmethod ifft-shift ((m matrix-base-dge))
+  (let* ((rows (rows m))
+         (fr (floor rows 2))
+         (cr (ceiling rows 2))
+         (cols (cols m))
+         (fc (floor cols 2))
+         (cc (ceiling cols 2))
+                 (store-m (matrix-store m))
+         (m2 (mcreate m))
+         (store-m2 (matrix-store m2)))
+         
+    (declare (type type-blas-store store-m store-m2)
+             (type type-blas-idx fr cr fc cc rows cols))
+    (dotimes (i rows)
+      (dotimes (j cols)
+        (setf (aref store-m2 (column-major-idx i j rows))
+              (aref store-m (column-major-idx  (if (< i cr) (+ i fr) (- i cr))
+                                               (if (< j cc) (+ j fc) (- j cc))
+                                               rows)))))
+    m2))
+
 ;;;; The implementing methods
 
 (defmethod fft1! ((x matrix-lisp-zge) &key)
@@ -119,10 +164,15 @@
 		 (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)))
+	      (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))
-		(setf (ref-blas-complex-store2 ftx (truly-the type-blas-idx (+ b dual)) start step) 
+		(setf (ref-blas-complex-store2 ftx (truly-the type-blas-idx (+ b dual)) 
+					       start 
+					       step) 
 		      (- (ref-blas-complex-store2 ftx b start step) wd))
 		(incf (ref-blas-complex-store2 ftx b start step) 
 		      wd)))   
@@ -151,8 +201,10 @@
 	(declare (type-blas-idx i k))
 	(when (< i j)
 	  (let ((tmp (ref-blas-complex-store2 vec i start step)))
-	    (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)))
+	    (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 type-blas-idx (- j k))
 		k (floor k 2)))

Modified: src/matrix/level2-constructors.lisp
==============================================================================
--- src/matrix/level2-constructors.lisp	(original)
+++ src/matrix/level2-constructors.lisp	Sun Dec 13 09:19:15 2009
@@ -113,27 +113,41 @@
 
 For example: (drange 4 0 1) -> 0 1 2 3, while 
              (drange 4 0 1 0.5) -> 0.5 1.5 2.5 3.5."
-  (let ((x (dnew 0 n 1))
-	(dx (./ (.- to from) N)))
-    (dotimes (i n)
-      (setf (vref x i) (.+ from (.* dx (.+ i shift)))))
-    x))
+  (let* ((x (dnew 0 n 1))
+	 (store (matrix-store x))
+	 (dx (/ (- (to-df to) (to-df from)) (to-df N)))
+	 (shift (* dx (to-df shift))))
+    (declare (type type-blas-store store)
+	     (type type-blas-idx n)
+	     (type double-float dx shift))
+    (do ((i 0 (1+ i))
+	 (v (to-df from) (+ v dx)))
+	((>= i n) x)
+      (declare (type double-float v))
+      (setf (aref store i) (+ shift v)))))
 
 (defun dgrid (xv yv)
   "Creates grid matrices from input vectors. Input are the x and y vectors 
 and outputs are a list of x and y matrices. The input vectors are  
 typically created with drange."
-  (let* ((r (size xv))
-	 (c (size yv))
-	 (x (dnew 0 r c))
-	 (y (dnew 0 r c)))
+   (let* ((r (size xv))
+	  (c (size yv))
+	  (x (dnew 0 r c))
+	  (y (dnew 0 r c))	  
+	  (xv* (matrix-store xv))
+	  (yv* (matrix-store yv))
+	  (x* (matrix-store x))
+	  (y* (matrix-store y)))
+    (declare (type type-blas-store xv* yv* x* y*)
+             (type type-blas-idx r c))
     (dotimes (i r)
-      (dotimes (j c)
-	(setf (mref x i j) (vref xv i)
-	      (mref y i j) (vref yv j))))
+      (dotimes (j c)    
+        (let ((k (column-major-idx i j r)))
+          (declare (type type-blas-idx k))
+          (setf (aref x* k) (aref xv* i)
+                (aref y* k) (aref yv* j)))))
     (list x y)))
 
-
 ;;; Constructors for matrix-zge
 
 (defmacro zmat (&body args)