Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-127-gcbb7c2b
|
From: Akshay S. <ak...@us...> - 2012-07-12 10:48:23
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, tensor has been updated
via cbb7c2bfaa2dedc65e56be04c1469e46be975801 (commit)
from 63d6b10a662cb7b8ad0b3dfd288db7a5f921abff (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit cbb7c2bfaa2dedc65e56be04c1469e46be975801
Author: Akshay Srinivasan <aks...@gm...>
Date: Thu Jul 12 16:13:06 2012 +0530
o Added a fortran-call-lower-bound to axpy and copy.
o Added a new num-axpy! function generating macro, added methods to axpy!
diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp
index edef2d3..2de1898 100644
--- a/src/level-1/axpy.lisp
+++ b/src/level-1/axpy.lisp
@@ -28,7 +28,7 @@
(in-package #:matlisp)
-(defmacro generate-typed-axpy! (func (tensor-class blas-func))
+(defmacro generate-typed-axpy! (func (tensor-class blas-func fortran-lb))
;;Be very careful when using functions generated by this macro.
;;Indexes can be tricky and this has no safety net
;;Use only after checking the arguments for compatibility.
@@ -37,29 +37,88 @@
`(defun ,func (alpha from to)
(declare (type ,tensor-class from to)
(type ,(getf opt :element-type) alpha))
- (if-let (strd-p (blas-copyable-p from to))
- (,blas-func (number-of-elements from) alpha (store from) (first strd-p) (store to) (second strd-p) (head from) (head to))
- (let ((f-sto (store from))
- (t-sto (store to)))
- (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto))
- (very-quickly
- ;;One would question the wisdom in calling the Fortran method here.
- ;;Simple benchmarks proved that SBCL is as quick as or better than
- ;;OpenBLAS's methods
- (mod-dotimes (idx (dimensions from))
- with (linear-sums
- (f-of (strides from) (head from))
- (t-of (strides to) (head to)))
- do (let ((f-val ,(funcall (getf opt :reader) 'f-sto 'f-of))
- (t-val ,(funcall (getf opt :reader) 't-sto 't-of)))
- (declare (type ,(getf opt :element-type) f-val t-val))
- (let ((t-new (+ (* f-val alpha) t-val)))
- (declare (type ,(getf opt :element-type) t-new))
- ,(funcall (getf opt :value-writer) 't-new 't-sto 't-of)))))))
+ (let ((strd-p (blas-copyable-p from to))
+ (call-fortran? (> (number-of-elements to)
+ ,fortran-lb)))
+ (cond
+ ((and strd-p call-fortran?)
+ (,blas-func (number-of-elements from) alpha
+ (store from) (first strd-p)
+ (store to) (second strd-p)
+ (head from) (head to)))
+ (t
+ (let ((f-sto (store from))
+ (t-sto (store to)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto))
+ (very-quickly
+ ;;One would question the wisdom in calling the Fortran method here.
+ ;;Simple benchmarks proved that SBCL is as quick as or better than
+ ;;OpenBLAS's methods
+ (mod-dotimes (idx (dimensions from))
+ with (linear-sums
+ (f-of (strides from) (head from))
+ (t-of (strides to) (head to)))
+ do (let ((f-val ,(funcall (getf opt :reader) 'f-sto 'f-of))
+ (t-val ,(funcall (getf opt :reader) 't-sto 't-of)))
+ (declare (type ,(getf opt :element-type) f-val t-val))
+ (let ((t-new (+ (* f-val alpha) t-val)))
+ (declare (type ,(getf opt :element-type) t-new))
+ ,(funcall (getf opt :value-writer) 't-new 't-sto 't-of)))))))))
to)))
-(generate-typed-axpy! real-typed-axpy! (real-tensor daxpy))
-(generate-typed-axpy! complex-typed-axpy! (complex-tensor zaxpy))
+(defmacro generate-typed-num-axpy! (func (tensor-class blas-func fortran-lb))
+ ;;Be very careful when using functions generated by this macro.
+ ;;Indexes can be tricky and this has no safety net
+ ;;(you don't see a matrix-ref do you ?)
+ ;;Use only after checking the arguments for compatibility.
+ (let* ((opt (get-tensor-class-optimization tensor-class)))
+ (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
+ `(defun ,func (num-from to)
+ (declare (type ,tensor-class to)
+ (type ,(getf opt :element-type) num-from))
+ (let ((min-strd (consecutive-store-p to))
+ (call-fortran? (> (number-of-elements to) ,fortran-lb)))
+ (cond
+ ((and min-strd call-fortran?)
+ (let ((num-array (,(getf opt :store-allocator) 1)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) num-array))
+ ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) 1) 'num-array 0)
+ (,blas-func (number-of-elements to) num-from
+ num-array 0
+ (store to) min-strd
+ 0 (head to))))
+ (t
+ (let-typed
+ ((t-sto (store to) :type ,(linear-array-type (getf opt :store-type))))
+ (very-quickly
+ (mod-dotimes (idx (dimensions to))
+ with (linear-sums
+ (t-of (strides to) (head to)))
+ do (let-typed
+ ((val ,(funcall (getf opt :reader) 't-sto 't-of) :type ,(getf opt :element-type)))
+ ,(funcall (getf opt :value-writer) '(+ num-from val) 't-sto 't-of))))))))
+ to)))
+
+;;Tweakable
+(defparameter *real-axpy-fortran-call-lower-bound* 20000
+ "If the size of the array is less than this parameter, the
+ lisp version of axpy is called in order to avoid FFI overheads")
+(generate-typed-axpy! real-typed-axpy! (real-tensor
+ daxpy
+ *real-axpy-fortran-call-lower-bound*))
+(generate-typed-num-axpy! real-typed-num-axpy! (real-tensor
+ daxpy
+ *real-axpy-fortran-call-lower-bound*))
+;;Tweakable
+(defparameter *complex-axpy-fortran-call-lower-bound* 10000
+ "If the size of the array is less than this parameter, the
+ lisp version of axpy is called in order to avoid FFI overheads")
+(generate-typed-axpy! complex-typed-axpy! (complex-tensor
+ zaxpy
+ *complex-axpy-fortran-call-lower-bound*))
+(generate-typed-num-axpy! complex-typed-num-axpy! (complex-tensor
+ zaxpy
+ *complex-axpy-fortran-call-lower-bound*))
;;---------------------------------------------------------------;;
(defgeneric axpy! (alpha x y)
@@ -71,6 +130,10 @@
Y <- alpha * x + y
+ If x is T, then
+
+ Y <- alpha + y
+
Purpose
=======
Same as AXPY except that the result
@@ -82,6 +145,12 @@
(:method ((alpha number) (x complex-tensor) (y real-tensor))
(error 'coercion-error :from 'complex-tensor :to 'real-tensor)))
+(defmethod axpy! ((alpha number) (x (eql t)) (y real-tensor))
+ (real-typed-num-axpy! (coerce-real alpha) y))
+
+(defmethod axpy! ((alpha number) (x (eql t)) (y complex-tensor))
+ (complex-typed-num-axpy! (coerce-complex alpha) y))
+
(defmethod axpy! ((alpha number) (x real-tensor) (y real-tensor))
(real-typed-axpy! (coerce-real alpha) x y))
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index 317c011..24bb2d3 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -28,7 +28,7 @@
(in-package #:matlisp)
-(defmacro generate-typed-copy! (func (tensor-class blas-func))
+(defmacro generate-typed-copy! (func (tensor-class blas-func fortran-lb))
;;Be very careful when using functions generated by this macro.
;;Indexes can be tricky and this has no safety net
;;Use only after checking the arguments for compatibility.
@@ -36,23 +36,30 @@
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
`(defun ,func (from to)
(declare (type ,tensor-class from to))
- (if-let (strd-p (blas-copyable-p from to))
- (,blas-func (number-of-elements from) (store from) (first strd-p) (store to) (second strd-p) (head from) (head to))
- (let ((f-sto (store from))
- (t-sto (store to)))
- (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto))
- (very-quickly
- ;;Can possibly make this faster (x2) by using ,blas-func in one of
- ;;the inner loops, but this is to me messy and as of now unnecessary.
- ;;SBCL can already achieve Fortran-ish speed inside this loop.
- (mod-dotimes (idx (dimensions from))
- with (linear-sums
- (f-of (strides from) (head from))
- (t-of (strides to) (head to)))
- do ,(funcall (getf opt :reader-writer) 'f-sto 'f-of 't-sto 't-of)))))
+ (let ((strd-p (blas-copyable-p from to))
+ (call-fortran? (> (number-of-elements to) ,fortran-lb)))
+ (cond
+ ((and strd-p call-fortran?)
+ (,blas-func (number-of-elements from)
+ (store from) (first strd-p)
+ (store to) (second strd-p)
+ (head from) (head to)))
+ (t
+ (let ((f-sto (store from))
+ (t-sto (store to)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto))
+ (very-quickly
+ ;;Can possibly make this faster (x2) by using ,blas-func in one of
+ ;;the inner loops, but this is to me messy and as of now unnecessary.
+ ;;SBCL can already achieve Fortran-ish speed inside this loop.
+ (mod-dotimes (idx (dimensions from))
+ with (linear-sums
+ (f-of (strides from) (head from))
+ (t-of (strides to) (head to)))
+ do ,(funcall (getf opt :reader-writer) 'f-sto 'f-of 't-sto 't-of)))))))
to)))
-(defmacro generate-typed-num-copy! (func (tensor-class blas-func))
+(defmacro generate-typed-num-copy! (func (tensor-class blas-func fortran-lb))
;;Be very careful when using functions generated by this macro.
;;Indexes can be tricky and this has no safety net
;;(you don't see a matrix-ref do you ?)
@@ -62,35 +69,63 @@
`(defun ,func (num-from to)
(declare (type ,tensor-class to)
(type ,(getf opt :element-type) num-from))
- (let ((t-dims (dimensions to))
- (t-stds (strides to))
- (t-sto (store to))
- (t-hd (head to)))
- (declare (type (index-array *) t-dims t-stds)
- (type index-type t-hd)
- (type ,(linear-array-type (getf opt :store-type)) t-sto))
- (if-let (min-stride (consecutive-store-p to))
- (let ((num-array (,(getf opt :store-allocator) 1)))
- (declare (type ,(linear-array-type (getf opt :store-type)) num-array))
- ,(funcall (getf opt :value-writer) 'num-from 'num-array 0)
- (,blas-func (number-of-elements to) num-array 0 t-sto min-stride 0 t-hd))
- (very-quickly
- ;;Can possibly make this faster (x2) by using ,blas-func in one of
- ;;the inner loops, but this is to me messy and as of now unnecessary.
- ;;SBCL can already achieve Fortran-ish speed inside this loop.
- (mod-dotimes (idx t-dims)
- with (linear-sums
- (t-of t-stds t-hd))
- do ,(funcall (getf opt :value-writer) 'num-from 't-sto 't-of))))
- to))))
-
-(generate-typed-copy! real-typed-copy! (real-tensor dcopy))
-(generate-typed-num-copy! real-typed-num-copy! (real-tensor dcopy))
-
-(generate-typed-copy! complex-typed-copy! (complex-tensor zcopy))
-(generate-typed-num-copy! complex-typed-num-copy! (complex-tensor zcopy))
+ (let ((min-stride (consecutive-store-p to))
+ (call-fortran? (> (number-of-elements to) ,fortran-lb)))
+ (cond
+ ((and min-stride call-fortran?)
+ (let ((num-array (,(getf opt :store-allocator) 1)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) num-array))
+ ,(funcall (getf opt :value-writer) 'num-from 'num-array 0)
+ (,blas-func (number-of-elements to)
+ num-array 0
+ (store to) min-stride
+ 0 (head to))))
+ (t
+ (let-typed
+ ((t-sto (store to) :type ,(linear-array-type (getf opt :store-type))))
+ (very-quickly
+ (mod-dotimes (idx (dimensions to))
+ with (linear-sums
+ (t-of (strides to) (head to)))
+ do ,(funcall (getf opt :value-writer) 'num-from 't-sto 't-of)))))))
+ to)))
+
+
+;;Tweakable
+(defparameter *real-copy-fortran-call-lower-bound* 20000
+ "
+ If the dimension of the arguments is less than this parameter,
+ then the Lisp version of copy is used. Default set with SBCL running
+ on x86-64 linux. A reasonable value would be something about 1000.")
+(generate-typed-copy! real-typed-copy! (real-tensor
+ dcopy
+ *real-copy-fortran-call-lower-bound*))
+(generate-typed-num-copy! real-typed-num-copy! (real-tensor
+ dcopy
+ *real-copy-fortran-call-lower-bound*))
+
+;;Tweakable
+(defparameter *complex-copy-fortran-call-lower-bound* 10000
+ "
+ If the dimension of the arguments is less than this parameter,
+ then the Lisp version of copy is used. Default set with SBCL
+ running on x86-64 linux. A reasonable value would be something
+ above 1000.")
+
+(generate-typed-copy! complex-typed-copy! (complex-tensor
+ zcopy
+ *complex-copy-fortran-call-lower-bound*))
+(generate-typed-num-copy! complex-typed-num-copy! (complex-tensor
+ zcopy
+ *complex-copy-fortran-call-lower-bound*))
;;---------------------------------------------------------------;;
+(defun test-copy (n r)
+ (let ((x (make-real-tensor n)))
+ (time (dotimes (i r)
+ (copy! pi x)))
+ t))
+
(defgeneric copy! (from-tensor to-tensor)
(:documentation
"
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index c67af94..911759e 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -1,9 +1,9 @@
(in-package #:matlisp)
(defmacro generate-typed-gemv! (func
- (matrix-class vector-class)
- blas-gemv-func
- fortran-call-lb)
+ (matrix-class vector-class
+ blas-gemv-func
+ fortran-call-lb))
;;Be very careful when using functions generated by this macro.
;;Indexes can be tricky and this has no safety net.
;;Use only after checking the arguments for compatibility.
@@ -68,9 +68,9 @@
MM with small matrices.
Default set with SBCL on x86-64 linux. A reasonable value
is something between 800 and 2000.")
-(generate-typed-gemv! real-typed-gemv!
- (real-matrix real-vector) dgemv
- *real-gemv-fortran-call-lower-bound*)
+(generate-typed-gemv! real-typed-gemv! (real-matrix real-vector
+ dgemv
+ *real-gemv-fortran-call-lower-bound*))
;;Tweakable
(defparameter *complex-gemv-fortran-call-lower-bound* 600
@@ -81,9 +81,9 @@
MM with small matrices.
Default set with SBCL on x86-64 linux. A reasonable value
is something between 400 and 1000.")
-(generate-typed-gemv! complex-typed-gemv!
- (complex-matrix complex-vector) zgemv
- *complex-gemv-fortran-call-lower-bound*)
+(generate-typed-gemv! complex-typed-gemv! (complex-matrix complex-vector
+ zgemv
+ *complex-gemv-fortran-call-lower-bound*))
;;---------------------------------------------------------------;;
;;Can't support "C" because the dual isn't supported by BLAS.
diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp
index d56732e..4cae6f3 100644
--- a/src/level-3/gemm.lisp
+++ b/src/level-3/gemm.lisp
@@ -28,7 +28,7 @@
(in-package #:matlisp)
-(defmacro generate-typed-gemm! (func (matrix-class) (blas-gemm-func blas-gemv-func) fortran-lb-parameter)
+(defmacro generate-typed-gemm! (func (matrix-class blas-gemm-func blas-gemv-func fortran-lb-parameter))
(let* ((opt (get-tensor-class-optimization matrix-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class)
`(defun ,func (alpha A B beta C job)
@@ -161,8 +161,9 @@
MM with small matrices.
Default set with SBCL on x86-64 linux. A reasonable value
is something between 20 and 200.")
-(generate-typed-gemm! real-typed-gemm! (real-matrix) (dgemm dgemv)
- *real-gemm-fortran-call-lower-bound*)
+(generate-typed-gemm! real-typed-gemm! (real-matrix
+ dgemm dgemv
+ *real-gemm-fortran-call-lower-bound*))
;;Tweakable
(defparameter *complex-gemm-fortran-call-lower-bound* 60
@@ -173,8 +174,9 @@
MM with small matrices.
Default set with SBCL on x86-64 linux. A reasonable value
is something between 20 and 200.")
-(generate-typed-gemm! complex-typed-gemm! (complex-matrix) (zgemm zgemv)
- *complex-gemm-fortran-call-lower-bound*)
+(generate-typed-gemm! complex-typed-gemm! (complex-matrix
+ zgemm zgemv
+ *complex-gemm-fortran-call-lower-bound*))
;;---------------------------------------------------------------;;
(defgeneric gemm! (alpha A B beta C &optional job)
diff --git a/src/old/complex-matrix.lisp b/src/old/complex-matrix.lisp
deleted file mode 100644
index a7bc004..0000000
--- a/src/old/complex-matrix.lisp
+++ /dev/null
@@ -1,290 +0,0 @@
-;;; Definitions of COMPLEX-MATRIX.
-
-(in-package :matlisp)
-
-(eval-when (load eval compile)
-(deftype complex-matrix-element-type ()
- "The type of the elements stored in a COMPLEX-MATRIX"
- 'double-float)
-
-(deftype complex-matrix-store-type (size)
- "The type of the storage structure for a COMPLEX-MATRIX"
- `(simple-array double-float (,size)))
-
-(deftype complex-double-float ()
- '(cl:complex (double-float * *)))
-)
-
-;;
-(declaim (inline complex-coerce)
- (ftype (function (number) (complex complex-matrix-element-type))
- complex-coerce))
-(defun complex-coerce (val)
- "
- Syntax
- ======
- (COMPLEX-COERCE number)
-
- Purpose
- =======
- Coerce NUMBER to a complex number.
-"
- (declare (type number val))
- (typecase val
- ((complex complex-matrix-element-type) val)
- (complex (complex (coerce (realpart val) 'complex-matrix-element-type)
- (coerce (imagpart val) 'complex-matrix-element-type)))
- (t (complex (coerce val 'complex-matrix-element-type) 0.0d0))))
-
-;;
-(defclass complex-matrix (standard-matrix)
- ((store
- :initform nil
- :type (complex-matrix-store-type *)))
- (:documentation "A class of matrices with complex elements."))
-
-(defclass sub-complex-matrix (complex-matrix)
- ((parent-matrix
- :initarg :parent
- :accessor parent
- :type complex-matrix))
- (:documentation "A class of matrices with complex elements."))
-
-;;
-(defmethod initialize-instance ((matrix complex-matrix) &rest initargs)
- (setf (store-size matrix) (/ (length (getf :store initargs)) 2))
- (call-next-method))
-
-;;
-(defmethod matrix-ref-1d ((matrix complex-matrix) (idx fixnum))
- (let ((store (store matrix)))
- (declare (type (complex-matrix-store-type *) store))
- (complex (aref store (* 2 idx)) (aref store (+ 1 (* 2 idx))))))
-
-(defmethod (setf matrix-ref-1d) ((value number) (matrix complex-matrix) (idx fixnum))
- (let ((store (store matrix))
- (coerced-value (complex-coerce value)))
- (declare (type (complex-matrix-store-type *) store))
- (setf (aref store (* 2 idx)) (realpart coerced-value)
- (aref store (+ 1 (* 2 idx))) (imagpart coerced-value))))
-
-;;
-(declaim (inline allocate-complex-store))
-(defun allocate-complex-store (size)
- (make-array (* 2 size) :element-type 'complex-matrix-element-type
- :initial-element (coerce 0 'complex-matrix-element-type)))
-
-;;
-(defmethod fill-matrix ((matrix complex-matrix) (fill number))
- (copy! fill matrix))
-
-;;
-(defun make-complex-matrix-dim (n m &key (fill #c(0.0d0 0.0d0)) (order :row-major))
- "
- Syntax
- ======
- (MAKE-COMPLEX-MATRIX-DIM n m {fill-element #C(0d0 0d0)} {order :row-major})
-
- Purpose
- =======
- Creates an NxM COMPLEX-MATRIX with initial contents FILL-ELEMENT,
- the default #c(0.0d0 0.0d0), in the row-major order by default.
-
- See MAKE-COMPLEX-MATRIX.
-"
- (declare (type fixnum n m))
- (let* ((size (* n m))
- (store (allocate-complex-store size)))
- (multiple-value-bind (row-stride col-stride)
- (ecase order
- (:row-major (values m 1))
- (:col-major (values 1 n)))
- (let ((matrix
- (make-instance 'complex-matrix
- :nrows n :ncols m
- :row-stride row-stride :col-stride col-stride
- :store store)))
- (fill-matrix matrix fill)
- matrix))))
-
-;;
-(defun make-complex-matrix-array (array &key (order :row-major))
- "
- Syntax
- ======
- (MAKE-COMPLEX-MATRIX-ARRAY array {order :row-major})
-
- Purpose
- =======
- Creates a COMPLEX-MATRIX with the same contents as ARRAY,
- in row-major order by default.
-"
- (let* ((n (array-dimension array 0))
- (m (array-dimension array 1))
- (size (* n m))
- (store (allocate-complex-store size)))
- (declare (type fixnum n m size)
- (type (complex-matrix-store-type *) store))
- (multiple-value-bind (row-stride col-stride)
- (ecase order
- (:row-major (values m 1))
- (:col-major (values 1 n)))
- (dotimes (i n)
- (declare (type fixnum i))
- (dotimes (j m)
- (declare (type fixnum j))
- (let* ((val (complex-coerce (aref array i j)))
- (realpart (realpart val))
- (imagpart (imagpart val))
- (index (* 2 (store-indexing i j 0 row-stride col-stride))))
- (declare (type complex-matrix-element-type realpart imagpart)
- (type (complex complex-matrix-element-type) val)
- (type fixnum index))
- (setf (aref store index) realpart)
- (setf (aref store (1+ index)) imagpart))))
- (make-instance 'complex-matrix
- :nrows n :ncols m
- :row-stride row-stride :col-stride col-stride
- :store store))))
-
-;;
-(defun make-complex-matrix-seq-of-seq (seq &key (order :row-major))
- (let* ((n (length seq))
- (m (length (elt seq 0)))
- (size (* n m))
- (store (allocate-complex-store size)))
- (declare (type fixnum n m size)
- (type (complex-matrix-store-type *) store))
- (multiple-value-bind (row-stride col-stride)
- (ecase order
- (:row-major (values m 1))
- (:col-major (values 1 n)))
- (dotimes (i n)
- (declare (type fixnum i))
- (let ((this-row (elt seq i)))
- (unless (= (length this-row) m)
- (error "Number of columns is not the same for all rows!"))
- (dotimes (j m)
- (declare (type fixnum j))
- (let* ((val (complex-coerce (elt this-row j)))
- (realpart (realpart val))
- (imagpart (imagpart val))
- (index (* 2 (store-indexing i j 0 row-stride col-stride))))
- (declare (type complex-matrix-element-type realpart imagpart)
- (type (complex complex-matrix-element-type) val)
- (type fixnum index))
- (setf (aref store index) realpart)
- (setf (aref store (1+ index)) imagpart)))))
- (make-instance 'complex-matrix
- :nrows n :ncols m
- :row-stride row-stride :col-stride col-stride
- :store store))))
-
-;;
-(defun make-complex-matrix-seq (seq &key (order :row-major))
- (let* ((n (length seq))
- (store (allocate-complex-store n)))
- (declare (type fixnum n)
- (type (complex-matrix-store-type *) store))
- (dotimes (k n)
- (declare (type fixnum k))
- (let* ((val (complex-coerce (elt seq k)))
- (realpart (realpart val))
- (imagpart (imagpart val))
- (index (* 2 k)))
- (declare (type complex-matrix-element-type realpart imagpart)
- (type (complex complex-matrix-element-type) val)
- (type fixnum index))
- (setf (aref store index) realpart)
- (setf (aref store (1+ index)) imagpart)))
-
- (ecase order
- (:row-major (make-instance 'complex-matrix
- :nrows 1 :ncols n
- :row-stride n :col-stride 1
- :store store))
- (:col-major (make-instance 'complex-matrix
- :nrows n :ncols 1
- :row-stride 1 :col-stride n
- :store store)))))
-
-;;
-(defun make-complex-matrix-sequence (seq &key (order :row-major))
- (cond ((or (listp seq) (vectorp seq))
- (let ((peek (elt seq 0)))
- (cond ((or (listp peek) (vectorp peek))
- ;; We have a seq of seqs
- (make-complex-matrix-seq-of-seq seq :order order))
- (t
- ;; Assume a simple sequence
- (make-complex-matrix-seq seq :order order)))))
- ((arrayp seq)
- (make-complex-matrix-array seq :order order))))
-
-;;
-(defun make-complex-matrix (&rest args)
- "
- Syntax
- ======
- (MAKE-COMPLEX-MATRIX {arg}*)
-
- Purpose
- =======
- Create a FLOAT-MATRIX.
-
- Examples
- ========
-
- (make-complex-matrix n)
- square NxN matrix
- (make-complex-matrix n m)
- NxM matrix
- (make-complex-matrix '((1 2 3) (4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
-
- (make-complex-matrix #((1 2 3) (4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
-
- (make-complex-matrix #((1 2 3) #(4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
-
- (make-complex-matrix #2a((1 2 3) (4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
-
-"
- (let ((nargs (length args)))
- (case nargs
- (1
- (let ((arg (first args)))
- (typecase arg
- (integer
- (assert (not (minusp arg)) nil
- "matrix dimension must be non-negative, not ~A" arg)
- (make-complex-matrix-dim arg arg))
- (sequence
- (make-complex-matrix-sequence arg))
- ((array * (* *))
- (make-complex-matrix-array arg))
- (t (error "don't know how to make matrix from ~a" arg)))))
- (2
- (destructuring-bind (n m)
- args
- (assert (and (typep n '(integer 0))
- (typep n '(integer 0)))
- nil
- "cannot make a ~A x ~A matrix" n m)
- (make-complex-matrix-dim n m)))
- (t
- (error "require 1 or 2 arguments to make a matrix")))))
diff --git a/src/old/mplus.lisp b/src/old/mplus.lisp
index 3e916f7..bbe1229 100644
--- a/src/old/mplus.lisp
+++ b/src/old/mplus.lisp
@@ -25,52 +25,7 @@
;;; ENHANCEMENTS, OR MODIFICATIONS.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Originally written by Raymond Toy.
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; $Id: mplus.lisp,v 1.7 2011/01/25 18:36:56 rtoy Exp $
-;;;
-;;; $Log: mplus.lisp,v $
-;;; Revision 1.7 2011/01/25 18:36:56 rtoy
-;;; Merge changes from automake-snapshot-2011-01-25-1327 to get the new
-;;; automake build infrastructure.
-;;;
-;;; Revision 1.6.2.1 2011/01/25 18:16:53 rtoy
-;;; Use cl:real instead of real.
-;;;
-;;; Revision 1.6 2004/05/24 16:34:22 rtoy
-;;; More SBCL support from Robert Sedgewick. The previous SBCL support
-;;; was incomplete.
-;;;
-;;; Revision 1.5 2002/07/29 01:06:59 rtoy
-;;; Don't use *1x1-real-array*.
-;;;
-;;; Revision 1.4 2000/07/11 18:02:03 simsek
-;;; o Added credits
-;;;
-;;; Revision 1.3 2000/07/11 02:11:56 simsek
-;;; o Added support for Allegro CL
-;;;
-;;; Revision 1.2 2000/05/08 17:19:18 rtoy
-;;; Changes to the STANDARD-MATRIX class:
-;;; o The slots N, M, and NXM have changed names.
-;;; o The accessors of these slots have changed:
-;;; NROWS, NCOLS, NUMBER-OF-ELEMENTS
-;;; The old names aren't available anymore.
-;;; o The initargs of these slots have changed:
-;;; :nrows, :ncols, :nels
-;;;
-;;; Revision 1.1 2000/04/14 00:11:12 simsek
-;;; o This file is adapted from obsolete files 'matrix-float.lisp'
-;;; 'matrix-complex.lisp' and 'matrix-extra.lisp'
-;;; o Initial revision.
-;;;
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(in-package "MATLISP")
+(in-package #:matlisp)
(defgeneric m+ (a b)
(:documentation
@@ -84,50 +39,51 @@
Create a new matrix which is the sum of A and B.
A or B (but not both) may be a scalar, in which
case the addition is element-wise.
-"))
-
-(defgeneric m+! (a b)
- (:documentation
- "
+")
+ (:method ((a number) (b number))
+ (+ a b))
+ (:method ((a standard-tensor) (b standard-tensor))
+ (axpy 1 a b))
+ (:method (
+
+(definline m.+ (a b)
+ "
Syntax
======
- (M+! a b)
+ (M.+ a b)
Purpose
=======
- Desctructive version of M+:
-
- B <- A + B
-"))
+ Same as M+
+"
+ (m+ a b))
-(defgeneric m.+ (a b)
+(defgeneric m+! (a b)
(:documentation
"
Syntax
======
- (M.+ a b)
+ (M+! a b)
Purpose
=======
- Same as M+
-"))
+ Desctructive version of M+:
-(defmethod m.+ (a b)
- (m+ a b))
+ B <- A + B
+")
+ (:method ((a number) (b number))
+ (+ a b)))
-(defgeneric m.+! (a b)
- (:documentation
- "
+(definline m.+! (a b)
+ "
Syntax
======
(M.+! a b)
Purpose
=======
- Same as M.+!
-"))
-
-(defmethod m.+! (a b)
+ Same as M+!
+"
(m+! a b))
(defmethod m+ :before ((a standard-matrix) (b standard-matrix))
diff --git a/src/old/real-matrix.lisp b/src/old/real-matrix.lisp
deleted file mode 100644
index 38ad1f4..0000000
--- a/src/old/real-matrix.lisp
+++ /dev/null
@@ -1,228 +0,0 @@
-
-;;
-(defmethod matrix-ref-1d ((matrix real-matrix) (idx fixnum))
- (let ((store (store matrix)))
- (declare (type (real-matrix-store-type *) store))
- (aref store idx)))
-
-(defmethod (setf matrix-ref-1d) ((value cl:real) (matrix real-matrix) (idx fixnum))
- (let ((store (store matrix)))
- (declare (type (real-matrix-store-type *) store))
- (setf (aref store idx) (coerce value 'double-float))))
-
-;;
-(declaim (inline allocate-real-store))
-(defun allocate-real-store (size &optional (initial-element 0d0))
- (make-array size :element-type 'real-matrix-element-type
- :initial-element (coerce initial-element 'real-matrix-element-type)))
-
-;;
-(defmethod fill-matrix ((matrix real-matrix) (fill cl:real))
- (copy! fill matrix))
-
-(defmethod fill-matrix ((matrix real-matrix) (fill complex))
- (error "cannot fill a real matrix with a complex number,
-don't know how to coerce COMPLEX to REAL"))
-
-;;
-
-;;
-
-;;
-(defun make-real-matrix-dim (n m &key (fill 0.0d0) (order :row-major))
- "
- Syntax
- ======
- (MAKE-REAL-MATRIX-DIM n m [fill-element])
-
- Purpose
- =======
- Creates an NxM REAL-MATRIX with initial contents FILL-ELEMENT,
- the default 0.0d0
-
- See MAKE-REAL-MATRIX.
-"
- (declare (type fixnum n m))
- (let ((casted-fill
- (typecase fill
- (real-matrix-element-type fill)
- (cl:real (coerce fill 'real-matrix-element-type))
- (t (error "argument FILL-ELEMENT to MAKE-REAL-MATRIX-DIM must be a REAL")))))
- (declare (type real-matrix-element-type casted-fill))
- (multiple-value-bind (row-stride col-stride)
- (ecase order
- (:row-major (values m 1))
- (:col-major (values 1 n)))
- (make-instance 'real-matrix
- :nrows n :ncols m
- :row-stride row-stride :col-stride col-stride
- :store (allocate-real-store (* n m) casted-fill)))))
-
-;;; Make a matrix from a 2-D Lisp array
-(defun make-real-matrix-array (array &key (order :row-major))
- "
- Syntax
- ======
- (MAKE-REAL-MATRIX-ARRAY array)
-
- Purpose
- =======
- Creates a REAL-MATRIX with the same contents as ARRAY.
-"
- (let* ((n (array-dimension array 0))
- (m (array-dimension array 1))
- (size (* n m))
- (store (allocate-real-store size)))
- (declare (type fixnum n m size)
- (type (real-matrix-store-type *) store))
- (multiple-value-bind (row-stride col-stride)
- (ecase order
- (:row-major (values m 1))
- (:col-major (values 1 n)))
- (dotimes (i n)
- (declare (type fixnum i))
- (dotimes (j m)
- (declare (type fixnum j))
- (setf (aref store (store-indexing i j 0 row-stride col-stride))
- (coerce (aref array i j) 'real-matrix-element-type))))
- (make-instance 'real-matrix
- :nrows n :ncols m
- :row-stride row-stride :col-stride col-stride
- :store store))))
-
-;;
-(defun make-real-matrix-seq-of-seq (seq &key (order :row-major))
- (let* ((n (length seq))
- (m (length (elt seq 0)))
- (size (* n m))
- (store (allocate-real-store size)))
- (declare (type fixnum n m size)
- (type (real-matrix-store-type *) store))
- (multiple-value-bind (row-stride col-stride)
- (ecase order
- (:row-major (values m 1))
- (:col-major (values 1 n)))
- (dotimes (i n)
- (declare (type fixnum i))
- (let ((this-row (elt seq i)))
- (unless (= (length this-row) m)
- (error "Number of columns is not the same for all rows!"))
- (dotimes (j m)
- (declare (type fixnum j))
- (setf (aref store (store-indexing i j 0 row-stride col-stride))
- (coerce (elt this-row j) 'real-matrix-element-type)))))
- (make-instance 'real-matrix
- :nrows n :ncols m
- :row-stride row-stride :col-stride col-stride
- :store store))))
-
-;;
-(defun make-real-matrix-seq (seq &key (order :row-major))
- (let* ((n (length seq))
- (store (allocate-real-store n)))
- (declare (type fixnum n))
- (dotimes (k n)
- (declare (type fixnum k))
- (setf (aref store k) (coerce (elt seq k) 'real-matrix-element-type)))
- (ecase order
- (:row-major (make-instance 'real-matrix
- :nrows 1 :ncols n
- :row-stride n :col-stride 1
- :store store))
- (:col-major (make-instance 'real-matrix
- :nrows n :ncols 1
- :row-stride 1 :col-stride n
- :store store)))))
-
-
-;;
-(defun make-real-matrix-sequence (seq &key (order :row-major))
- (cond ((or (listp seq) (vectorp seq))
- (let ((peek (elt seq 0)))
- (cond ((or (listp peek) (vectorp peek))
- ;; We have a seq of seqs
- (make-real-matrix-seq-of-seq seq :order order))
- (t
- ;; Assume a simple sequence
- (make-real-matrix-seq seq :order order)))))
- ((arrayp seq)
- (make-real-matrix-array seq :order order))))
-
-;;
-(defun make-real-matrix (&rest args)
- "
- Syntax
- ======
- (MAKE-REAL-MATRIX {arg}*)
-
- Purpose
- =======
- Create a REAL-MATRIX.
-
- Examples
- ========
-
- (make-real-matrix n)
- square NxN matrix
- (make-real-matrix n m)
- NxM matrix
- (make-real-matrix '((1 2 3) (4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
-
- (make-real-matrix #((1 2 3) (4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
-
- (make-real-matrix #((1 2 3) #(4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
-
- (make-real-matrix #2a((1 2 3) (4 5 6)))
- 2x3 matrix:
-
- 1 2 3
- 4 5 6
- (make-real-matrix #(1 2 3 4))
- 4x1 matrix (column vector)
-
- 1
- 2
- 3
- 4
-
- (make-real-matrix #((1 2 3 4))
- 1x4 matrix (row vector)
-
- 1 2 3 4
-"
- (let ((nargs (length args)))
- (case nargs
- (1
- (let ((arg (first args)))
- (typecase arg
- (integer
- (assert (not (minusp arg)) nil
- "matrix dimension must be positive, not ~A" arg)
- (make-real-matrix-dim arg arg))
- (sequence
- (make-real-matrix-sequence arg))
- ((array * (* *))
- (make-real-matrix-array arg))
- (t (error "don't know how to make matrix from ~a" arg)))))
- (2
- (destructuring-bind (n m)
- args
- (assert (and (typep n '(integer 0))
- (typep n '(integer 0)))
- nil
- "cannot make a ~A x ~A matrix" n m)
- (make-real-matrix-dim n m)))
- (t
- (error "require 1 or 2 arguments to make a matrix")))))
-----------------------------------------------------------------------
Summary of changes:
src/level-1/axpy.lisp | 113 ++++++++++++++----
src/level-1/copy.lisp | 117 +++++++++++------
src/level-2/gemv.lisp | 18 ++--
src/level-3/gemm.lisp | 12 +-
src/old/complex-matrix.lisp | 290 -------------------------------------------
src/old/mplus.lisp | 92 ++++----------
src/old/real-matrix.lisp | 228 ---------------------------------
7 files changed, 207 insertions(+), 663 deletions(-)
delete mode 100644 src/old/complex-matrix.lisp
delete mode 100644 src/old/real-matrix.lisp
hooks/post-receive
--
matlisp
|
