[lisplab-cvs] r107 - src/matrix
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Mon Nov 2 20:06:51 UTC 2009
Author: jivestgarden
Date: Mon Nov 2 15:06:51 2009
New Revision: 107
Log:
In the middle of refactoring. May not work
Added:
src/matrix/store-operators.lisp
src/matrix/store-ordinary-functions.lisp
Modified:
lisplab.asd
src/matrix/level1-util.lisp
src/matrix/level2-matrix-dge.lisp
src/matrix/level2-matrix-zge.lisp
Modified: lisplab.asd
==============================================================================
--- lisplab.asd (original)
+++ lisplab.asd Mon Nov 2 15:06:51 2009
@@ -62,7 +62,12 @@
:components
(
(:file "level1-interface")
- (:file "level1-util")
+
+ ;; These three should be independent of the rest
+ (:file "level1-util")
+ (:file "store-operators")
+ (:file "store-ordinary-functions")
+
(:file "level1-classes")
(:file "level1-constructors")
(:file "level1-matrix")
Modified: src/matrix/level1-util.lisp
==============================================================================
--- src/matrix/level1-util.lisp (original)
+++ src/matrix/level1-util.lisp Mon Nov 2 15:06:51 2009
@@ -19,9 +19,19 @@
;;; 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
+;;; many classes such as matrix-base-dge and matrix-base-ddi, etc.
+;;;
+;;; The content of this file must be highly optimized
+;;; and should not depend anything exept Common Lisp itself.
(in-package :lisplab)
+;;; Things that are common both for real and complex stores
+
(deftype type-blas-store ()
'(simple-array double-float (*)))
@@ -32,10 +42,7 @@
#-:sbcl (deftype type-blas-idx ()
'fixnum)
-
(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
@@ -44,6 +51,21 @@
type-blas-idx)
column-major-idx))
+(defun column-major-idx (i j rows)
+ (truly-the type-blas-idx (+ i (truly-the type-blas-idx (* j rows)))))
+
+(defun copy-matrix-stores (a b)
+ (let ((len (length a)))
+ (declare (type type-blas-store a b)
+ (type type-blas-idx len))
+ (dotimes (i len)
+ (setf (aref b i) (aref a i))))
+ b)
+
+;;;; The real store
+
+(declaim (inline ref-blas-real-store (setf ref-blas-real-store)))
+
(declaim (ftype (function
(type-blas-store
type-blas-idx
@@ -52,6 +74,14 @@
double-float)
ref-blas-real-store))
+(defun ref-blas-real-store (store row col rows)
+ "Matrix accessor for the real blas store"
+ (aref (truly-the type-blas-store store)
+ (truly-the type-blas-idx
+ (column-major-idx (truly-the type-blas-idx row)
+ (truly-the type-blas-idx col)
+ rows))))
+
(declaim (ftype (function
(double-float
type-blas-store
@@ -61,35 +91,6 @@
)
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 (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"
- (aref (truly-the type-blas-store store)
- (truly-the type-blas-idx
- (column-major-idx (truly-the type-blas-idx row)
- (truly-the type-blas-idx col)
- rows))))
-
(defun (setf ref-blas-real-store) (value store row col rows)
(setf (aref (truly-the type-blas-store store)
@@ -101,11 +102,13 @@
value)
(defun allocate-real-store (size &optional (initial-element 0.0))
+ ;; All matrix double and complex double constructors
+ ;; should call this one
(let ((x (coerce initial-element 'double-float)))
(declare (type double-float x)
(type type-blas-idx size))
- ;; Stupid efficiency hack, on SBCL. All matrix double and complex double
- ;; should call this one
+ ;; Stupid efficiency hack for SBCL. Allocations of arrays with zeros
+ ;; is significantly faster than others!
(if (= x 0.0)
(make-array size
:element-type 'double-float
@@ -113,9 +116,55 @@
(make-array size
:element-type 'double-float
:initial-element x))))
-
+
+;;; Hopfully helpful iterators
+
+(defmacro with-indices-df-1 (a m idx &body body)
+ "Iterats over all the indices of one array"
+ `(let ((,m ,a))
+ (declare (type type-blas-store ,m))
+ (dotimes (,idx (length ,m))
+ (declare (type type-blas-idx ,idx))
+ , at body)))
+
+(defmacro with-elements-df-1 (a elm &body body)
+ "Iterats over all the elements of one array"
+ (let ((m (gensym))
+ (idx (gensym)))
+ `(let ((,m ,a))
+ (declare (type type-blas-store ,m))
+ (dotimes (,idx (length ,m))
+ (declare (type type-blas-idx ,idx))
+ (let ((,elm (aref ,m ,idx)))
+ (declare (type double-float ,elm))
+ , at body)))))
+
+;;;; The complex store
+
+(defun allocate-complex-store (size &optional (value 0.0))
+ (let* ((2size (* 2 size))
+ (rv (coerce (realpart value) 'double-float))
+ (iv (coerce (imagpart value) 'double-float))
+ (store (allocate-real-store 2size iv)))
+ (declare (type type-blas-idx 2size)
+ (type double-float rv iv))
+ (when (/= rv iv)
+ (loop for i from 0 below 2size by 2 do
+ (setf (aref store i) rv)))
+ store))
+
+(declaim (inline ref-blas-complex-store (setf ref-blas-complex-store)))
+
+(declaim (ftype (function
+ (type-blas-store
+ type-blas-idx
+ type-blas-idx
+ type-blas-idx)
+ (complex double-float))
+ ref-blas-complex-store))
+
(defun ref-blas-complex-store (store row col rows)
- "Accessor for the complet blas store"
+ "Matrix accessor for the complet blas store"
(let ((idx (truly-the type-blas-idx
(* 2 (column-major-idx (truly-the type-blas-idx row)
(truly-the type-blas-idx col)
@@ -124,6 +173,16 @@
(complex (aref store idx)
(aref store (1+ idx)))))
+(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 (setf ref-blas-complex-store) (value store row col rows)
(let ((idx (truly-the type-blas-idx
(* 2 (column-major-idx (truly-the type-blas-idx row)
@@ -134,22 +193,31 @@
(aref store (1+ idx)) (imagpart value))
value))
-(defun allocate-complex-store (size &optional (value 0.0))
- (let* ((2size (* 2 size))
- (rv (coerce (realpart value) 'double-float))
- (iv (coerce (imagpart value) 'double-float))
- (store (allocate-real-store 2size iv)))
- (declare (type type-blas-idx 2size)
- (type double-float rv iv))
- (when (/= rv iv)
- (loop for i from 0 below 2size by 2 do
- (setf (aref store i) rv)))
- store))
+(declaim (ftype (function
+ (type-blas-store
+ type-blas-idx)
+ (complex double-float))
+ vref-blas-complex-store))
+
+(defun vref-blas-complex-store (store idx)
+ "Matrix accessor for the complex blas store"
+ (let ((idx2 (truly-the type-blas-idx (* 2 idx))))
+ (declare (type-blas-idx idx2))
+ (complex (aref store idx2)
+ (aref store (1+ idx2)))))
+
+(declaim (ftype (function
+ ((complex double-float)
+ type-blas-store
+ type-blas-idx
+ )
+ (complex double-float))
+ (setf vref-blas-complex-store)))
+
+(defun (setf vref-blas-complex-store) (value store idx)
+ (let ((idx2 (truly-the type-blas-idx (* 2 idx))))
+ (declare (type-blas-idx idx2))
+ (setf (aref store idx2) (realpart value)
+ (aref store (1+ idx2)) (imagpart value))
+ value))
-(defun copy-matrix-stores (a b)
- (let ((len (length a)))
- (declare (type type-blas-store a b)
- (type type-blas-idx len))
- (dotimes (i len)
- (setf (aref b i) (aref a i))))
- b)
\ 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 Mon Nov 2 15:06:51 2009
@@ -48,12 +48,9 @@
out)
(defmethod msum ((m matrix-base-dge))
- (let ((sum 0.0)
- (m0 (matrix-store m)))
- (declare (type double-float sum)
- (type type-blas-store m0))
- (dotimes (i (length m0))
- (incf sum (aref m0 i)))
+ (let ((sum 0.0))
+ (with-elements-df-1 (matrix-store m) x
+ (incf sum x))
sum))
(defmethod .some (pred (a matrix-base-dge) &rest args)
@@ -64,62 +61,128 @@
(let ((stores (mapcar #'matrix-store (cons a args))))
(apply #'every pred stores)))
-(defmacro def-binary-op-matrix-base-dge (new old)
+
+;;; Matrix and real
+
+(define-constant +generic-function-dfa-df-map+
+ '((.add . +_dfa-df)
+ (.sub . -_dfa-df)
+ (.mul . *_dfa-df)
+ (.div . /_dfa-df)
+ (.expt . ^_dfa-df)
+ (.max . max_dfa-df)
+ (.min . min_dfa-df)))
+
+(defmacro defmethod-dfa-df (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)))))
+ (c (gensym "c")))
+ `(defmethod ,name ((,a matrix-base-dge) (,b real))
+ (let ((,c (mcreate ,a)))
+ (,underlying-function (matrix-store ,a) (coerce ,b 'double-float) (matrix-store ,c))
+ ,c))))
-(def-binary-op-matrix-base-dge .add +)
+(defmacro expand-generic-function-dfa-df-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-dfa-df ,(car x) ,(cdr x)))
+ +generic-function-dfa-df-map+)))
+
+(expand-generic-function-dfa-df-map)
+
+;;; Real and matrix
+
+(define-constant +generic-function-df-dfa-map+
+ '((.add . +_df-dfa)
+ (.sub . -_df-dfa)
+ (.mul . *_df-dfa)
+ (.div . /_df-dfa)
+ (.expt . ^_df-dfa)
+ (.max . max_df-dfa)
+ (.min . min_df-dfa)))
-(def-binary-op-matrix-base-dge .sub -)
+(defmacro defmethod-df-dfa (name underlying-function)
+ (let ((a (gensym "a"))
+ (b (gensym "b"))
+ (c (gensym "c")))
+ `(defmethod ,name ((,a real) (,b matrix-base-dge))
+ (let ((,c (mcreate ,b)))
+ (,underlying-function (coerce ,a 'double-float) (matrix-store ,b) (matrix-store ,c))
+ ,c))))
-(def-binary-op-matrix-base-dge .mul *)
+(defmacro expand-generic-function-df-dfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-df-dfa ,(car x) ,(cdr x)))
+ +generic-function-df-dfa-map+)))
+
+(expand-generic-function-df-dfa-map)
+
+;;;; Matrix and matrix
+
+(define-constant +generic-function-dfa-dfa-map+
+ '((.add . +_dfa-dfa)
+ (.sub . -_dfa-dfa)
+ (.mul . *_dfa-dfa)
+ (.div . /_dfa-dfa)
+ (.expt . ^_dfa-dfa)
+ (.max . max_dfa-dfa)
+ (.min . min_dfa-dfa)))
-(def-binary-op-matrix-base-dge .div /)
+(defmacro defmethod-dfa-dfa (name underlying-function)
+ (let ((a (gensym "a"))
+ (b (gensym "b"))
+ (c (gensym "c")))
+ `(defmethod ,name ((,a matrix-base-dge) (,b matrix-base-dge))
+ (let ((,c (mcreate ,a)))
+ (,underlying-function (matrix-store ,a) (matrix-store ,b) (matrix-store ,c))
+ ,c))))
-(def-binary-op-matrix-base-dge .expt expt)
+(defmacro expand-generic-function-dfa-dfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-dfa-dfa ,(car x) ,(cdr x)))
+ +generic-function-dfa-dfa-map+)))
-(def-binary-op-matrix-base-dge .min min)
+(expand-generic-function-dfa-dfa-map)
-(def-binary-op-matrix-base-dge .max max)
+;;; The ordinary functions
+
+;;;; Matrix and matrix
+
+(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) (.sqrt . sqrt_dfat) (.conjugate . conjugate_dfa)
+ (.realpart . realpart_dfa) (.imagpart . imagpart_dfa) (.abs . abs_dfa)))
+
+(defmacro defmethod-dfa-to-dfa (name underlying-function)
+ (let ((a (gensym "a"))
+ (b (gensym "b")))
+ `(defmethod ,name ((,a matrix-base-dge))
+ (let ((,b (mcreate ,a)))
+ (,underlying-function (matrix-store ,a) (matrix-store ,b) )
+ ,b))))
+
+(defmacro expand-generic-function-dfa-to-dfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-dfa-to-dfa ,(car x) ,(cdr x)))
+ +generic-function-dfa-to-dfa-map+)))
+(expand-generic-function-dfa-to-dfa-map)
+
+;;; The reset 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
@@ -196,27 +259,66 @@
(defmethod .besh2 (n (x matrix-base-dge))
(each-matrix-element-df-to-complex-df x (.besh2 n x)))
+|#
#|
+;;; Old code
-(defmethod .imagpart ((a matrix-base-dge))
- (mcreate a 0))
+(defmacro def-binary-op-matrix-base-dge (new old)
+ (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)))))
-(defmethod .realpart ((a matrix-base-dge))
- (copy a))
+(def-binary-op-matrix-base-dge .add +)
-(defmethod .abs ((a matrix-base-dge))
- (let ((b (mcreate a)))
- (copy-contents a b #'abs)
- b))
+(def-binary-op-matrix-base-dge .sub -)
+(def-binary-op-matrix-base-dge .mul *)
-(defmacro expand-on-matrix-dge-lisplab-two-argument-functions-alist ()
- (cons 'progn
- (mapcar (lambda (name)
- `(def-binary-op-matrix-base-dge ,(car name) ,(cdr name)))
- +lisplab-two-argument-functions-alist+)))
+(def-binary-op-matrix-base-dge .div /)
+
+(def-binary-op-matrix-base-dge .expt expt)
+
+(def-binary-op-matrix-base-dge .min min)
+
+(def-binary-op-matrix-base-dge .max max)
-(expand-on-matrix-dge-lisplab-two-argument-functions-alist)
|#
\ No newline at end of file
Modified: src/matrix/level2-matrix-zge.lisp
==============================================================================
--- src/matrix/level2-matrix-zge.lisp (original)
+++ src/matrix/level2-matrix-zge.lisp Mon Nov 2 15:06:51 2009
@@ -107,6 +107,8 @@
(sqrt (+ (* x x) (* y y))))))
b))
+;;; Old code
+
(defmacro def-binary-op-matrix-base-zge (new old)
;;; TODO speed up for real numbers. Is it worth the work?
(let ((a (gensym "a"))
@@ -191,6 +193,97 @@
(def-binary-op-matrix-base-zge .expt expt)
+;;;; new code
+
+ ;;; Matrix and complex
+
+(define-constant +generic-function-cdfa-cdf-map+
+ '((.add . +_cdfa-cdf)
+ (.sub . -_cdfa-cdf)
+ (.mul . *_cdfa-cdf)
+ (.div . /_cdfa-cdf)
+ (.expt . ^_cdfa-cdf)
+ (.max . max_cdfa-cdf)
+ (.min . min_cdfa-cdf)))
+
+(defmacro defmethod-cdfa-cdf (name underlying-function)
+ (let ((a (gensym "a"))
+ (b (gensym "b"))
+ (c (gensym "c")))
+ `(defmethod ,name ((,a matrix-base-zge) (,b number))
+ (let ((,c (mcreate ,a)))
+ (,underlying-function (matrix-store ,a)
+ (coerce ,b '(complex double-float))
+ (matrix-store ,c))
+ ,c))))
+
+(defmacro expand-generic-function-cdfa-cdf-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-cdfa-cdf ,(car x) ,(cdr x)))
+ +generic-function-cdfa-cdf-map+)))
+
+(expand-generic-function-cdfa-cdf-map)
+
+;;; Real and matrix
+
+(define-constant +generic-function-cdf-cdfa-map+
+ '((.add . +_cdf-cdfa)
+ (.sub . -_cdf-cdfa)
+ (.mul . *_cdf-cdfa)
+ (.div . /_cdf-cdfa)
+ (.expt . ^_cdf-cdfa)
+ (.max . max_cdf-cdfa)
+ (.min . min_cdf-cdfa)))
+
+(defmacro defmethod-cdf-cdfa (name underlying-function)
+ (let ((a (gensym "a"))
+ (b (gensym "b"))
+ (c (gensym "c")))
+ `(defmethod ,name ((,a number) (,b matrix-base-zge))
+ (let ((,c (mcreate ,b)))
+ (,underlying-function (coerce ,a '(complex double-float))
+ (matrix-store ,b)
+ (matrix-store ,c))
+ ,c))))
+
+(defmacro expand-generic-function-cdf-cdfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-cdf-cdfa ,(car x) ,(cdr x)))
+ +generic-function-cdf-cdfa-map+)))
+
+(expand-generic-function-cdf-cdfa-map)
+
+;;;; Matrix and matrix
+
+(define-constant +generic-function-cdfa-cdfa-map+
+ '((.add . +_cdfa-cdfa)
+ (.sub . -_cdfa-cdfa)
+ (.mul . *_cdfa-cdfa)
+ (.div . /_cdfa-cdfa)
+ (.expt . ^_cdfa-cdfa)
+ (.max . max_cdfa-cdfa)
+ (.min . min_cdfa-cdfa)))
+
+(defmacro defmethod-cdfa-cdfa (name underlying-function)
+ (let ((a (gensym "a"))
+ (b (gensym "b"))
+ (c (gensym "c")))
+ `(defmethod ,name ((,a matrix-base-zge) (,b matrix-base-zge))
+ (let ((,c (mcreate ,a)))
+ (,underlying-function (matrix-store ,a) (matrix-store ,b) (matrix-store ,c))
+ ,c))))
+
+(defmacro expand-generic-function-cdfa-cdfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defmethod-cdfa-cdfa ,(car x) ,(cdr x)))
+ +generic-function-cdfa-cdfa-map+)))
+
+(expand-generic-function-cdfa-cdfa-map)
+
+
(defmacro each-element-function-matrix-base-zge (x form)
"Applies a form on each element of an matrix-base-zge."
Added: src/matrix/store-operators.lisp
==============================================================================
--- (empty file)
+++ src/matrix/store-operators.lisp Mon Nov 2 15:06:51 2009
@@ -0,0 +1,249 @@
+;;; Lisplab, store-operators.lisp
+;;; Double float and complex double float operators (such as +,-,*, etc) on
+;;; simple arrays. Used by the matrix implementations.
+;;;
+
+;;; 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.
+
+;;; 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
+;;; many classes such as matrix-base-dge and matrix-base-ddi, etc.
+;;;
+;;; The content of this file must be highly optimized
+;;; and should not depend anything exept Common Lisp itself.
+
+(in-package :lisplab)
+
+;;; 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 reason for generating ordinary functions and not using methods,
+;;; is that the real and complex stores have the same type!
+
+;;; 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
+;;; be taken. This is also the reason why it's hard to generate more code
+;;; automatically.
+
+;;; The below code generates ordinary lisp functions
+;;; for element-wise operations on simple double-float arrays.
+;;; 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)
+
+
+;;; Array and number
+
+(define-constant +operators-dfa-df-map+
+ '((+ . +_dfa-df)
+ (- . -_dfa-df)
+ (* . *_dfa-df)
+ (/ . /_dfa-df)
+ (expt . ^_dfa-df)
+ (max . max_dfa-df)
+ (min . min_dfa-df)
+ (log . log_dfa-df) ;; Note the log here
+ ))
+
+(defmacro defun-dfa-df (name op)
+ (let ((a (gensym))
+ (b (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(defun ,name (,a ,b ,out)
+ (declare (type type-blas-store ,a ,out)
+ (type double-float ,b))
+ (dotimes (,i (length ,a))
+ (setf (aref ,out ,i) (,op (aref ,a ,i) ,b)))
+ ,out)))
+
+(defmacro expand-operators-dfa-df-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-dfa-df ,(cdr x) ,(car x)))
+ +operators-dfa-df-map+)))
+
+(expand-operators-dfa-df-map)
+
+;;; The three parts of code below contains some common strucutre that could
+;;; in principle be joined, and there is also some clumsy code,
+;;; but I have not the energy to find out how to clean up.
+
+;;; Number and array
+
+(define-constant +operators-df-dfa-map+
+ '((+ . +_df-dfa)
+ (- . -_df-dfa)
+ (* . *_df-dfa)
+ (/ . /_df-dfa)
+ (expt . ^_df-dfa)
+ (max . max_df-dfa)
+ (min . min_df-dfa)
+ (log . log_df-dfa)))
+
+(defmacro defun-df-dfa (name op)
+ (let ((a (gensym))
+ (b (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(defun ,name (,a ,b ,out)
+ (declare (type type-blas-store ,b ,out)
+ (type double-float ,a))
+ (dotimes (,i (length ,b))
+ (setf (aref ,out ,i) (,op ,a (aref ,b ,i))))
+ ,out)))
+
+(defmacro expand-operators-df-dfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-df-dfa ,(cdr x) ,(car x)))
+ +operators-df-dfa-map+)))
+
+(expand-operators-df-dfa-map)
+
+;;; Array and array
+
+(define-constant +operators-dfa-dfa-map+
+ '((+ . +_dfa-dfa)
+ (- . -_dfa-dfa)
+ (* . *_dfa-dfa)
+ (/ . /_dfa-dfa)
+ (expt . ^_dfa-dfa)
+ (max . max_dfa-dfa)
+ (min . min_dfa-dfa)
+ (log . log_dfa-dfa)))
+
+(defmacro defun-dfa-dfa (name op)
+ (let ((a (gensym))
+ (b (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(defun ,name (,a ,b ,out)
+ (declare (type type-blas-store ,a ,b ,out))
+ (dotimes (,i (length ,a))
+ (setf (aref ,out ,i) (,op (aref ,a ,i) (aref ,b ,i))))
+ ,out)))
+
+(defmacro expand-operators-dfa-dfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-dfa-dfa ,(cdr x) ,(car x)))
+ +operators-dfa-dfa-map+)))
+
+(expand-operators-dfa-dfa-map)
+
+
+;;; Now the complex operators
+
+
+;;; Array and number
+
+(define-constant +operators-cdfa-cdf-map+
+ '((+ . +_cdfa-cdf)
+ (- . -_cdfa-cdf)
+ (* . *_cdfa-cdf)
+ (/ . /_cdfa-cdf)
+ (expt . ^_cdfa-cdf)
+ ))
+
+(defmacro defun-cdfa-cdf (name op)
+ (let ((a (gensym))
+ (b (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(defun ,name (,a ,b ,out)
+ (declare (type type-blas-store ,a ,out)
+ (type (complex double-float) ,b))
+ (dotimes (,i (floor (length ,a) 2))
+ (setf (vref-blas-complex-store ,out ,i)
+ (coerce (,op (vref-blas-complex-store ,a ,i) ,b) '(complex double-float))))
+ ,out)))
+
+(defmacro expand-operators-cdfa-cdf-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-cdfa-cdf ,(cdr x) ,(car x)))
+ +operators-cdfa-cdf-map+)))
+
+(expand-operators-cdfa-cdf-map)
+
+;;; The three parts of code below contains some common strucutre that could
+;;; in principle be joined, and there is also some clumsy code,
+;;; but I have not the energy to find out how to clean up.
+
+;;; Number and array
+
+(define-constant +operators-cdf-cdfa-map+
+ '((+ . +_cdf-cdfa)
+ (- . -_cdf-cdfa)
+ (* . *_cdf-cdfa)
+ (/ . /_cdf-cdfa)
+ (expt . ^_cdf-cdfa)))
+
+(defmacro defun-cdf-cdfa (name op)
+ (let ((a (gensym))
+ (b (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(defun ,name (,a ,b ,out)
+ (declare (type type-blas-store ,b ,out)
+ (type (complex double-float) ,a))
+ (dotimes (,i (floor (length ,b) 2))
+ (setf (vref-blas-complex-store ,out ,i) (,op ,a (vref-blas-complex-store ,b ,i))))
+ ,out)))
+
+(defmacro expand-operators-cdf-cdfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-cdf-cdfa ,(cdr x) ,(car x)))
+ +operators-cdf-cdfa-map+)))
+
+(expand-operators-cdf-cdfa-map)
+
+;;; Array and array
+
+(define-constant +operators-cdfa-cdfa-map+
+ '((+ . +_cdfa-cdfa)
+ (- . -_cdfa-cdfa)
+ (* . *_cdfa-cdfa)
+ (/ . /_cdfa-cdfa)
+ (expt . ^_cdfa-cdfa)))
+
+(defmacro defun-cdfa-cdfa (name op)
+ (let ((a (gensym))
+ (b (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(defun ,name (,a ,b ,out)
+ (declare (type type-blas-store ,a ,b ,out))
+ (dotimes (,i (floor (length ,a) 2))
+ (setf (vref-blas-complex-store ,out ,i)
+ (,op (vref-blas-complex-store ,a ,i) (vref-blas-complex-store ,b ,i))))
+ ,out)))
+
+(defmacro expand-operators-cdfa-cdfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-cdfa-cdfa ,(cdr x) ,(car x)))
+ +operators-cdfa-cdfa-map+)))
+
+(expand-operators-cdfa-cdfa-map)
\ No newline at end of file
Added: src/matrix/store-ordinary-functions.lisp
==============================================================================
--- (empty file)
+++ src/matrix/store-ordinary-functions.lisp Mon Nov 2 15:06:51 2009
@@ -0,0 +1,136 @@
+;;; Lisplab, store-ordinary-functions.lisp
+;;; Double float and complex double float ordinary functions (such as sin, log, etc) on
+;;; simple arrays. Used by the matrix implementations.
+;;;
+
+;;; 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.
+
+;;; 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
+;;; many classes such as matrix-base-dge and matrix-base-ddi, etc.
+;;;
+;;; The content of this file must be highly optimized
+;;; and should not depend anything exept Common Lisp itself.
+
+(in-package :lisplab)
+
+;;; Now the ordinary functions
+
+;;; 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) (sqrt . sqrt_dfat) (conjugate . conjugate_dfa)
+ (realpart . realpart_dfa) (imagpart . imagpart_dfa) (abs . abs_dfa)))
+
+(defmacro defun-dfa-to-dfa (name op)
+ (let ((a (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(defun ,name (,a ,out)
+ (declare (type type-blas-store ,a ,out))
+ (dotimes (,i (length ,a))
+ (setf (aref ,out ,i) (,op (aref ,a ,i))))
+ ,out)))
+
+(defmacro expand-ordinary-functions-dfa-to-dfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-dfa ,(cdr x) ,(car x)))
+ +ordinary-functions-dfa-map+)))
+
+(expand-ordinary-functions-dfa-to-dfa-map)
+
+;;; double float to complex double float
+
+(define-constant +ordinary-functions-dfa-to-cdfa-map+
+ '((asin . asin_dfa) (acos . acos_dfa) (atanh . atanh_dfa)))
+
+(defmacro defun-dfa-to-cdfa (name op)
+ (let ((a (gensym))
+ (out (gensym))
+ (i (gensym)))
+ `(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))))
+ ,out)))
+
+(defmacro expand-ordinary-functions-dfa-to-cdfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-dfa ,(cdr x) ,(car x)))
+ +ordinary-functions-dfa-to-cdfa-map+)))
+
+(expand-ordinary-functions-dfa-to-cdfa-map)
+
+;;; Complex double float to double float
+
+(define-constant +ordinary-functions-cdfa-to-dfa-map+
+ '((realpart . realpart_cdfa) (imagpart . imagpart_cdfa) (abs . cabs_dfa)))
+
+(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-dfa ,(cdr x) ,(car x)))
+ +ordinary-functions-cdfa-to-dfa-map+)))
+
+(expand-ordinary-functions-cdfa-to-dfa-map)
+
+;;; Complex double float to complex double float
+
+(define-constant +ordinary-functions-cdfa-to-cdfa-map+
+ '((sin . sin_cdfa) (cos . cos_cdfa) (tan . tan_cdfa)
+ (atan . atan_cdfa)
+ (sinh . sinh_cdfa) (cosh . cosh_cdfa) (tanh . tanh_cdfa)
+ (asinh . asinh_cdfa) (acosh . acosh_cdfa)
+ (exp . exp_cdfa) (sqrt . sqrt_cdfat) (conjugate . conjugate_cdfa)
+ (asin . asin_cdfa) (acos . acos_cdfa) (atanh . atanh_cdfa)))
+
+(defmacro defun-cdfa-to-cdfa (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 (vref-blas-complex-store ,out ,i) (,op (vref-blas-complex-store ,a ,i))))
+ ,out)))
+
+(defmacro expand-ordinary-functions-cdfa-to-cdfa-map ()
+ (cons 'progn
+ (mapcar (lambda (x)
+ `(defun-dfa ,(cdr x) ,(car x)))
+ +ordinary-functions-cdfa-to-cdfa-map+)))
+
+(expand-ordinary-functions-cdfa-to-cdfa-map)
\ No newline at end of file
More information about the lisplab-cvs
mailing list