Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch matlisp-cffi updated. 2012-02-23-40-g2cbbcb6
|
From: Akshay S. <ak...@us...> - 2012-03-12 17:28:26
|
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, matlisp-cffi has been updated
via 2cbbcb64c997e7ed95c2e9a344347ca9606fce66 (commit)
via b0f1cb2dd42338c9189c83cbcbcb177eaf1c7845 (commit)
via 98cdecd68f57b4a561ac8f68a0ede4a0374a6a95 (commit)
via 097c9251aaa15f702cc8b9701832dbdc1d9bcb13 (commit)
via b0415f8f7e3a4af5682f22aed1885e20f3484188 (commit)
from 96f2abfd9395d6f25520c3b828c4c69a0f35d8a3 (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 2cbbcb64c997e7ed95c2e9a344347ca9606fce66
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Mar 12 22:48:55 2012 +0530
-> Support for copy(!) and scal(!) are reasonably complete.
-> Restricted the store type of {real/complex}-matrix to one-dimension.
-> Usual churning in utilities.lisp, made mlet* even more incomprehensible (but nicer!).
diff --git a/packages.lisp b/packages.lisp
index 383af89..566dbd6 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -298,6 +298,7 @@
"*PRINT-MATRIX*"
"AXPY!"
"AXPY"
+ "BLAS-COPYABLE-P"
"COL-VECTOR-P"
"COMPLEX-COERCE"
"COMPLEX-MATRIX"
@@ -420,6 +421,7 @@
"SWAP!"
"TR"
"TRANSPOSE"
+ "TRANSPOSE!"
"VEC"
"UNLOAD-BLAS-&-LAPACK-LIBRARIES"
"ZEROS"
diff --git a/src/axpy.lisp b/src/axpy.lisp
index c3a6baf..57104cf 100644
--- a/src/axpy.lisp
+++ b/src/axpy.lisp
@@ -150,8 +150,8 @@
(store-y (store y))
(store-result (store result)))
(declare (type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store-y)
- (type (complex-matrix-store-type (*)) store-x store-result))
+ (type (real-matrix-store-type *) store-y)
+ (type (complex-matrix-store-type *) store-x store-result))
(zcopy nxm store-x 1 store-result 1) ;; same as (COPY! x result)
(zdscal nxm alpha store-result 1) ;; same as (SCAL! alpha result)
@@ -192,8 +192,8 @@
(c-alpha (complex-coerce alpha)))
(declare (type complex-double-float c-alpha)
(type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store-x)
- (type (complex-matrix-store-type (*)) store-y store-result))
+ (type (real-matrix-store-type *) store-x)
+ (type (complex-matrix-store-type *) store-y store-result))
(dcopy nxm store-x 1 store-result 2)
(zscal nxm c-alpha store-result 1)
@@ -209,7 +209,7 @@
(c-alpha (complex-coerce alpha)))
(declare (type complex-double-float c-alpha)
(type fixnum nxm)
- (type (complex-matrix-store-type (*)) store-result))
+ (type (complex-matrix-store-type *) store-result))
(zscal nxm c-alpha store-result 1)
(daxpy nxm 1.0d0 (store y) 1 store-result 2)
@@ -283,8 +283,8 @@ don't know how to coerce COMPLEX to REAL"))
(c-alpha (complex-coerce alpha)))
(declare (type complex-double-float c-alpha)
(type fixnum nxm)
- (type (real-matrix-store-type (*)) store-x)
- (type (complex-matrix-store-type (*)) store-y))
+ (type (real-matrix-store-type *) store-x)
+ (type (complex-matrix-store-type *) store-y))
(daxpy nxm (realpart c-alpha) store-x 1 store-y 2)
(with-vector-data-addresses ((addr-y store-y)
diff --git a/src/blas.lisp b/src/blas.lisp
index 999b389..964f7b0 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -1,4 +1,4 @@
-t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
@@ -204,7 +204,7 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
"
(n :integer :input)
(da :double-float :input)
- (dx (* :double-float :inc head-dx) :output)
+ (dx (* :double-float :inc head-x) :output)
(incx :integer :input)
)
@@ -363,7 +363,7 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
"
(n :integer :input)
(da :double-float :input)
- (zx (* :complex-double-float :inc head-zx) :output)
+ (zx (* :complex-double-float :inc head-x) :output)
(incx :integer :input)
)
@@ -405,7 +405,7 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
"
(n :integer :input)
(za :complex-double-float)
- (zx (* :complex-double-float) :output)
+ (zx (* :complex-double-float :inc head-x) :output)
(incx :integer :input)
)
@@ -530,7 +530,8 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
considered in the operations are:
Y(0),Y(2*INCY), ... , Y(2*(N-1)*INCY)
-" (n :integer :input)
+"
+ (n :integer :input)
(zx (* :complex-double-float) :input)
(incx :integer :input)
(zy (* :complex-double-float) :input)
@@ -547,6 +548,11 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
(def-fortran-routine dasum :double-float
"
+ Purpose
+ =======
+
+ Takes the sum of the absolute values.
+
"
(n :integer :input)
(dx (* :double-float) :input)
diff --git a/src/compat.lisp b/src/compat.lisp
index 4b0f329..b30383e 100644
--- a/src/compat.lisp
+++ b/src/compat.lisp
@@ -93,11 +93,11 @@
(deftype float-matrix-array-type (size)
"Defines the same type as (REAL-MATRIX-STORE-TYPE (*))"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
(deftype complex-matrix-array-type (size)
"Defines the same type as (COMPLEX-MATRIX-STORE-TYPE (*))"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
(deftype float-matrix ()
"Defines the same type as REAL-MATRIX"
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 48e7d0f..64382d7 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -9,7 +9,7 @@
(deftype complex-matrix-store-type (size)
"The type of the storage structure for a COMPLEX-MATRIX"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
)
;;
@@ -37,7 +37,7 @@
(defclass complex-matrix (standard-matrix)
((store
:initform nil
- :type (simple-array complex-matrix-element-type (*))))
+ :type (complex-matrix-store-type *)))
(:documentation "A class of matrices with complex elements."))
;;
@@ -55,13 +55,13 @@
;;
(defmethod matrix-ref-1d ((matrix complex-matrix) (idx fixnum))
(let ((store (store matrix)))
- (declare (type (complex-matrix-store-type (*)) store))
+ (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))
+ (declare (type (complex-matrix-store-type *) store))
(setf (aref store (* 2 idx)) (realpart coerced-value)
(aref store (+ 1 (* 2 idx))) (imagpart coerced-value))))
@@ -121,7 +121,7 @@
(size (* n m))
(store (allocate-complex-store size)))
(declare (type fixnum n m size)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
@@ -151,7 +151,7 @@
(size (* n m))
(store (allocate-complex-store size)))
(declare (type fixnum n m size)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
@@ -182,7 +182,7 @@
(let* ((n (length seq))
(store (allocate-complex-store n)))
(declare (type fixnum n)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(dotimes (k n)
(declare (type fixnum k))
(let* ((val (complex-coerce (elt seq k)))
diff --git a/src/copy.lisp b/src/copy.lisp
index 586aaf0..68dff10 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -79,29 +79,6 @@
(in-package "MATLISP")
;;
-(defun blas-copyable-p (matrix)
- (declare (optimize (safety 0) (speed 3))
- (type (or real-matrix complex-matrix) matrix))
- (mlet* ((nr (nrows matrix) :type fixnum)
- (nc (ncols matrix) :type fixnum)
- (rs (row-stride matrix) :type fixnum)
- (cs (col-stride matrix) :type fixnum)
- (ne (number-of-elements matrix) :type fixnum))
- (cond
- ((= nc 1) (values t rs ne))
- ((= nr 1) (values t cs ne))
- ((= rs (* nc cs)) (values t cs ne))
- ((= cs (* nr rs)) (values t rs ne))
- (t (values nil -1 -1)))))
-
-;;
-(defmacro with-transpose! (matlst &rest body)
- `(progn
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
- ,@body
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
-
-;;
(defmacro generate-typed-copy!-func (func store-type matrix-type blas-func)
`(defun ,func (mat-a mat-b)
(declare (type ,matrix-type mat-a mat-b)
@@ -111,7 +88,7 @@
((hd-a st-a sz) (slot-values mat-a '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum))
((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *))))
(if (and cp-a cp-b)
- (,blas-func sz (store mat-a) inc-a (store mat-b) inc-b :head-x hd-a :head-y hd-b)
+ (,blas-func sz st-a inc-a st-b inc-b :head-x hd-a :head-y hd-b)
(symbol-macrolet
((common-code
(mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
@@ -128,8 +105,33 @@
mat-b)))
;;
-(defvar *1x1-real-array* (make-array 1 :element-type 'double-float))
-(defvar *1x1-complex-array* (make-array 2 :element-type 'double-float))
+(defmacro generate-typed-num-copy!-func (func element-type store-type matrix-type blas-func
+ array-decl)
+ (let ((num-arg (car array-decl)))
+ (destructuring-bind (var var-maker-form var-setf-form &key type) (cadr array-decl)
+ `(mlet* (((,var) (,@var-maker-form) ,@(if type
+ `(:type (,type)))))
+ (defun ,func (,num-arg mat-x)
+ (declare (type ,element-type ,num-arg)
+ (type ,matrix-type mat-x)
+ (optimize (safety 0) (speed 3)))
+ (,@var-setf-form)
+ (mlet* (((cp-x inc-x sz-x) (blas-copyable-p mat-x) :type (boolean fixnum nil))
+ ((hd-x st-x sz) (slot-values mat-x '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum)))
+ (if cp-x
+ (,blas-func sz ,var 0 st-x inc-x :head-y hd-x)
+ (symbol-macrolet
+ ((common-code
+ (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum)))
+ (loop for i from 0 below nr-x
+ do (,blas-func nc-x ,var 0 st-x cs-x :head-y (+ hd-x (* i rs-x)))))))
+ ;;Choose the smaller of the loops
+ (if (> (nrows mat-x) (ncols mat-x))
+ (with-transpose! (mat-x)
+ common-code)
+ common-code)))
+ mat-x))))))
;;
(defgeneric copy! (matrix new-matrix)
@@ -157,102 +159,72 @@
"))
(defmethod copy! :before ((x standard-matrix) (y standard-matrix))
- (let ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y)))
- (declare (type fixnum nxm-x nxm-y))
- (if (not (= nxm-x nxm-y))
- (warn "arguments X,Y to COPY! are of different sizes"))))
+ (mlet* (((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols)) :type (fixnum fixnum))
+ ((nr-y nc-y) (slot-values y '(number-of-rows number-of-cols)) :type (fixnum fixnum)))
+ (unless (and (= nr-x nr-y) (= nc-x nc-y))
+ (error "Arguments X,Y to COPY! are of different dimensions."))))
;;
-(generate-typed-copy!-func real-double-copy!-typed real-matrix-store-type real-matrix blas:dcopy)
-
-(defmethod copy! ((x real-matrix) (y real-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
- (dcopy nxm (store x) 1 (store y) 1)
+(defmethod copy! ((x standard-matrix) (y standard-matrix))
+ (mlet* (((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols))
+ :type (fixnum fixnum)))
+ (dotimes (i nr-x)
+ (dotimes (j nc-x)
+ (declare (type fixnum i j))
+ (setf (matrix-ref-2d y i j) (matrix-ref-2d x i j))))
y))
-(defmethod copy! ((x real-matrix) (y complex-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
+;;
+(generate-typed-copy!-func real-double-copy!-typed real-matrix-store-type real-matrix blas:dcopy)
- ;; Set the imaginary parts of Y to zero.
- (zdscal nxm 0.0d0 (store y) 1)
+(generate-typed-num-copy!-func real-double-num-copy!-typed
+ double-float real-matrix-store-type real-matrix
+ blas:dcopy
+ (num
+ (1x1-array
+ (allocate-real-store 1)
+ (setf (aref 1x1-array 0) num)
+ :type (real-matrix-store-type 1))))
- ;; Copy the elements of X to the real parts of Y.
- (dcopy nxm (store x) 1 (store y) 2)
- y))
-
(defmethod copy! ((x complex-matrix) (y real-matrix))
- (error "cannot copy a COMPLEX-MATRIX into a REAL-MATRIX,
+ (error "Cannot copy a COMPLEX-MATRIX into a REAL-MATRIX,
don't know how to coerce a COMPLEX to a REAL"))
+(defmethod copy! ((x complex) (y real-matrix))
+ (error "Cannot copy ~a to ~a, don't know how to coerce COMPLEX to REAL"
+ x y))
-;;
-(generate-typed-copy!-func complex-double-copy!-typed complex-matrix-store-type complex-matrix blas:zcopy)
-
-(defmethod copy! ((x complex-matrix) (y complex-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
- (dcopy (* 2 nxm) (store x) 1 (store y) 1)
- y))
-
-(defmethod copy! ((x standard-matrix) (y standard-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
-
- (dotimes (i nxm)
- (declare (type fixnum i))
- (setf (matrix-ref y i) (matrix-ref x i)))
-
- y))
-
-(defmethod copy! ((x #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (y real-matrix))
- (let ((nxm (number-of-elements y)))
- (setf (aref *1x1-real-array* 0) x)
- (dcopy nxm *1x1-real-array* 0 (store y) 1)
- y))
+(defmethod copy! ((x real-matrix) (y real-matrix))
+ (real-double-copy!-typed x y))
(defmethod copy! ((x cl:real) (y real-matrix))
- (let ((nxm (number-of-elements y)))
- (setf (aref *1x1-real-array* 0) (coerce x 'real-matrix-element-type))
- (dcopy nxm *1x1-real-array* 0 (store y) 1)
- y))
+ (real-double-num-copy!-typed (coerce x 'double-float) y))
-(defmethod copy! ((x complex) (y real-matrix))
- (error "cannot copy ~a to ~a, don't know how to coerce COMPLEX to REAL"
- x
- y))
+;;
+(generate-typed-copy!-func complex-double-copy!-typed complex-matrix-store-type complex-matrix blas:zcopy)
-(defmethod copy! ((x #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (y complex-matrix))
- (let ((nxm (number-of-elements y)))
+(generate-typed-num-copy!-func complex-double-num-copy!-typed
+ (complex (double-float * *)) complex-matrix-store-type complex-matrix
+ blas:zcopy
+ (num
+ (1x1-z-array
+ (allocate-complex-store 1)
+ (setf (aref 1x1-z-array 0) (realpart num)
+ (aref 1x1-z-array 1) (imagpart num))
+ :type (complex-matrix-store-type 2))))
- #-(or cmu sbcl) (setq x (complex-coerce x))
+(defmethod copy! ((x complex-matrix) (y complex-matrix))
+ (complex-double-copy!-typed x y))
- (setf (aref *1x1-complex-array* 0) (realpart x))
- (setf (aref *1x1-complex-array* 1) (imagpart x))
- (zcopy nxm *1x1-complex-array* 0 (store y) 1)
- y))
+(defmethod copy! ((x real-matrix) (y complex-matrix))
+ (real-double-copy!-typed x (realpart! y))
+ (scal! 0d0 (imagpart! y))
+ y)
(defmethod copy! ((x number) (y complex-matrix))
- (let ((nxm (number-of-elements y)))
- (setq x (complex-coerce x))
- (setf (aref *1x1-complex-array* 0) (realpart x))
- (setf (aref *1x1-complex-array* 1) (imagpart x))
- (zcopy nxm *1x1-complex-array* 0 (store y) 1)
- y))
+ (complex-double-num-copy!-typed (complex-coerce x) y))
-;;
+;;;;
(defgeneric copy (matrix)
(:documentation
"
@@ -264,32 +236,23 @@ don't know how to coerce a COMPLEX to a REAL"))
=======
Return a copy of the matrix X"))
-(defmethod copy ((matrix standard-matrix))
- (make-instance 'standard-matrix :nrows (nrows matrix) :ncols (ncols matrix) :store (copy-seq (store matrix))))
-
(defmethod copy ((matrix real-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
- (m (ncols matrix))
+ (let* ((n (nrows matrix))
+ (m (nrows matrix))
(result (make-real-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:dcopy size (store matrix) 1 (store result) 1)
- result))
+ (declare (type fixnum n m))
+ (copy! matrix result)))
(defmethod copy ((matrix complex-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
+ (let* ((n (nrows matrix))
(m (ncols matrix))
(result (make-complex-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:zcopy size (store matrix) 1 (store result) 1)
- result))
+ (declare (type fixnum n m))
+ (copy! matrix result)))
(defmethod copy ((matrix number))
matrix)
-
-
;;
(defgeneric convert-to-lisp-array (matrix)
(:documentation
diff --git a/src/dot.lisp b/src/dot.lisp
index 4b643c2..badff89 100644
--- a/src/dot.lisp
+++ b/src/dot.lisp
@@ -133,8 +133,8 @@
(store-x (store x))
(store-y (store y)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store-x)
- (type (complex-matrix-store-type (*)) store-y))
+ (type (real-matrix-store-type *) store-x)
+ (type (complex-matrix-store-type *) store-y))
(let ((realpart (ddot nxm store-x 1 store-y 2))
(imagpart (with-vector-data-addresses ((addr-x store-x)
@@ -186,8 +186,8 @@
(store-x (store x))
(store-y (store y)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store-y)
- (type (complex-matrix-store-type (*)) store-x))
+ (type (real-matrix-store-type *) store-y)
+ (type (complex-matrix-store-type *) store-x))
(let ((realpart (ddot nxm store-x 2 store-y 1))
(imagpart (with-vector-data-addresses ((addr-x store-x)
diff --git a/src/map.lisp b/src/map.lisp
index 8206e55..5b2afa4 100644
--- a/src/map.lisp
+++ b/src/map.lisp
@@ -124,7 +124,7 @@
(let ((nxm (number-of-elements mat))
(store (store mat)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (k nxm mat)
(declare (type fixnum k))
(setf (aref store k) (funcall func (aref store k))))))
@@ -143,7 +143,7 @@
(let ((nxm (number-of-elements mat))
(store (store mat)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (k nxm mat)
(declare (type fixnum k))
(setf (aref store k) (funcall func (aref store k))))))
diff --git a/src/norm.lisp b/src/norm.lisp
index b03b512..7087dbe 100644
--- a/src/norm.lisp
+++ b/src/norm.lisp
@@ -128,7 +128,7 @@
(nxm (number-of-elements a))
(store (store a)))
(declare (type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(if (row-or-col-vector-p a)
(case p
@@ -199,7 +199,7 @@
(nxm (number-of-elements a))
(store (store a)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(if (row-or-col-vector-p a)
(case p
diff --git a/src/real-matrix.lisp b/src/real-matrix.lisp
index f6a0e0e..87b6684 100644
--- a/src/real-matrix.lisp
+++ b/src/real-matrix.lisp
@@ -9,14 +9,13 @@
(deftype real-matrix-store-type (size)
"The type of the storage structure for a REAL-MATRIX"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
)
-
;;
(defclass real-matrix (standard-matrix)
((store
:initform nil
- :type (simple-array real-matrix-element-type (*))))
+ :type (real-matrix-store-type *)))
(:documentation "A class of matrices with real elements."))
;;
@@ -34,13 +33,13 @@
;;
(defmethod matrix-ref-1d ((matrix real-matrix) (idx fixnum))
(let ((store (store matrix)))
- (declare (type (real-matrix-store-type (*)) store))
+ (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))
+ (declare (type (real-matrix-store-type *) store))
(setf (aref store idx) value)))
;;
@@ -103,7 +102,7 @@ don't know how to coerce COMPLEX to REAL"))
(size (* n m))
(store (allocate-real-store size)))
(declare (type fixnum n m size)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
@@ -126,7 +125,7 @@ don't know how to coerce COMPLEX to REAL"))
(size (* n m))
(store (allocate-real-store size)))
(declare (type fixnum n m size)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
diff --git a/src/realimag.lisp b/src/realimag.lisp
index 0b7ae11..2b35319 100644
--- a/src/realimag.lisp
+++ b/src/realimag.lisp
@@ -112,8 +112,8 @@
(store (store mat))
(new-store (allocate-real-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store)
- (type (real-matrix-store-type (*)) new-store))
+ (type (complex-matrix-store-type *) store)
+ (type (real-matrix-store-type *) new-store))
(dcopy nxm store 2 new-store 1)
@@ -140,8 +140,8 @@ its element types are unknown"))
(store (store mat))
(new-store (allocate-real-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store)
- (type (real-matrix-store-type (*)) new-store))
+ (type (complex-matrix-store-type *) store)
+ (type (real-matrix-store-type *) new-store))
(with-vector-data-addresses ((addr-store store)
(addr-new-store new-store))
diff --git a/src/ref.lisp b/src/ref.lisp
index ce21f79..c5619ef 100644
--- a/src/ref.lisp
+++ b/src/ref.lisp
@@ -152,7 +152,7 @@
(store (allocate-real-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store mat-store store))
+ (type (real-matrix-store-type *) idx-store mat-store store))
(if (and (> n 1)
(> m 1))
@@ -183,7 +183,7 @@
(store (allocate-real-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) mat-store store))
+ (type (real-matrix-store-type *) mat-store store))
(if (and (> n 1)
(> m 1))
@@ -216,7 +216,7 @@
(mat-store (store mat))
(store (store slice)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store
mat-store
@@ -247,7 +247,7 @@
(store (store slice)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
mat-store
store))
@@ -269,7 +269,7 @@
(store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (i n)
(declare (type fixnum i))
@@ -303,7 +303,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -472,7 +472,7 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store new-store mat-store))
+ (type (real-matrix-store-type *) idx-store new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -496,7 +496,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m)
- (type (real-matrix-store-type (*)) new-store mat-store))
+ (type (real-matrix-store-type *) new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -519,7 +519,7 @@
(new-store (store new))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store
new-store
@@ -545,7 +545,7 @@
(mat-store (store mat)))
(declare (type fixnum l n)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
new-store
mat-store))
@@ -572,7 +572,7 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store mat-store))
+ (type (real-matrix-store-type *) idx-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -594,7 +594,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m)
- (type (real-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -618,7 +618,7 @@
(col-idx-store (store col-idx))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store
mat-store))
@@ -642,7 +642,7 @@
(mat-store (store mat)))
(declare (type fixnum l)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
mat-store))
@@ -676,7 +676,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -775,7 +775,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -838,7 +838,7 @@
(new-store (store new)))
(declare (type fixnum n)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(setf (aref store (fortran-matrix-indexing i 0 n))
(aref new-store 0))))
@@ -888,7 +888,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -1025,7 +1025,7 @@
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1100,7 +1100,7 @@
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1134,7 +1134,7 @@
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1235,8 +1235,8 @@
(store (allocate-complex-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store)
- (type (complex-matrix-store-type (*)) mat-store store))
+ (type (real-matrix-store-type *) idx-store)
+ (type (complex-matrix-store-type *) mat-store store))
(if (and (> n 1)
(> m 1))
@@ -1272,7 +1272,7 @@
(store (allocate-complex-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (complex-matrix-store-type (*)) mat-store store))
+ (type (complex-matrix-store-type *) mat-store store))
(if (and (> n 1)
(> m 1))
@@ -1311,10 +1311,10 @@
(mat-store (store mat))
(store (store slice)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store
store))
@@ -1347,7 +1347,7 @@
(store (store slice)))
(declare (type fixnum l n m)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store
store))
@@ -1375,7 +1375,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1452,7 +1452,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1498,7 +1498,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1595,8 +1595,8 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store)
- (type (complex-matrix-store-type (*)) new-store mat-store))
+ (type (real-matrix-store-type *) idx-store)
+ (type (complex-matrix-store-type *) new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -1624,7 +1624,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m) ;; k)
- (type (complex-matrix-store-type (*)) new-store mat-store))
+ (type (complex-matrix-store-type *) new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -1651,10 +1651,10 @@
(new-store (store new))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
new-store
mat-store))
@@ -1683,7 +1683,7 @@
(mat-store (store mat)))
(declare (type fixnum l n) ;;m)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
new-store
mat-store))
@@ -1715,8 +1715,8 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store new-store)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) idx-store new-store)
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1744,8 +1744,8 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m) ;; k)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1772,11 +1772,11 @@
(new-store (store new))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
new-store
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
(dotimes (i n)
@@ -1804,9 +1804,9 @@
(mat-store (store mat)))
(declare (type fixnum l n) ;; m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
new-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
(let ((i -1))
@@ -1836,8 +1836,8 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) idx-store)
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1864,7 +1864,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m) ;; k)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1892,10 +1892,10 @@
(col-idx-store (store col-idx))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
(dotimes (i n)
@@ -1922,7 +1922,7 @@
(mat-store (store mat)))
(declare (type fixnum l) ;; n m)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
;; Hmm, another redundant i,j (see above)
@@ -1954,7 +1954,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -2059,7 +2059,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -2126,7 +2126,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -2239,8 +2239,8 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) store))
(let ((p (if (integerp i)
@@ -2344,8 +2344,8 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) store))
(let ((p (if (integerp i)
@@ -2395,8 +2395,8 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) store))
(let ((p (if (integerp i)
@@ -2494,7 +2494,7 @@
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
#+(or :ccl :allegro) (setq new (complex-coerce new))
@@ -2583,7 +2583,7 @@
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
#+(or :ccl :allegro) (setq new (complex-coerce new))
@@ -2626,7 +2626,7 @@
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
#+(or :ccl :allegro) (setq new (complex-coerce new))
@@ -2811,7 +2811,7 @@
(m (ncols item))
(store (store item)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (i n)
(declare (type fixnum i))
(dotimes (j m)
diff --git a/src/reshape.lisp b/src/reshape.lisp
index 7321600..9681990 100644
--- a/src/reshape.lisp
+++ b/src/reshape.lisp
@@ -115,7 +115,7 @@
(new-size (* new-nrows new-ncols))
(new-store (allocate-real-store new-size)))
(declare (fixnum old-size new-size)
- (type (real-matrix-store-type (*)) new-store))
+ (type (real-matrix-store-type *) new-store))
(dcopy (min old-size new-size) (store mat) 1 new-store 1)
@@ -132,7 +132,7 @@
;; it's bigger than the old size. Allocate it and copy the
;; elements over.
(let ((new-store (allocate-real-store new-size)))
- (declare (type (real-matrix-store-type (*)) new-store))
+ (declare (type (real-matrix-store-type *) new-store))
(dcopy new-size (store mat) 1 new-store 1)
(setf (slot-value mat 'store) new-store)))
@@ -148,7 +148,7 @@
(new-size (* new-nrows new-ncols))
(new-store (allocate-complex-store new-size)))
(declare (fixnum old-size new-size)
- (type (complex-matrix-store-type (*)) new-store))
+ (type (complex-matrix-store-type *) new-store))
(zcopy (min old-size new-size) (store mat) 1 new-store 1)
(make-instance 'complex-matrix :nrows new-nrows :ncols new-ncols :store new-store)))
@@ -164,7 +164,7 @@
;; it's bigger than the old size. Allocate it and copy the
;; elements over.
(let ((new-store (allocate-complex-store new-size)))
- (declare (type (complex-matrix-store-type (*)) new-store))
+ (declare (type (complex-matrix-store-type *) new-store))
(zcopy new-size (store mat) 1 new-store 1)
(setf (slot-value mat 'store) new-store)))
diff --git a/src/scal.lisp b/src/scal.lisp
index 6d2026e..8548a26 100644
--- a/src/scal.lisp
+++ b/src/scal.lisp
@@ -69,34 +69,34 @@
(in-package "MATLISP")
-#+nil (use-package "BLAS")
-#+nil (use-package "LAPACK")
-#+nil (use-package "FORTRAN-FFI-ACCESSORS")
-
-#+nil (export '(scal!
- scal))
-
-(defgeneric scal (alpha x)
- (:documentation
-"
- Sytnax
- ======
- (SCAL alpha x)
-
- Purpose
- =======
- Computes and returns a new matrix equal to
-
- alpha * X
-
- where alpha is a scalar and X is a matrix.
-
-"))
-
+(defmacro generate-typed-scal!-func (func element-type store-type matrix-type blas-func)
+ `(defun ,func (alpha mat-x)
+ (declare (type ,matrix-type mat-x)
+ (type ,element-type alpha)
+ (optimize (safety 0) (speed 3)))
+ (mlet* (((cp-x inc-x sz-x) (blas-copyable-p mat-x)
+ :type (boolean fixnum fixnum))
+ ((hd-x st-x) (slot-values mat-x '(head store))
+ :type (fixnum (,store-type *))))
+ (if cp-x
+ (,blas-func sz-x alpha st-x inc-x :head-x hd-x)
+ (symbol-macrolet
+ ((common-code
+ (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum)))
+ (loop for i from 0 below nr-x
+ do (,blas-func nc-x alpha st-x cs-x :head-x (+ hd-x (* i rs-x)))))))
+ (if (> (nrows mat-x) (ncols mat-x))
+ (with-transpose! (mat-x)
+ common-code)
+ common-code)))
+ mat-x)))
+
+;;
(defgeneric scal! (alpha x)
(:documentation
"
- Sytnax
+ Syntax
======
(SCAL! alpha x)
@@ -106,117 +106,66 @@
stored in X.
"))
-(defmethod scal ((alpha number) (x number))
- (* alpha x))
-
-(defmethod scal ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix))
- (let ((nxm (number-of-elements x))
- (result (copy x)))
- (declare (type fixnum nxm))
-
- (dscal nxm alpha (store result) 1)
- result))
-
-(defmethod scal ((alpha cl:real) (x real-matrix))
- (scal (coerce alpha 'real-matrix-element-type) x))
-
-(defmethod scal ((alpha #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (x real-matrix))
- (let* ((nxm (number-of-elements x))
- (n (nrows x))
- (m (ncols x))
- (result (make-complex-matrix-dim n m)))
- (declare (type fixnum n m nxm))
-
- #-(or cmu sbcl) (setq alpha (complex-coerce alpha))
-
- (copy! x result)
- (setf (aref *1x1-complex-array* 0) (realpart alpha))
- (setf (aref *1x1-complex-array* 1) (imagpart alpha))
- (zscal nxm *1x1-complex-array* (store result) 1)
-
- result))
-
-#+(or :cmu :sbcl)
-(defmethod scal ((alpha complex) (x real-matrix))
- (scal (complex-coerce alpha) x))
-
-(defmethod scal ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix))
- (let ((nxm (number-of-elements x))
- (result (copy x)))
- (declare (type fixnum nxm))
- (zdscal nxm alpha (store result) 1)
-
- result))
-
-(defmethod scal ((alpha cl:real) (x complex-matrix))
- (scal (coerce alpha 'real-matrix-element-type) x))
-
-(defmethod scal ((alpha #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (x complex-matrix))
- (let ((nxm (number-of-elements x))
- (result (copy x)))
- (declare (type fixnum nxm))
-
- #-(or cmu sbcl) (setq alpha (complex-coerce alpha))
-
- (setf (aref *1x1-complex-array* 0) (realpart alpha))
- (setf (aref *1x1-complex-array* 1) (imagpart alpha))
- (zscal nxm *1x1-complex-array* (store result) 1)
-
- result))
-
-#+(or :cmu :sbcl)
-(defmethod scal ((alpha complex) (x complex-matrix))
- (scal (complex-coerce alpha) x))
-
+;;
+(generate-typed-scal!-func real-double-dscal!-typed double-float real-matrix-store-type real-matrix blas:dscal)
(defmethod scal! ((alpha number) (x number))
- (error "cannot SCAL! two scalars, arg X must
+ (error "Cannot SCAL! two scalars, arg X must
be a matrix to SCAL!"))
-(defmethod scal! ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
-
- (dscal nxm alpha (store x) 1)
- x))
+(defmethod scal! ((alpha complex) (x real-matrix))
+ (error "Cannot SCAL! a REAL-MATRIX by a COMPLEX, don't know
+how to coerce COMPLEX to REAL"))
(defmethod scal! ((alpha cl:real) (x real-matrix))
- (scal! (coerce alpha 'real-matrix-element-type) x))
+ (real-double-dscal!-typed (coerce alpha 'double-float) x))
-(defmethod scal! ((alpha complex) (x real-matrix))
- (error "cannot SCAL! a REAL-MATRIX by a COMPLEX, don't know
-how to coerce COMPLEX to REAL"))
+;;
+(generate-typed-scal!-func complex-double-dscal!-typed double-float complex-matrix-store-type complex-matrix blas:zdscal)
-(defmethod scal! ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
- (zdscal nxm alpha (store x) 1)
-
- x))
+(generate-typed-scal!-func complex-double-zscal!-typed (complex (double-float * *)) complex-matrix-store-type complex-matrix blas:zscal)
(defmethod scal! ((alpha cl:real) (x complex-matrix))
- (scal! (coerce alpha 'real-matrix-element-type) x))
+ (complex-double-dscal!-typed (coerce alpha 'double-float) x))
+
+(defmethod scal! ((alpha complex) (x complex-matrix))
+ (complex-double-zscal!-typed (complex-coerce alpha) x))
+
+;;;;
+(defgeneric scal (alpha x)
+ (:documentation
+"
+ Syntax
+ ======
+ (SCAL alpha x)
-(defmethod scal! ((alpha #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (x complex-matrix))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
+ Purpose
+ =======
+ Computes and returns a new matrix equal to
- #-(or cmu sbcl) (setq alpha (complex-coerce alpha))
+ alpha * X
- (setf (aref *1x1-complex-array* 0) (realpart alpha))
- (setf (aref *1x1-complex-array* 1) (imagpart alpha))
- (zscal nxm *1x1-complex-array* (store x) 1)
+ where alpha is a scalar and X is a matrix.
- x))
+"))
-#+(or :cmu :sbcl)
-(defmethod scal! ((alpha complex) (x complex-matrix))
- (scal! (complex-coerce alpha) x))
+(defmethod scal ((alpha number) (x number))
+ (* alpha x))
+
+;;
+(defmethod scal ((alpha cl:real) (x real-matrix))
+ (let ((result (copy x)))
+ (scal! alpha result)))
+(defmethod scal ((alpha complex) (x real-matrix))
+ (let* ((n (nrows x))
+ (m (ncols x))
+ (result (make-complex-matrix-dim n m)))
+ (declare (type fixnum n m))
+ (copy! x result)
+ (scal! alpha result)))
+;;
+(defmethod scal ((alpha number) (x complex-matrix))
+ (let ((result (copy x)))
+ (scal! alpha result)))
\ No newline at end of file
diff --git a/src/special.lisp b/src/special.lisp
index 93ff726..20138bc 100644
--- a/src/special.lisp
+++ b/src/special.lisp
@@ -153,7 +153,7 @@
(unity #.(coerce 1 'real-matrix-element-type)))
(declare (fixnum size)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (k size)
(declare (fixnum k))
(setf (aref store k) (random unity state)))
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index 142639e..47c4ff9 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -307,4 +307,27 @@ matrix and a number"))
:nrows nrows :ncols ncols
:store (store matrix)
:head (store-indexing i j hd rs cs)
- :row-stride rs :col-stride cs)))
\ No newline at end of file
+ :row-stride rs :col-stride cs)))
+
+;;
+(defmacro with-transpose! (matlst &rest body)
+ `(progn
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
+ ,@body
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
+
+;;
+(defun blas-copyable-p (matrix)
+ (declare (optimize (safety 0) (speed 3))
+ (type (or real-matrix complex-matrix) matrix))
+ (mlet* ((nr (nrows matrix) :type fixnum)
+ (nc (ncols matrix) :type fixnum)
+ (rs (row-stride matrix) :type fixnum)
+ (cs (col-stride matrix) :type fixnum)
+ (ne (number-of-elements matrix) :type fixnum))
+ (cond
+ ((= nc 1) (values t rs ne))
+ ((= nr 1) (values t cs ne))
+ ((= rs (* nc cs)) (values t cs ne))
+ ((= cs (* nr rs)) (values t rs ne))
+ (t (values nil -1 -1)))))
\ No newline at end of file
diff --git a/src/sum.lisp b/src/sum.lisp
index ad05d8f..49d0907 100644
--- a/src/sum.lisp
+++ b/src/sum.lisp
@@ -87,7 +87,7 @@
(nxm (number-of-elements a)))
(declare (type fixnum nxm)
(type real-matrix-element-type result)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (i nxm)
(declare (type fixnum i))
(incf result (aref store i)))
@@ -99,7 +99,7 @@
(store-a (store a))
(store-result (store result)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store-a store-result))
+ (type (real-matrix-store-type *) store-a store-result))
(dotimes (i n)
(declare (type fixnum i))
(setf (aref store-result i)
@@ -118,7 +118,7 @@
(nxm (number-of-elements a)))
(declare (type fixnum nxm)
(type complex-matrix-element-type imagpart realpart)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(dotimes (i nxm)
(declare (type fixnum i))
(let ((k (* 2 i)))
@@ -140,7 +140,7 @@
(store-a (store a))
(store-result (store result)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store-a store-result))
+ (type (complex-matrix-store-type *) store-a store-result))
(dotimes (i n)
(declare (type fixnum i))
(let ((realpart 0.0d0)
diff --git a/src/trans.lisp b/src/trans.lisp
index 48b6496..192d548 100644
--- a/src/trans.lisp
+++ b/src/trans.lisp
@@ -141,7 +141,7 @@
(new-store (allocate-real-store nxm)))
(declare (type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(dotimes (i n)
(declare (type fixnum i))
@@ -160,7 +160,7 @@
(new-store (allocate-complex-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(dotimes (i n)
(declare (type fixnum i))
@@ -201,7 +201,7 @@ a STANDARD-MATRIX, element types are not known"))
(new-store (allocate-complex-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(dotimes (i n)
(declare (type fixnum i))
diff --git a/src/utilities.lisp b/src/utilities.lisp
index 2dcbf64..76caf8c 100644
--- a/src/utilities.lisp
+++ b/src/utilities.lisp
@@ -1,6 +1,63 @@
(in-package :utilities)
;;
+(defmacro mlet* (decls &rest body)
+" mlet* ({ {(var*) | var} values-form &keyform declare type}*) form*
+
+ o var is just one symbol -> expands into let
+ o var is a list -> expands into multiple-value-bind
+
+ This macro also handles type declarations.
+
+ Example:
+ (mlet* ((x 2 :type fixnum :declare ((optimize (safety 0) (speed 3))))
+ ((a b) (floor 3) :type (nil fixnum)))
+ (+ x b))
+
+ expands into:
+
+ (let ((x 2))
+ (declare (optimize (safety 0) (speed 3))
+ (type fixnum x))
+ (multiple-value-bind (a b)
+ (floor 3)
+ (declare (ignore a)
+ (type fixnum b))
+ (+ x b)))
+"
+ (labels ((mlet-decl (vars type decls)
+ (when (or type decls)
+ `((declare ,@decls
+ ,@(when type
+ (mapcar #'(lambda (tv) (if (null (first tv))
+ `(ignore ,(second tv))
+ `(type ,(first tv) ,(second tv))))
+ (map 'list #'list type vars)))))))
+ ;;
+ (mlet-transform (elst nest-code)
+ (destructuring-bind (vars form &key declare type) elst
+ `(,(append (cond
+ ;;If there is only one element use let
+ ;;instead of multiple-value-bind
+ ((or (symbolp vars) (null (cdr vars)))
+ `(let ((,(car (ensure-list vars)) ,form))))
+ ;;
+ (t
+ `(multiple-value-bind (,@vars) ,form)))
+ (mlet-decl (ensure-list vars) (ensure-list type) declare)
+ nest-code))))
+ ;;
+ (mlet-walk (elst body)
+ (if (null elst)
+ `(,@body)
+ (mlet-transform (car elst) (mlet-walk (cdr elst) body)))))
+ ;;
+ (if decls
+ (car (mlet-walk decls body))
+ `(progn
+ ,@body))))
+
+;;
(defmacro with-gensyms (symlist &body body)
`(let ,(mapcar #'(lambda (sym)
`(,sym (gensym ,(symbol-name sym))))
@@ -55,64 +112,6 @@
lst
`(,lst)))
-(defmacro mlet* (decls &rest body)
-"
-mlet* ({ {(var*) | var} values-form &keyform declare type}*) form*
-
-o var is just one symbol -> expands into let
-o var is a list -> expands into multiple-value-bind
-
-This macro also handles type declarations.
-
-Example:
-(mlet* ((x 2 :type fixnum :declare ((optimize (safety 0) (speed 3))))
- ((a b) (floor 3) :type (nil fixnum)))
- (+ x b))
-
-expands into:
-
-(let ((x 2))
- (declare (optimize (safety 0) (speed 3))
- (type fixnum x))
- (multiple-value-bind (a b)
- (floor 3)
- (declare (ignore a)
- (type fixnum b))
- (+ x b)))
-"
- (labels ((mlet-decl (vars type decls)
- (when (or type decls)
- `((declare ,@decls
- ,@(when type
- (mapcar #'(lambda (tv) (if (null (first tv))
- `(ignore ,(second tv))
- `(type ,(first tv) ,(second tv))))
- (map 'list #'list type vars)))))))
- ;;
- (mlet-transform (elst nest-code)
- (destructuring-bind (vars form &key declare type) elst
- `(,(append (cond
- ;;If there is only one element use let
- ;;instead of multiple-value-bind
- ((or (symbolp vars) (null (cdr vars)))
- `(let ((,(car (ensure-list vars)) ,form))))
- ;;
- (t
- `(multiple-value-bind (,@vars) ,form)))
- (mlet-decl (ensure-list vars) (ensure-list type) declare)
- nest-code))))
- ;;
- (mlet-walk (elst body)
- (if (null elst)
- `(,@body)
- (mlet-transform (car elst) (mlet-walk (cdr elst) body)))))
- ;;
- (if decls
- (car (mlet-walk decls body))
- `(progn
- ,@body))))
-
-;;
(defmacro if-ret (form &rest else-body)
"if-ret (form &rest else-body)
Evaluate form, and if the form is not nil, then return it,
@@ -155,4 +154,6 @@ else run else-body"
,@(app-equal (cdr lst))))))
`(let ((,key-eval ,keyform))
(cond
- ,@(app-equal body))))))
\ No newline at end of file
+ ,@(app-equal body))))))
+
+;;
commit b0f1cb2dd42338c9189c83cbcbcb177eaf1c7845
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Mar 12 10:07:12 2012 +0530
Added slot-values to the list of symbols exported by utilities
diff --git a/packages.lisp b/packages.lisp
index f95b7ac..383af89 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -162,6 +162,7 @@
#:get-arg
#:nconsc
#:with-gensyms
+ #:slot-values
#:mlet*))
(defpackage :fortran-ffi-accessors
commit 98cdecd68f57b4a561ac8f68a0ede4a0374a6a95
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Mar 11 23:58:49 2012 +0530
-> Wrote macros to generate a typed copy! function.
Haven't yet adapted the generic methods yet.
diff --git a/src/blas.lisp b/src/blas.lisp
index f8adeed..999b389 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -151,9 +151,9 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
Y(0),Y(INCY), ... , Y((N-1)*INCY)
"
(n :integer :input)
- (dx (* :double-float :inc head-dx))
+ (dx (* :double-float :inc head-x))
(incx :integer :input)
- (dy (* :double-float :inc head-dy) :output)
+ (dy (* :double-float :inc head-y) :output)
(incy :integer :input)
)
@@ -328,9 +328,9 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
Y(0),Y(2*INCY), ... , Y(2*(N-1)*INCY)
"
(n :integer :input)
- (zx (* :complex-double-float :inc head-zx))
+ (zx (* :complex-double-float :inc head-x))
(incx :integer :input)
- (zy (* :complex-double-float :inc head-zy) :output)
+ (zy (* :complex-double-float :inc head-y) :output)
(incy :integer :input)
)
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 28785bc..48e7d0f 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -284,4 +284,21 @@
"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")))))
\ No newline at end of file
+ (error "require 1 or 2 arguments to make a matrix")))))
+
+;;
+
+(defun realpart! (mat)
+ (cond
+ ((typep mat 'real-matrix) mat)
+ ((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (* 2 (head mat))))))
+
+(defun imagpart! (mat)
+ (cond
+ ((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (+ 1 (* 2 (head mat)))))))
\ No newline at end of file
diff --git a/src/copy.lisp b/src/copy.lisp
index 264ce46..586aaf0 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -78,6 +78,7 @@
(in-package "MATLISP")
+;;
(defun blas-copyable-p (matrix)
(declare (optimize (safety 0) (speed 3))
(type (or real-matrix complex-matrix) matrix))
@@ -91,54 +92,45 @@
((= nr 1) (values t cs ne))
((= rs (* nc cs)) (values t cs ne))
((= cs (* nr rs)) (values t rs ne))
- (t nil))))
-
+ (t (values nil -1 -1)))))
+;;
+(defmacro with-transpose! (matlst &rest body)
+ `(progn
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
+ ,@body
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
-(defmacro generate-typed-copy!-func (func)
- `(defun ,func (matrix-a matrix-b)
- (declare (optimize (safety 0) (speed 3))
- (type (or ,matrix-type matrix-a matrix-b)))))
+;;
+(defmacro generate-typed-copy!-func (func store-type matrix-type blas-func)
+ `(defun ,func (mat-a mat-b)
+ (declare (type ,matrix-type mat-a mat-b)
+ (optimize (safety 0) (speed 3)))
+ (mlet* (((cp-a inc-a sz-a) (blas-copyable-p mat-a) :type (boolean fixnum nil))
+ ((cp-b inc-b sz-b) (blas-copyable-p mat-b) :type (boolean fixnum nil))
+ ((hd-a st-a sz) (slot-values mat-a '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum))
+ ((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *))))
+ (if (and cp-a cp-b)
+ (,blas-func sz (store mat-a) inc-a (store mat-b) inc-b :head-x hd-a :head-y hd-b)
+ (symbol-macrolet
+ ((common-code
+ (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum))
+ ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride))
+ :type (fixnum fixnum)))
+ (loop for i from 0 below nr-a
+ do (,blas-func nc-a st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
+ ;;Choose the smaller of the loops
+ (if (> (nrows mat-a) (ncols mat-a))
+ (with-transpose! (mat-a mat-b)
+ common-code)
+ common-code)))
+ mat-b)))
;;
(defvar *1x1-real-array* (make-array 1 :element-type 'double-float))
(defvar *1x1-complex-array* (make-array 2 :element-type 'double-float))
-(defgeneric copy (matrix)
- (:documentation
- "
- Syntax
- ======
- (COPY x)
-
- Purpose
- =======
- Return a copy of the matrix X"))
-
-(defmethod copy ((matrix standard-matrix))
- (make-instance 'standard-matrix :nrows (nrows matrix) :ncols (ncols matrix) :store (copy-seq (store matrix))))
-
-(defmethod copy ((matrix real-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
- (m (ncols matrix))
- (result (make-real-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:dcopy size (store matrix) 1 (store result) 1)
- result))
-
-(defmethod copy ((matrix complex-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
- (m (ncols matrix))
- (result (make-complex-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:zcopy size (store matrix) 1 (store result) 1)
- result))
-
-(defmethod copy ((matrix number))
- matrix)
-
;;
(defgeneric copy! (matrix new-matrix)
(:documentation
@@ -164,7 +156,6 @@
REAL-MATRIX but the converse is possible.
"))
-
(defmethod copy! :before ((x standard-matrix) (y standard-matrix))
(let ((nxm-x (number-of-elements x))
(nxm-y (number-of-elements y)))
@@ -172,6 +163,9 @@
(if (not (= nxm-x nxm-y))
(warn "arguments X,Y to COPY! are of different sizes"))))
+;;
+(generate-typed-copy!-func real-double-copy!-typed real-matrix-store-type real-matrix blas:dcopy)
+
(defmethod copy! ((x real-matrix) (y real-matrix))
(let* ((nxm-x (number-of-elements x))
(nxm-y (number-of-elements y))
@@ -197,6 +191,10 @@
(error "cannot copy a COMPLEX-MATRIX into a REAL-MATRIX,
don't know how to coerce a COMPLEX to a REAL"))
+
+;;
+(genera...
[truncated message content] |
