[lisplab-cvs] r118 - in src: core matrix
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Sat Dec 12 19:31:02 UTC 2009
Author: jivestgarden
Date: Sat Dec 12 14:31:02 2009
New Revision: 118
Log:
optimized and cleaned up functions and operators
Modified:
src/core/level0-default.lisp
src/matrix/level2-function.lisp
src/matrix/level2-matrix-dge.lisp
src/matrix/level2-matrix-zge.lisp
src/matrix/store-operators.lisp
src/matrix/store-ordinary-functions.lisp
Modified: src/core/level0-default.lisp
==============================================================================
--- src/core/level0-default.lisp (original)
+++ src/core/level0-default.lisp Sat Dec 12 14:31:02 2009
@@ -42,6 +42,12 @@
;;; Some known exepctions
+(defmethod function-output-type-spec ((fun (eql '.ln)) (input-spec (eql :d)))
+ :z)
+
+(defmethod function-output-type-spec ((fun (eql '.sqrt)) (input-spec (eql :d)))
+ :z)
+
(defmethod function-output-type-spec ((fun (eql '.asin)) (input-spec (eql :d)))
:z)
Modified: src/matrix/level2-function.lisp
==============================================================================
--- src/matrix/level2-function.lisp (original)
+++ src/matrix/level2-function.lisp Sat Dec 12 14:31:02 2009
@@ -40,10 +40,8 @@
`(def-each-element-function ,name))
+ordinary-functions-number-to-number-list+ )))
-
(expand-each-element-ordinary-functions)
-
;;; Some special functions. Should maybe be separated out.
@@ -103,53 +101,3 @@
a))
-
-#|
-
-
-(defmacro each-element-function-matrix-ge (x form)
- "Applies a form on each element of an matrix-ge."
- (let ((i (gensym))
- (y (gensym)))
- `(let* ((,y (copy ,x)))
- (dotimes (,i (size ,y))
- (let ((,x (vref ,y ,i)))
- (setf (vref ,y ,i)
- ,form)))
- ,y)))
-
-(defmacro expand-matrix-ge-num-num ()
- (cons 'progn
- (mapcar (lambda (name)
- ;; Note: not using the (cdr name) , which is only valid
- ;; for build in lisp types.
- `(defmethod ,(car name) ((x matrix-ge))
- (each-element-function-matrix-ge x (,(car name) x))))
- +functions-real-to-real+)))
-
-(expand-matrix-ge-num-num)
-
-(defmethod .log ((x matrix-ge) &optional base)
- (each-element-function-matrix-ge x (.log x base)))
-
-;;; Bessel functions
-
-(defmethod .besj (n (x matrix-ge))
- (each-element-function-matrix-ge x (.besj n x)))
-
-(defmethod .besy (n (x matrix-ge))
- (each-element-function-matrix-ge x (.besy n x)))
-
-(defmethod .besi (n (x matrix-ge))
- (each-element-function-matrix-ge x (.besi n x)))
-
-(defmethod .besk (n (x matrix-ge))
- (each-element-function-matrix-ge x (.besk n x)))
-
-(defmethod .besh1 (n (x matrix-ge))
- (each-element-function-matrix-ge x (.besh1 n x)))
-
-(defmethod .besh2 (n (x matrix-ge))
- (each-element-function-matrix-ge x (.besh2 n x)))
-
-|#
\ No newline at end of file
Modified: src/matrix/level2-matrix-dge.lisp
==============================================================================
--- src/matrix/level2-matrix-dge.lisp (original)
+++ src/matrix/level2-matrix-dge.lisp Sat Dec 12 14:31:02 2009
@@ -26,13 +26,16 @@
(fill store x)))
(defmethod copy ((matrix matrix-base-dge))
- (make-instance (class-name (class-of matrix))
- :store (copy-seq (the type-blas-store (matrix-store matrix)))
- :dim (dim matrix)))
+ (let ((store (matrix-store matrix)))
+ (declare (type type-blas-store store))
+ (make-instance (class-name (class-of matrix))
+ :store store
+ :dim (dim matrix))))
(defmethod copy-contents ((a matrix-base-dge) (b matrix-base-dge) &optional (converter nil))
(let ((store-a (matrix-store a))
(store-b (matrix-store b)))
+ (declare (type type-blas-store store-a store-b))
(if converter
(map-into store-b converter store-a)
(copy-matrix-stores store-a store-b)))
@@ -47,12 +50,18 @@
out)
(defmethod mmap ((type (eql nil)) f (a matrix-base-dge) &rest args)
- (apply #'map
- nil
- (lambda (&rest args)
- (coerce (apply f args) 'double-float))
- (matrix-store a) (mapcar #'matrix-store args))
- nil)
+ "No output. Called for side effects."
+ (declare (type (function (double-float) t) f))
+ (if args
+ (apply #'map
+ nil
+ (lambda (&rest args)
+ (apply f args))
+ (matrix-store a) (mapcar #'matrix-store args))
+ (let ((store (matrix-store a)))
+ (declare (type type-blas-store store))
+ (map nil f store)))
+ nil)
(defmethod msum ((m matrix-base-dge))
(let ((sum 0.0))
@@ -153,17 +162,21 @@
(expand-generic-function-dfa-dfa-map)
-;;; The ordinary functions
-;;;; Ordinary functions
+;;;; Ordinary functions that map real to real
(define-constant +generic-function-dfa-to-dfa-map+ ;really bad name
- '((.sin . sin_dfa) (.cos . cos_dfa) (.tan . tan_dfa)
- (.atan . atan_dfa)
- (.sinh . sinh_dfa) (.cosh . cosh_dfa) (.tanh . tanh_dfa)
- (.asinh . asinh_dfa) (.acosh . acosh_dfa)
- (.exp . exp_dfa) (.ln . log_dfa) (.sqrt . sqrt_dfat) (.conjugate . conjugate_dfa)
- (.re . realpart_dfa) (.im . imagpart_dfa) (.abs . abs_dfa)))
+ '((.sin . sin_dfa-to-dfa) (.cos . cos_dfa-to-dfa) (.tan . tan_dfa-to-dfa)
+ (.atan . atan_dfa-to-dfa)
+ (.sinh . sinh_dfa-to-dfa) (.cosh . cosh_dfa-to-dfa) (.tanh . tanh_dfa-to-dfa)
+ (.asinh . asinh_dfa-to-dfa) (.acosh . acosh_dfa-to-dfa)
+ (.exp . exp_dfa-to-dfa)
+ #+nil (.ln . log_dfa)
+ #+nil (.sqrt . sqrt_dfat)
+ #+nil (.conj . conjugate_dfa)
+ #+nil (.re . realpart_dfa)
+ #+nil (.im . imagpart_dfa)
+ (.abs . abs_dfa-to-dfa)))
(defmacro defmethod-dfa-to-dfa (name underlying-function)
(let ((a (gensym "a"))
@@ -181,144 +194,42 @@
(expand-generic-function-dfa-to-dfa-map)
-;;; The rest must wait until tomorrow
-
-
-
-
-
-;;;; Old code
-
-#|
-
-(defmacro each-matrix-element-df-to-df (x form)
- "Applies a form on each element of an matrix-dge. The form must
-make real output for real arguments"
- (let ((i (gensym))
- (store (gensym)))
- `(let* ((,x (copy ,x))
- (,store (matrix-store ,x)))
- (declare (type type-blas-store ,store))
- (dotimes (,i (length ,store))
- (let ((,x (aref ,store ,i)))
- (declare (type type-blas-idx ,i)
- (type double-float ,x))
- (setf (aref ,store ,i)
- ,form)))
- ,x)))
-
-(defmacro expand-matrix-dge-num-num ()
- (cons 'progn
- (mapcar (lambda (name)
- `(defmethod ,(car name) ((x matrix-base-dge))
- (each-matrix-element-df-to-df x (,(cdr name) x))))
- +functions-real-to-real+)))
-
-(expand-matrix-dge-num-num)
-
-;;; Bessel functions
+;;; Some trivialities
-(defmethod .besj (n (x matrix-base-dge))
- (each-matrix-element-df-to-df x (.besj n x)))
-
-(defmethod .besy (n (x matrix-base-dge))
- (each-matrix-element-df-to-df x (.besy n x)))
+(defmethod .re ((a matrix-base-dge))
+ (copy a))
-(defmethod .besi (n (x matrix-base-dge))
- (each-matrix-element-df-to-df x (.besi n x)))
+(defmethod .im ((a matrix-base-dge))
+ (mcreate a))
-(defmethod .besk (n (x matrix-base-dge))
- (each-matrix-element-df-to-df x (.besk n x)))
+(defmethod .conj ((a matrix-base-dge))
+ (copy a))
-(defmacro each-matrix-element-df-to-complex-df (x form)
- "Applies a form on each element of an matrix-dge. The form must
-make complex output for real arguments. TODO optimize? Probably no need. The
-Hankel functions are slow anyway."
- (let ((i (gensym))
- (b (gensym))
- (spec-b (gensym)))
- `(let* ((,spec-b (create-matrix-description ,x :et :z))
- (,b (convert ,x ,spec-b) ))
- (dotimes (,i (size ,x))
- (let ((,x (vref ,x ,i)))
- (setf (vref ,b ,i) ,form)))
- ,b)))
+;;; Now these sad cases where output might be complex for real input
-(defmethod .asin ((x matrix-base-dge))
- (each-matrix-element-df-to-complex-df x (asin x)))
+(define-constant +generic-function-dfa-to-cdfa-map+ ;really bad name
+ '((.sqrt . sqrt_dfa)
+ (.ln . log_dfa)
+ (.asin . asin_dfa)
+ (.acos . acos_dfa)
+ (.atanh . atanh_dfa)
+ ;; Some more?
+ ))
-(defmethod .acos ((x matrix-base-dge))
- (each-matrix-element-df-to-complex-df x (asin x)))
-
-(defmethod .atanh ((x matrix-base-dge))
- (each-matrix-element-df-to-complex-df x (asin x)))
-
-(defmethod .besh1 (n (x matrix-base-dge))
- (each-matrix-element-df-to-complex-df x (.besh1 n x)))
-
-(defmethod .besh2 (n (x matrix-base-dge))
- (each-matrix-element-df-to-complex-df x (.besh2 n x)))
-
-|#
-
-#|
-
-;;; Old code
-
-(defmacro def-binary-op-matrix-base-dge (new old)
+(defmacro defmethod-dfa-to-cdfa (name underlying-function)
(let ((a (gensym "a"))
(b (gensym "b"))
- (len (gensym "len"))
- (store (gensym "store"))
- (store2 (gensym "store2"))
- (i (gensym "i")))
- `(progn
- (defmethod ,new ((,a matrix-base-dge) (,b real))
- (let* ((,a (copy ,a))
- (,store (matrix-store ,a))
- (,b (coerce ,b 'double-float))
- (,len (size ,a)))
- (declare (type double-float ,b)
- (type type-blas-store ,store)
- (type type-blas-idx ,len))
- (dotimes (,i ,len)
- (setf (aref ,store ,i) (,old (aref ,store ,i) ,b)))
- ,a))
- (defmethod ,new ((,a real) (,b matrix-base-dge))
- (let* ((,b (copy ,b))
- (,store (matrix-store ,b))
- (,a (coerce ,a 'double-float))
- (,len (size ,b)))
- (declare (type double-float ,a)
- (type type-blas-store ,store)
- (type type-blas-idx ,len))
- (dotimes (,i ,len)
- (setf (aref ,store ,i) (,old ,a (aref ,store ,i))))
- ,b))
- (defmethod ,new ((,a matrix-base-dge) (,b matrix-base-dge))
- (let* ((,a (copy ,a))
- (,store (matrix-store ,a))
- (,store2 (matrix-store ,b))
- (,len (size ,a)))
- (declare (type type-blas-store ,store)
- (type type-blas-store ,store2)
- (type type-blas-idx ,len))
- (dotimes (,i ,len)
- (setf (aref ,store ,i) (,old (aref ,store ,i) (aref ,store2 ,i))))
- ,a)))))
-
-(def-binary-op-matrix-base-dge .add +)
-
-(def-binary-op-matrix-base-dge .sub -)
-
-(def-binary-op-matrix-base-dge .mul *)
-
-(def-binary-op-matrix-base-dge .div /)
-
-(def-binary-op-matrix-base-dge .expt expt)
-
-(def-binary-op-matrix-base-dge .min min)
+ (spec (gensym "spec")))
+ `(defmethod ,name ((,a matrix-base-dge))
+ (let* ((,spec (cons :z (cdr (type-spec ,a))))
+ (,b (make-matrix-instance ,spec (dim ,a) 0)))
+ (,underlying-function (matrix-store ,a) (matrix-store ,b) )
+ ,b))))
-(def-binary-op-matrix-base-dge .max max)
+(defmacro expand-generic-function-dfa-to-cdfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-dfa-to-cdfa ,(car x) ,(cdr x)))
+ +generic-function-dfa-to-cdfa-map+)))
-|#
\ No newline at end of file
+(expand-generic-function-dfa-to-cdfa-map)
Modified: src/matrix/level2-matrix-zge.lisp
==============================================================================
--- src/matrix/level2-matrix-zge.lisp (original)
+++ src/matrix/level2-matrix-zge.lisp Sat Dec 12 14:31:02 2009
@@ -50,7 +50,7 @@
(declare (type type-blas-store store-a store-b)
(type type-blas-idx len))
(dotimes (i len)
- (setf (aref store-b (* 2 i)) (aref store-a i)))
+ (setf (aref store-b (truly-the type-blas-idx (* 2 i))) (aref store-a i)))
to)))
(defmethod msum ((m matrix-base-zge))
@@ -263,7 +263,7 @@
(.asin . asin_cdfa) (.acos . acos_cdfa) (.atan . atan_cdfa)
(.sinh . sinh_cdfa) (.cosh . cosh_cdfa) (.tanh . tanh_cdfa)
(.asinh . asinh_cdfa) (.acosh . acosh_cdfa) (.atanh . atanh_cdfa)
- (.exp . exp_cdfa) (.ln . log_cdfa) (.sqrt . sqrt_cdfat) (.conjugate . conjugate_cdfa)))
+ (.exp . exp_cdfa) (.ln . log_cdfa) (.sqrt . sqrt_cdfat) (.conj . conjugate_cdfa)))
(defmacro defmethod-cdfa-to-cdfa (name underlying-function)
(let ((a (gensym "a"))
@@ -284,26 +284,20 @@
;;;; Exceptions, in that output is real for complex input
(defmethod .im ((a matrix-base-zge))
- (let ((out (make-matrix-instance
- (function-output-type-spec '.im (type-spec a))
- (dim a)
- 0)))
+ (let* ((spec-out (cons :d (cdr (type-spec a))))
+ (out (make-matrix-instance spec-out (dim a) 0)))
(imagpart_cdfa (matrix-store a) (matrix-store out))
out))
(defmethod .re ((a matrix-base-zge))
- (let ((out (make-matrix-instance
- (function-output-type-spec '.re (type-spec a))
- (dim a)
- 0)))
+ (let* ((spec-out (cons :d (cdr (type-spec a))))
+ (out (make-matrix-instance spec-out (dim a) 0)))
(realpart_cdfa (matrix-store a) (matrix-store out))
out))
(defmethod .abs ((a matrix-base-zge))
- (let ((out (make-matrix-instance
- (function-output-type-spec '.abs (type-spec a))
- (dim a)
- 0)))
+ (let* ((spec-out (cons :d (cdr (type-spec a))))
+ (out (make-matrix-instance spec-out (dim a) 0)))
(abs_cdfa (matrix-store a) (matrix-store out))
out))
Modified: src/matrix/store-operators.lisp
==============================================================================
--- src/matrix/store-operators.lisp (original)
+++ src/matrix/store-operators.lisp Sat Dec 12 14:31:02 2009
@@ -33,10 +33,11 @@
;;; TODO: there must be some easier way to generate the code in this file,
;;; but I have not the energy to do it. I do, however, think that
-;;; the basic idea of having a layer of ordinary functions is correct.
+;;; the basic idea of having a layer of ordinary functions is a good one.
;;; The reason for generating ordinary functions and not using methods,
-;;; is that the real and complex stores have the same type!
+;;; is that the real and complex stores have the same type. The fortran-compatible
+;;; complex arrays are just subsequent real and complex double-floats.
;;; The reason for having both real and complex in the same file is that
;;; not all operators function on both real and complex arguments. Care must
@@ -48,8 +49,126 @@
;;; They use a naming conventions, which should be pretty easy to
;;; guess, such as df = double float and cdfa = complex double float array.
;;;
-;;; (The last one should for consistnt naming be changed to zdfa, but its not really important)
+;;; (The convention for complex should for consistnt naming be changed to zdfa,
+;;; but its not really important)
+;;;
+;;; I use map-into when its performance is equal or better to the iterations.
+;;; The iterative version for all operations are still in the file, since other lisps
+;;; than sbcl might have a slower map-into, so that they can be needed later.
+;;; For real numbers, map-into can be used for all operations, while for complex
+;;; number only + and - (*, / and expt mix the real and complex parts)
+
+;;; The below operations are based on map-into. It is hope that some
+;;; clever lisp machine can be very fast on them (On my SBCL 32 they are exactly
+;;; the same speed as the iterative versions
+(defun +_dfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (map-into c #'+ a b))
+
+(defun +_dfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (map-into c (lambda (x) (+ x b)) a))
+
+(defun +_df-dfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (map-into c (lambda (x) (+ a x)) b))
+
+(defun -_dfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (map-into c #'- a b))
+
+(defun -_dfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (map-into c (lambda (x) (- x b)) a))
+
+(defun -_df-dfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (map-into c (lambda (x) (- a x)) b))
+
+(defun *_dfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (map-into c #'* a b))
+
+(defun *_dfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (map-into c (lambda (x) (* x b)) a))
+
+(defun *_df-dfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (map-into c (lambda (x) (* a x)) b))
+
+(defun /_dfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (map-into c #'/ a b))
+
+(defun /_dfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (map-into c (lambda (x) (/ x b)) a))
+
+(defun /_df-dfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (map-into c (lambda (x) (/ a x)) b))
+
+(defun ^_dfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (map-into c #'expt a b))
+
+(defun ^_dfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (map-into c (lambda (x) (expt x b)) a))
+
+(defun ^_df-dfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (map-into c (lambda (x) (expt a x)) b))
+
+(defun max_dfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (map-into c #'max a b))
+
+(defun max_dfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (map-into c (lambda (x) (max x b)) a))
+
+(defun max_df-dfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (map-into c (lambda (x) (max a x)) b))
+
+(defun min_dfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (map-into c #'min a b))
+
+(defun min_dfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (map-into c (lambda (x) (min x b)) a))
+
+(defun min_df-dfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (map-into c (lambda (x) (min a x)) b))
+
+
+
+;;; Complex map-into-operations
+
+(defun +_cdfa-cdfa (a b c)
+ (+_dfa-dfa a b c))
+
+(defun -_cdfa-cdfa (a b c)
+ (-_dfa-dfa a b c))
;;; Array and number
@@ -82,7 +201,7 @@
`(defun-dfa-df ,(cdr x) ,(car x)))
+operators-dfa-df-map+)))
-(expand-operators-dfa-df-map)
+#+nil (expand-operators-dfa-df-map) ; These are the iterative versions
;;; The three parts of code below contains some common strucutre that could
;;; in principle be joined, and there is also some clumsy code,
@@ -118,7 +237,7 @@
`(defun-df-dfa ,(cdr x) ,(car x)))
+operators-df-dfa-map+)))
-(expand-operators-df-dfa-map)
+#+nil (expand-operators-df-dfa-map) ; These are the iterative versions
;;; Array and array
@@ -149,7 +268,7 @@
`(defun-dfa-dfa ,(cdr x) ,(car x)))
+operators-dfa-dfa-map+)))
-(expand-operators-dfa-dfa-map)
+#+nil (expand-operators-dfa-dfa-map) ; These are the iterative versions
;;; Now the complex operators
@@ -158,8 +277,8 @@
;;; Array and number
(define-constant +operators-cdfa-cdf-map+
- '((+ . +_cdfa-cdf)
- (- . -_cdfa-cdf)
+ '((+ . +_cdfa-cdf) ; iterative version
+ (- . -_cdfa-cdf) ; iterative version
(* . *_cdfa-cdf)
(/ . /_cdfa-cdf)
(expt . ^_cdfa-cdf)
@@ -222,8 +341,8 @@
;;; Array and array
(define-constant +operators-cdfa-cdfa-map+
- '((+ . +_cdfa-cdfa)
- (- . -_cdfa-cdfa)
+ '(; (+ . +_cdfa-cdfa)
+ ; (- . -_cdfa-cdfa)
(* . *_cdfa-cdfa)
(/ . /_cdfa-cdfa)
(expt . ^_cdfa-cdfa)))
Modified: src/matrix/store-ordinary-functions.lisp
==============================================================================
--- src/matrix/store-ordinary-functions.lisp (original)
+++ src/matrix/store-ordinary-functions.lisp Sat Dec 12 14:31:02 2009
@@ -19,8 +19,6 @@
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-;;; TODO: change name of this to something about blas store
-;;;
;;; This file contains manipulations of simple double-float arrays
;;; and should be called by the spesialized matrix methods.
;;; The purpose of this layer is that it can be used by
@@ -28,6 +26,14 @@
;;;
;;; The content of this file must be highly optimized
;;; and should not depend anything exept Common Lisp itself.
+;;;
+;;; I use map-into for the real functions (impossible for the complex functions)
+;;; but keep the iterative versions still in the file since future versions
+;;; of lisplab might use them.
+
+;;; Generate more real-to-real functions. With some kind of input these will
+;;; fail and give complex output but for speed it can be ok to have them
+
(in-package :lisplab)
@@ -36,14 +42,18 @@
;;; Double-float to double float
(define-constant +ordinary-functions-dfa-to-dfa-map+
- '((sin . sin_dfa) (cos . cos_dfa) (tan . tan_dfa)
- (atan . atan_dfa)
- (sinh . sinh_dfa) (cosh . cosh_dfa) (tanh . tanh_dfa)
- (asinh . asinh_dfa) (acosh . acosh_dfa)
- (exp . exp_dfa) (log . log_dfa) (sqrt . sqrt_dfat) (conjugate . conjugate_dfa)
- (realpart . realpart_dfa) (imagpart . imagpart_dfa) (abs . abs_dfa)))
+ ;; List of functions that should map real to real. The list
+ ;; is not that long, since many functions, such as sqrt, have
+ ;; potentially complex output
+ '((sin . sin_dfa-to-dfa) (cos . cos_dfa-to-dfa) (tan . tan_dfa-to-dfa)
+ (atan . atan_dfa-to-dfa)
+ (sinh . sinh_dfa-to-dfa) (cosh . cosh_dfa-to-dfa) (tanh . tanh_dfa-to-dfa)
+ (asinh . asinh_dfa-to-dfa) (acosh . acosh_dfa-to-dfa)
+ (exp . exp_dfa-to-dfa)
+ (abs . abs_dfa-to-dfa)))
-(defmacro defun-dfa-to-dfa (name op)
+;; the iterative version
+#+lisplab-iterative (defmacro defun-dfa-to-dfa-iterative (name op)
(let ((a (gensym))
(out (gensym))
(i (gensym)))
@@ -53,18 +63,44 @@
(setf (aref ,out ,i) (,op (aref ,a ,i))))
,out)))
+;; the iterative version
+#+lisplab-iterative (defmacro expand-ordinary-functions-dfa-to-dfa-map-iterative ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-dfa-to-dfa-iterative ,(cdr x) ,(car x)))
+ +ordinary-functions-dfa-to-dfa-map+)))
+
+;; the iterative version
+#+lisplab-iterative (expand-ordinary-functions-dfa-to-dfa-map-iterative)
+
+;;; Non-iterative version
+#-lisplab-iterative
+(defmacro defun-dfa-to-dfa (name op)
+ (let ((a (gensym))
+ (out (gensym)))
+ `(defun ,name (,a ,out)
+ (declare (type type-blas-store ,a ,out))
+ (map-into ,out #',op ,a)
+ ,out)))
+
+;;; Non-iterative version
+#-lisplab-iterative
(defmacro expand-ordinary-functions-dfa-to-dfa-map ()
(cons 'progn
(mapcar (lambda (x)
`(defun-dfa-to-dfa ,(cdr x) ,(car x)))
+ordinary-functions-dfa-to-dfa-map+)))
-(expand-ordinary-functions-dfa-to-dfa-map)
+;;; Non-iterative version
+#-lisplab-iterative
+(expand-ordinary-functions-dfa-to-dfa-map) ; the iterative version
+
;;; double float to complex double float
(define-constant +ordinary-functions-dfa-to-cdfa-map+
- '((asin . asin_dfa) (acos . acos_dfa) (atanh . atanh_dfa)))
+ '((asin . asin_dfa) (acos . acos_dfa) (atanh . atanh_dfa)
+ (log . log_dfa) (sqrt . sqrt_dfa)))
(defmacro defun-dfa-to-cdfa (name op)
(let ((a (gensym))
@@ -73,7 +109,8 @@
`(defun ,name (,a ,out)
(declare (type type-blas-store ,a ,out))
(dotimes (,i (length ,a))
- (setf (vref-blas-complex-store ,out ,i) (,op (aref ,a ,i))))
+ (setf (vref-blas-complex-store ,out ,i)
+ (coerce (,op (aref ,a ,i)) '(complex double-float))))
,out)))
(defmacro expand-ordinary-functions-dfa-to-cdfa-map ()
@@ -86,26 +123,26 @@
;;; Complex double float to double float
-(define-constant +ordinary-functions-cdfa-to-dfa-map+
- '((realpart . realpart_cdfa) (imagpart . imagpart_cdfa) (abs . abs_cdfa)))
-
-(defmacro defun-cdfa-to-dfa (name op)
- (let ((a (gensym))
- (out (gensym))
- (i (gensym)))
- `(defun ,name (,a ,out)
- (declare (type type-blas-store ,a ,out))
- (dotimes (,i (floor (length ,a) 2))
- (setf (aref ,out ,i) (,op (vref-blas-complex-store ,a ,i))))
- ,out)))
-
-(defmacro expand-ordinary-functions-cdfa-to-dfa-map ()
- (cons 'progn
- (mapcar (lambda (x)
- `(defun-cdfa-to-dfa ,(cdr x) ,(car x)))
- +ordinary-functions-cdfa-to-dfa-map+)))
-
-(expand-ordinary-functions-cdfa-to-dfa-map)
+(defun abs_cdfa (a out)
+ (declare (type type-blas-store a out))
+ (dotimes (i (length out))
+ (setf (aref out i)
+ (abs (vref-blas-complex-store a i))))
+ out)
+
+(defun realpart_cdfa (a out)
+ (declare (type type-blas-store a out))
+ (dotimes (i (length out))
+ (setf (aref out i)
+ (aref a (truly-the type-blas-idx (* 2 i)))))
+ out)
+
+(defun imagpart_cdfa (a out)
+ (declare (type type-blas-store a out))
+ (dotimes (i (length out))
+ (setf (aref out i)
+ (aref a (truly-the type-blas-idx (1+ (* 2 i))))))
+ out)
;;; Complex double float to complex double float
@@ -114,7 +151,8 @@
(atan . atan_cdfa)
(sinh . sinh_cdfa) (cosh . cosh_cdfa) (tanh . tanh_cdfa)
(asinh . asinh_cdfa) (acosh . acosh_cdfa)
- (exp . exp_cdfa) (log . log_cdfa) (sqrt . sqrt_cdfat) (conjugate . conjugate_cdfa)
+ (exp . exp_cdfa) (log . log_cdfa) (sqrt . sqrt_cdfat)
+ #+nil (conjugate . conjugate_cdfa) ;; separate implementation!
(asin . asin_cdfa) (acos . acos_cdfa) (atanh . atanh_cdfa)))
(defmacro defun-cdfa-to-cdfa (name op)
@@ -133,4 +171,16 @@
`(defun-cdfa-to-cdfa ,(cdr x) ,(car x)))
+ordinary-functions-cdfa-to-cdfa-map+)))
-(expand-ordinary-functions-cdfa-to-cdfa-map)
\ No newline at end of file
+(expand-ordinary-functions-cdfa-to-cdfa-map)
+
+;; Conjugate
+
+(defun conjugate_cdfa (a out)
+ (declare (type type-blas-store a out))
+ (dotimes (i (floor (length a) 2))
+ (let* ((2i (* 2 i))
+ (2i+1 (1+ 2i)))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref out 2i) (aref a 2i)
+ (aref out 2i+1) (- (aref a 2i+1)))))
+ out)
\ No newline at end of file
More information about the lisplab-cvs
mailing list