Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-108-g71aca48
|
From: Akshay S. <ak...@us...> - 2012-07-05 11:51:51
|
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 71aca48b041b5be2cd4c6ab8d514b260bdc02b19 (commit)
via 174d27300595c21a466a330fa34ab66fa7131bdf (commit)
via f25a68740987eeac4539d48b7a58d189da5b28e7 (commit)
via 536e528120dcf8631dbb9a8d4efd9af5541e55e0 (commit)
from d5f7ad309ca59d41c6e405c512f9a3544be01ea2 (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 71aca48b041b5be2cd4c6ab8d514b260bdc02b19
Author: Akshay Srinivasan <aks...@gm...>
Date: Thu Jul 5 17:16:14 2012 +0530
o Loopy has a new macro to recurse through lists of lists of lists ...
o Added a tensor-maker macro.
o dot now works.
diff --git a/README.org b/README.org
index 071d4fc..05bc347 100644
--- a/README.org
+++ b/README.org
@@ -26,6 +26,8 @@ This is the development branch of Matlisp.
* ODEPACK: Add abstraction for DLSODE, and DLSODAR.
* Tensor contraction: Hard to do very quickly.
Might have to copy stuff into a contiguous array; like Femlisp.
+*** Syntactic sugar
+ * Add array slicing macros
*** Python-bridge
(C)Python has far too many things, that we cannot even begin to hope to replicate.
diff --git a/matlisp.asd b/matlisp.asd
index 5db4d00..ec882aa 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -80,12 +80,10 @@
:components ((:file "blas")
(:file "lapack")
(:file "dfftpack")))
- (:module "matlisp-essentials"
+ (:module "matlisp-base"
:pathname "src/"
- :depends-on ("foreign-interface"
- "foreign-functions")
+ :depends-on ("foreign-functions")
:components ((:file "conditions")
- ;;
(:file "standard-tensor"
:depends-on ("conditions"))
;;
@@ -94,28 +92,26 @@
(:file "permutation"
:depends-on ("standard-tensor"))
(:file "blas-helpers"
- :depends-on ("standard-tensor" "permutation"))
- ;;
- (:file "real-tensor"
- :depends-on ("standard-tensor"))
- (:file "complex-tensor"
- :depends-on ("standard-tensor"))
- (:file "standard-matrix"
- :depends-on ("standard-tensor" "real-tensor" "complex-tensor"))
- ;; (:file "real-matrix"
- ;; :depends-on ("standard-matrix"))
- ;; (:file "complex-matrix"
- ;; :depends-on ("standard-matrix"))
+ :depends-on ("standard-tensor" "permutation"))
(:file "print"
- :depends-on ("standard-tensor" "standard-matrix"))
- ;;Copy, Scal
- (:file "copy"
- :depends-on ("real-tensor" "complex-tensor" "loopy"))
+ :depends-on ("standard-tensor"))))
+ (:module "matlisp-classes"
+ :pathname "src/"
+ :depends-on ("matlisp-base")
+ :components ((:file "real-tensor")
+ (:file "complex-tensor")
+ (:file "matrix"
+ :depends-on ("real-tensor" "complex-tensor"))))
+ (:module "matlisp-level-1"
+ :pathname "src/"
+ :depends-on ("matlisp-base" "matlisp-classes" "foreign-functions")
+ :components ((:file "tensor-maker")
+ (:file "copy")
+ (:file "dot")
(:file "scal"
- :depends-on ("copy" "loopy"))
+ :depends-on ("copy"))
(:file "realimag"
- :depends-on ("real-tensor" "complex-tensor" "copy"))
- ))))
+ :depends-on ("copy"))))))
;; (defclass f2cl-cl-source-file (asdf:cl-source-file)
diff --git a/packages.lisp b/packages.lisp
index 5226d68..d202ce8 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -162,7 +162,7 @@
#:recursive-append #:unquote-args #:flatten
#:format-to-string #:string+
#:linear-array-type
- #:seq-max #:seq-min
+ #:list-dimensions
;;Macros
#:when-let #:if-let #:if-ret #:with-gensyms #:let-rec
#:mlet* #:make-array-allocator
diff --git a/src/complex-tensor.lisp b/src/complex-tensor.lisp
index 4a636cf..aecac89 100644
--- a/src/complex-tensor.lisp
+++ b/src/complex-tensor.lisp
@@ -1,24 +1,24 @@
(in-package :matlisp)
-(eval-when (load eval compile)
- (deftype complex-base-type ()
- "The type of the elements stored in a COMPLEX-MATRIX"
- 'double-float)
-
- (deftype complex-base-array (size)
- "The type of the storage structure for a COMPLEX-MATRIX"
- `(simple-array complex-base-type (,size)))
-
- (deftype complex-type ()
- "Complex number with Re, Im parts in complex-base-type."
- '(cl:complex complex-base-type))
- )
+(deftype complex-base-type ()
+ "The type of the elements stored in a COMPLEX-MATRIX"
+ 'double-float)
+
+(deftype complex-base-array (size)
+ "The type of the storage structure for a COMPLEX-MATRIX"
+ `(simple-array complex-base-type (,size)))
+
+(deftype complex-type ()
+ "Complex number with Re, Im parts in complex-base-type."
+ '(cl:complex complex-base-type))
;;
(definline allocate-complex-store (size)
-"(allocate-complex-store size)
-Allocates real storage of size (* SIZE 2).
-Default initial-element = 0d0."
+ "
+ (allocate-complex-store size)
+ Allocates real storage of size (* SIZE 2).
+ Default initial-element = 0d0.
+"
(make-array (* 2 size) :element-type 'complex-base-type
:initial-element (coerce 0 'complex-base-type)))
@@ -87,11 +87,3 @@ Cannot hold complex numbers."))
"~11,5,,,,,'Eg"
"#C(~11,4,,,,,'Ee ~11,4,,,,,'Ee)")
realpart imagpart)))
-
-;;
-(defun make-complex-tensor-dims (&rest subs)
- (let* ((dims (make-index-store subs))
- (ss (reduce #'* dims))
- (store (allocate-complex-store ss)))
- (make-instance 'complex-tensor :store store :dimensions dims)))
-
diff --git a/src/conditions.lisp b/src/conditions.lisp
index 27986e0..842624e 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -33,7 +33,7 @@
;;Generic conditions---------------------------------------------;;
(defcondition generic-error (error)
((message :reader message :initarg :message :initform ""))
- (:method print-method ((c generic-error) stream)
+ (:method print-object ((c generic-error) stream)
(format stream (message c))))
(defcondition invalid-type (generic-error)
@@ -64,9 +64,25 @@
(to :reader to :initarg :to))
(:documentation "Cannot coerce one type into another.")
(:method print-object ((c coercion-error) stream)
- (format stream "Cannot coerce ~a into ~a." (from c) (to c))
+ (format stream "Cannot coerce ~a into ~a.~%" (from c) (to c))
(call-next-method)))
+(defcondition out-of-bounds-error (generic-error)
+ ((requested :reader requested :initarg :requested)
+ (bound :reader bound :initarg :bound))
+ (:documentation "General out-of-bounds error")
+ (:method print-object ((c out-of-bounds-error) stream)
+ (format stream "Out-of-bounds error, requested index : ~a, bound : ~a.~%" (requested c) (bound c))
+ (call-next-method)))
+
+(defcondition non-uniform-bounds-error (generic-error)
+ ((assumed :reader assumed :initarg :assumed)
+ (found :reader found :initarg :found))
+ (:documentation "Bounds are not uniform")
+ (:method print-object ((c out-of-bounds-error) stream)
+ (format stream "The bounds are not uniform, assumed bound : ~a, now found to be : ~a.~%" (assumed c) (found c))
+ (call-next-method)))
+
;;Tensor conditions----------------------------------------------;;
(define-condition tensor-error (error)
;;Optional argument for error-handling.
@@ -92,6 +108,12 @@
(:report (lambda (c stream)
(format stream "Given tensor with rank ~A, is not a matrix." (rank c)))))
+(define-condition tensor-not-vector (tensor-error)
+ ((tensor-rank :reader rank :initarg :rank))
+ (:documentation "Given tensor is not a vector.")
+ (:report (lambda (c stream)
+ (format stream "Given tensor with rank ~A, is not a vector." (rank c)))))
+
(define-condition tensor-index-out-of-bounds (tensor-error)
((argument :reader argument :initarg :argument)
(index :reader index :initarg :index)
diff --git a/src/dot.lisp b/src/dot.lisp
index 8983caa..1affebc 100644
--- a/src/dot.lisp
+++ b/src/dot.lisp
@@ -97,149 +97,67 @@
otherwise.
"))
-(defmethod dot ((x number) (y number) &optional (conjugate-p nil))
+(defmethod dot ((x number) (y number) &optional (conjugate-p t))
(if conjugate-p
(* (conjugate x) y)
(* x y)))
-(defmethod dot :before ((x standard-matrix) (y standard-matrix) &optional conjugate-p)
+(defmethod dot :before ((x standard-tensor) (y standard-tensor) &optional (conjugate-p t))
(declare (ignore conjugate-p))
- (if (not (row-or-col-vector-p x))
- (error "argument X to DOT is not a row or column vector")
- (if (not (row-or-col-vector-p y))
- (error "argument Y to DOT is not a row or column vector")
- (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))
- (error "arguments X,Y to DOT are not of the same size"))))))
-
-(defmethod dot ((x real-matrix) (y real-matrix) &optional conjugate-p)
+ (unless (and (vector-p x) (vector-p y))
+ (error 'tensor-not-vector
+ :rank (cond
+ ((not (vector-p x))
+ (rank x))
+ ((not (vector-p y))
+ (rank y)))))
+ (unless (idx= (dimensions x) (dimensions y))
+ (error 'tensor-dimension-mismatch)))
+
+(defmethod dot ((x real-tensor) (y real-tensor) &optional (conjugate-p t))
(declare (ignore conjugate-p))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
- (ddot nxm (store x) 1 (store y) 1)))
+ (ddot (number-of-elements x)
+ (store x) (aref (strides x) 0)
+ (store y) (aref (strides y) 0)
+ (head x) (head y)))
-;;#+(or :cmu :sbcl)
-(defmethod dot ((x real-matrix) (y complex-matrix) &optional conjugate-p)
+(defmethod dot ((x real-tensor) (y complex-tensor) &optional (conjugate-p t))
(declare (ignore conjugate-p))
- (let ((nxm (number-of-elements x))
- (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))
-
- (let ((realpart (ddot nxm store-x 1 store-y 2))
- (imagpart (with-vector-data-addresses ((addr-x store-x)
- (addr-y store-y))
- (incf-sap addr-y :double-float)
- (ddot nxm addr-x 1 addr-y 2))))
-
- (declare (type complex-matrix-element-type realpart imagpart))
-
- #+:complex-arg-implies-complex-result
- (complex realpart imagpart)
- #-:complex-arg-implies-comples-result
- (if (zerop imagpart)
- realpart
- (complex realpart imagpart))
- )))
-
-
-;;#+:allegro
-#+nil
-(defmethod dot ((x real-matrix) (y complex-matrix) &optional conjugate-p)
- (declare (ignore conjugate-p))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
-
- (let ((realpart 0.0d0)
- (imagpart 0.0d0))
- (declare (type complex-matrix-element-type realpart imagpart))
-
- (dotimes (i nxm)
- (declare (type fixnum i))
- (let ((x-elt (matrix-ref x i))
- (y-elt (matrix-ref y i)))
- (incf realpart (+ x-elt (realpart y-elt)))
- (incf imagpart (+ x-elt (imagpart y-elt)))))
-
-
- #+:complex-arg-implies-complex-result
- (complex realpart imagpart)
- #-:complex-arg-implies-comples-result
- (if (zerop imagpart)
- realpart
- (complex realpart imagpart))
- )))
-
-;;#+(or :cmu :sbcl)
-(defmethod dot ((x complex-matrix) (y real-matrix) &optional (conjugate-p t))
- (let ((nxm (number-of-elements x))
- (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))
-
- (let ((realpart (ddot nxm store-x 2 store-y 1))
- (imagpart (with-vector-data-addresses ((addr-x store-x)
- (addr-y store-y))
- (incf-sap addr-x :double-float)
- (ddot nxm addr-x 2 addr-y 1))))
-
- (declare (type complex-matrix-element-type realpart imagpart))
-
- (if conjugate-p
- (setq imagpart (- imagpart)))
-
- #+:complex-arg-implies-complex-result
- (complex realpart imagpart)
- #-:complex-arg-implies-comples-result
- (if (zerop imagpart)
- realpart
- (complex realpart imagpart))
- )))
-
-#+nil
-(defmethod dot ((x complex-matrix) (y real-matrix) &optional (conjugate-p t))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
-
- (let ((realpart 0.0d0)
- (imagpart 0.0d0))
- (declare (type complex-matrix-element-type realpart imagpart))
-
- (dotimes (i nxm)
- (declare (type fixnum i))
- (let ((x-elt (matrix-ref x i))
- (y-elt (matrix-ref y i)))
- (incf realpart (+ y-elt (realpart x-elt)))
- (incf imagpart (+ y-elt (imagpart x-elt)))))
-
- (if conjugate-p
- (setq imagpart (- imagpart)))
-
- #+:complex-arg-implies-complex-result
- (complex realpart imagpart)
- #-:complex-arg-implies-comples-result
- (if (zerop imagpart)
- realpart
- (complex realpart imagpart))
- )))
-
-(defmethod dot ((x complex-matrix) (y complex-matrix) &optional (conjugate-p t))
- (let* ((nxm (number-of-elements x))
- (store-x (store x))
- (store-y (store y))
- (dot (if conjugate-p
- (zdotc nxm store-x 1 store-y 1)
- (zdotu nxm store-x 1 store-y 1))))
-
- #-:complex-arg-implies-complex-result
- (if (zerop (imagpart dot))
- (realpart dot)
- dot)
- #+:complex-arg-implies-complex-result
- dot))
+ (let ((nele (number-of-elements x))
+ (std-x (aref (strides x) 0))
+ (hd-x (head x))
+ (std-y (aref (strides y) 0))
+ (hd-y (head y)))
+ (declare (type index-type nele std-x std-y hd-x hd-y))
+ (let ((rpart (ddot nele (store x) std-x (store y) (* 2 std-y) hd-x (* 2 hd-y)))
+ (ipart (ddot nele (store x) std-x (store y) (* 2 std-y) hd-x (1+ (* 2 hd-y)))))
+ (declare (type complex-base-type rpart ipart))
+ (if (zerop ipart)
+ rpart
+ (complex rpart ipart)))))
+
+(defmethod dot ((x complex-tensor) (y real-tensor) &optional (conjugate-p t))
+ (let ((cres (dot y x)))
+ (if conjugate-p
+ (conjugate cres)
+ cres)))
+
+(defmethod dot ((x complex-tensor) (y complex-tensor) &optional (conjugate-p t))
+ (let ((nele (number-of-elements x))
+ (std-x (aref (strides x) 0))
+ (hd-x (head x))
+ (std-y (aref (strides y) 0))
+ (hd-y (head y)))
+ (declare (type index-type nele std-x hd-x std-y hd-y))
+ (let ((ret (if conjugate-p
+ (zdotc nele
+ (store x) std-x
+ (store y) std-y
+ hd-x hd-y)
+ (zdotu nele
+ (store x) std-x
+ (store y) std-y
+ hd-x hd-y))))
+ (if (zerop (imagpart ret))
+ (realpart ret)
+ ret))))
diff --git a/src/loopy.lisp b/src/loopy.lisp
index 32037bc..1828f66 100644
--- a/src/loopy.lisp
+++ b/src/loopy.lisp
@@ -8,38 +8,28 @@
{with (loop-order {:row-major :col-major})}
{with (linear-sums
{(offsets {stride-seq})}*)}
- {with (variables
- {(vars init &key type)}*)}
{do ({code}*)})
Examples:
- > (mod-dotimes (idx (vidx 2 2))
- with (linear-sums (of (vidx 2 1)))
+ > (mod-dotimes (idx (idxv 2 2))
+ with (linear-sums (of (idxv 2 1)))
do (format t \"~a ~a~%\" idx of))
#(0 0) 0
#(0 1) 1
#(1 0) 2
#(1 1) 3
- > (mod-dotimes (idx (vidx 2 2))
+ > (mod-dotimes (idx (idxv 2 2))
with (loop-order :col-major)
- with (linear-sums (of (vidx 2 1)))
+ with (linear-sums (of (idxv 2 1)))
do (format t \"~a ~a~%\" idx of))
#(0 0) 0
#(1 0) 2
#(0 1) 1
#(1 1) 3
- > (mod-dotimes (idx (vidx 2 2))
- with (variables (tmp 1d0 :type double-float))
- with (linear-sums (of (vidx 2 1)))
- do (progn
- (format t \"~a ~a ~a~%\" idx of tmp)
- (incf tmp)))
- #(0 0) 0 0d0
- #(0 1) 1 1d0
- #(1 0) 2 2d0
- #(1 1) 3 3d0
+ Make sure that \"do\" is specified at the end. Parser stops
+ at the first 'do it finds.
"
(check-type idx symbol)
(labels ((parse-code (body ret)
@@ -50,6 +40,11 @@
(multiple-value-bind (indic decl) (parse-with (cadr body))
(setf (getf ret indic) decl))
(parse-code (cddr body) ret))
+ ;;Let's not do too much
+ #+nil
+ ((eq (car body) 'finally)
+ (setf (getf ret :finally) (second body))
+ (parse-code (cddr body) ret))
((eq (car body) 'do)
(values (cadr body) ret))
(t (error 'unknown-token :token (car body) :message "Error in macro: mod-dotimes -> parse-code.~%"))))
@@ -66,6 +61,8 @@
((and (eq (car code) 'loop-order)
(member (cadr code) '(:row-major :col-major)))
(values :loop-order (second code)))
+ ;;Useless without a finally clause
+ #+nil
((eq (car code) 'variables)
(values :variables
(loop for decl in (cdr code)
@@ -119,4 +116,111 @@
(unless (= ,cstrd 0)
(incf ,(getf decl :offset-sym) ,cstrd)))))
(return t)))
- finally (return nil))))))))))))
+ finally (return nil))))
+ ,@(unless (null (getf sdecl :finally))
+ `(finally (,@(getf sdecl :finally))))))))))))
+
+(defmacro list-loop ((idx ele lst) &rest body)
+ "
+ (list-loop (idx ele {list}) compound-form*)
+
+ Examples:
+ > (list-loop (idx ele '((1 2) (4 5)))
+ with (linear-sums (of (idxv 2 1)))
+ do (format t \"~a ~a ~a~%\" idx of ele))
+ #(0 0) 0 1
+ #(0 1) 1 2
+ #(1 0) 2 4
+ #(1 1) 3 5
+"
+ (check-type idx symbol)
+ (check-type ele symbol)
+ (labels ((parse-code (body ret)
+ (cond
+ ((null body)
+ (values nil ret))
+ ((eq (car body) 'with)
+ (multiple-value-bind (indic decl) (parse-with (cadr body))
+ (setf (getf ret indic) decl))
+ (parse-code (cddr body) ret))
+ ;;Let's not do too much.
+ #+nil
+ ((eq (car body) 'finally)
+ (setf (getf ret :finally) (second body))
+ (parse-code (cddr body) ret))
+ ((eq (car body) 'do)
+ (values (cadr body) ret))
+ (t (error 'unknown-token :token (car body) :message "Error in macro: mod-dotimes -> parse-code.~%"))))
+ (parse-with (code)
+ (cond
+ ((eq (car code) 'linear-sums)
+ (values :linear-sums
+ (loop for decl in (cdr code)
+ collect (destructuring-bind (offst strds &optional (init 0)) decl
+ (list :offset-sym offst
+ :offset-init init
+ :stride-sym (gensym (string+ (symbol-name offst) "-stride"))
+ :stride-expr strds)))))
+ ;;Traversing the list the other way is far too inefficient and/or too hard to do.
+ #+nil
+ ((and (eq (car code) 'loop-order)
+ (member (cadr code) '(:row-major :col-major)))
+ (values :loop-order (second code)))
+ ;;Useless without a finally clause.
+ #+nil
+ ((eq (car code) 'variables)
+ (values :variables
+ (loop for decl in (cdr code)
+ collect (destructuring-bind (sym init &key type) decl
+ (list :variable sym
+ :init init
+ :type type)))))
+ (t (error 'unknown-token :token (car code) :message "Error in macro: mod-dotimes -> parse-with.~%")))))
+ (multiple-value-bind (code sdecl) (parse-code body nil)
+ (with-gensyms (lst-sym dims-sym rank-sym lst-rec-sym lst-rec-count-sym lst-rec-lst-sym)
+ `(let ((,lst-sym ,lst))
+ (declare (type list ,lst-sym))
+ (let ((,dims-sym (make-index-store (list-dimensions ,lst-sym))))
+ (declare (type (index-array *) ,dims-sym))
+ (let ((,rank-sym (array-dimension ,dims-sym 0)))
+ (declare (type index-type ,rank-sym))
+ (let ((,idx (allocate-index-store ,rank-sym))
+ ,@(mapcar #'(lambda (x) `(,(getf x :offset-sym) ,(getf x :offset-init))) (getf sdecl :linear-sums))
+ ,@(mapcar #'(lambda (x) `(,(getf x :stride-sym) ,(getf x :stride-expr))) (getf sdecl :linear-sums))
+ ,@(mapcar #'(lambda (x) `(,(getf x :variable) ,(getf x :init))) (getf sdecl :variables)))
+ (declare (type (index-array *) ,idx)
+ ,@(when (getf sdecl :linear-sums)
+ `((type (index-array *) ,@(mapcar #'(lambda (x) (getf x :stride-sym)) (getf sdecl :linear-sums)))
+ (type index-type ,@(mapcar #'(lambda (x) (getf x :offset-sym)) (getf sdecl :linear-sums)))))
+ ,@(loop for x in (getf sdecl :variables)
+ unless (null (getf x :type))
+ collect `(type ,(getf x :type) ,(getf x :variable))))
+ (labels ((,lst-rec-sym (,lst-rec-count-sym ,lst-rec-lst-sym)
+ (if (null ,lst-rec-lst-sym)
+ (progn
+ (unless (= (aref ,idx ,lst-rec-count-sym) (aref ,dims-sym ,lst-rec-count-sym))
+ (error 'non-uniform-bounds-error :assumed (aref ,dims-sym ,lst-rec-count-sym) :found ,lst-rec-count-sym
+ :message "Error in list-loop, given list is not uniform in dimensions."))
+ (setf (aref ,idx ,lst-rec-count-sym) 0)
+ ,@(loop
+ for decl in (getf sdecl :linear-sums)
+ collect `(decf ,(getf decl :offset-sym) (* (aref ,(getf decl :stride-sym) ,lst-rec-count-sym) (aref ,dims-sym ,lst-rec-count-sym))))
+ ,@(if (null (getf sdecl :finally))`(nil)
+ `((when (= ,lst-rec-count-sym 0)
+ ,(getf sdecl :finally)))))
+ (progn
+ ;;list-dimensions does not parse the entire list, just goes through caaa..r's to find out the
+ ;;dimensions if it is uniform.
+ (unless (< -1 (aref ,idx ,lst-rec-count-sym) (aref ,dims-sym ,lst-rec-count-sym))
+ (error 'out-of-bounds-error :requested ,lst-rec-count-sym :bound (aref ,dims-sym ,lst-rec-count-sym)
+ :message "Error in list-loop, given list is not uniform in dimensions."))
+ (if (consp (car ,lst-rec-lst-sym))
+ (,lst-rec-sym (1+ ,lst-rec-count-sym) (car ,lst-rec-lst-sym))
+ (let ((,ele (car ,lst-rec-lst-sym)))
+ ,code))
+ (incf (aref ,idx ,lst-rec-count-sym))
+ ,@(loop
+ for decl in (getf sdecl :linear-sums)
+ collect `(incf ,(getf decl :offset-sym) (aref ,(getf decl :stride-sym) ,lst-rec-count-sym)))
+ (,lst-rec-sym ,lst-rec-count-sym (cdr ,lst-rec-lst-sym))))))
+ (,lst-rec-sym 0 ,lst-sym))))))))))
diff --git a/src/matrix.lisp b/src/matrix.lisp
index 4454c6f..e8fa449 100644
--- a/src/matrix.lisp
+++ b/src/matrix.lisp
@@ -1,511 +1,87 @@
-;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Copyright (c) 2000 The Regents of the University of California.
-;;; All rights reserved.
-;;;
-;;; Permission is hereby granted, without written agreement and without
-;;; license or royalty fees, to use, copy, modify, and distribute this
-;;; software and its documentation for any purpose, provided that the
-;;; above copyright notice and the following two paragraphs appear in all
-;;; copies of this software.
-;;;
-;;; IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
-;;; FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
-;;; ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
-;;; THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
-;;; SUCH DAMAGE.
-;;;
-;;; THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
-;;; INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-;;; MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
-;;; PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
-;;; CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
-;;; ENHANCEMENTS, OR MODIFICATIONS.
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Originally written by Raymond Toy
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; $Id: matrix.lisp,v 1.16 2011/01/25 18:36:56 rtoy Exp $
-;;;
-;;; $Log: matrix.lisp,v $
-;;; Revision 1.16 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.15.2.1 2011/01/25 18:16:53 rtoy
-;;; Use cl:real instead of real.
-;;;
-;;; Revision 1.15 2010/12/12 02:07:31 rtoy
-;;; matrix.lisp:
-;;;
-;;; o Apply patch from Nicolas Neuss for matrices with 0 dimensions. (See
-;;; <http://sourceforge.net/mailarchive/message.php?msg_id=1124576>)
-;;;
-;;; print.lisp:
-;;; o Apparently the above patch to print was also applied previously. We
-;;; just fix a bug in printing 0xm matrices. Just exit early if the
-;;; matrix has no dimensions.
-;;;
-;;; Revision 1.14 2004/05/24 16:34:22 rtoy
-;;; More SBCL support from Robert Sedgewick. The previous SBCL support
-;;; was incomplete.
-;;;
-;;; Revision 1.13 2003/05/31 22:20:26 rtoy
-;;; o Add some support for CMUCL with Gerd's PCL so we can inline
-;;; accessors and such for the matrix classes.
-;;; o Only use one, system-independent, standard-matrix class.
-;;; o Try to declare the types of the slots of the matrix classes
-;;; o FORTRAN-MATRIX-INDEXING changed to use fixnum arithmetic.
-;;;
-;;; Revision 1.12 2003/03/09 14:26:30 rtoy
-;;; Forgot one more :type 'fixnum bug. From Gerd Moellmann.
-;;;
-;;; Revision 1.11 2003/02/19 21:59:52 rtoy
-;;; Correct the slot type declarations.
-;;;
-;;; Revision 1.10 2001/10/29 18:00:28 rtoy
-;;; Updates from M. Koerber to support QR routines with column pivoting:
-;;;
-;;; o Add an integer4 type and allocate-integer4-store routine.
-;;; o Add the necessary Fortran routines
-;;; o Add Lisp interface to the Fortran routines
-;;; o Update geqr for the new routines.
-;;;
-;;; Revision 1.9 2001/06/22 12:51:49 rtoy
-;;; o Added ALLOCATE-REAL-STORE and ALLOCATE-COMPLEX-STORE functions to
-;;; allocate appropriately sized arrays for holding real and complex
-;;; matrix elements.
-;;; o Use it to allocate space.
-;;;
-;;; Revision 1.8 2000/10/04 15:56:50 simsek
-;;; o Fixed bug in (MAKE-COMPLEX-MATRIX n)
-;;;
-;;; Revision 1.7 2000/07/11 18:02:03 simsek
-;;; o Added credits
-;;;
-;;; Revision 1.6 2000/07/11 02:11:56 simsek
-;;; o Added support for Allegro CL
-;;;
-;;; Revision 1.5 2000/05/11 18:28:10 rtoy
-;;; After the great standard-matrix renaming, row-vector-p and
-;;; col-vector-p were swapped.
-;;;
-;;; Revision 1.4 2000/05/11 18:02:55 rtoy
-;;; o After the great standard-matrix renaming, I forgot a few initargs
-;;; that needed to be changed
-;;; o MAKE-REAL-MATRIX and MAKE-COMPLEX-MATRIX didn't properly handle
-;;; things like #(1 2 3 4) and #((1 2 3 4)). Make them accept these.
-;;;
-;;; Revision 1.3 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.2 2000/05/05 21:35:16 simsek
-;;; o Fixed row-vector-p and col-vector-p
-;;;
-;;; 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)
-;;; Definitions of STANDARD-MATRIX, REAL-MATRIX, COMPLEX-MATRIX.
-
-(in-package "MATLISP")
-
-#+nil (export '(real-matrix
- complex-matrix
- standard-matrix
- real-matrix-element-type
- real-matrix-store-type
- complex-matrix-element-type
- complex-matrix-store-type
- #|
- n
- m
- nxm
- |#
- nrows
- ncols
- number-of-elements
- row-vector-p
- col-vector-p
- row-or-col-vector-p
- square-matrix-p
- size
- fortran-matrix-indexing
- fortran-complex-matrix-indexing
- complex-coerce
- fill-matrix
- make-real-matrix-dim
- make-real-matrix
- make-complex-matrix-dim
- make-complex-matrix))
-
-(eval-when (load eval compile)
-(deftype integer4-matrix-element-type ()
- '(signed-byte 32))
-
-(deftype real-matrix-element-type ()
- "The type of the elements stored in a REAL-MATRIX"
- 'double-float)
-
-(deftype real-matrix-store-type (size)
- "The type of the storage structure for a REAL-MATRIX"
- `(simple-array double-float ,size))
-
-(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))
-)
-
-(declaim (ftype (function (standard-matrix) fixnum)
- n
- m
- nxm
- store-size)
- (ftype (function (real-matrix) fixnum)
- n
- m
- nxm
- store-size)
- (ftype (function (complex-matrix) fixnum)
- n
- m
- nxm
- store-size)
- (ftype (function (real-matrix) (simple-array double-float (*)))
- store)
- (ftype (function (complex-matrix) (simple-array double-float (*)))
- store))
-
-
-#|
-(defgeneric n (matrix)
- (:documentation
-"
- Syntax
- ======
- (N matrix)
-
- Purpose
- =======
- Returns the number of rows of MATRIX.
-"))
-
-(defgeneric m (matrix)
- (:documentation
-"
- Syntax
- ======
- (M matrix)
-
- Purpose
- =======
- Returns the number of columns of MATRIX.
-"))
-
-(defgeneric nxm (matrix)
- (:documentation
-"
- Syntax
- ======
- (NxM matrix)
-
- Purpose
- =======
- Returns the number of elements of MATRIX;
- which is number of rows * number of columns.
-"))
-|#
-
-(defgeneric store-size (matrix)
- (:documentation
-"
- Syntax
- ======
- (STORE-SIZE matrix)
-
- Purpose
- =======
- Total number of elements needed to store the matrix. (Usually
- the same as (NxM matrix), but not necessarily so!
-"))
-
-(defgeneric store (matrix)
- (:documentation
-"
- Syntax
- ======
- (STORE matrix)
-
- Purpose
- =======
-The actual storage for the matrix. It is typically a one dimensional
-array but not necessarily so. The float and complex matrices do use
-1-D arrays. The complex matrix actually stores the real and imaginary
-parts in successive elements of the matrix because Fortran stores them
-that way.
-"))
+;;
+(defclass standard-matrix (standard-tensor)
+ ((rank
+ :accessor rank
+ :type index-type
+ :initform 2
+ :documentation "For a matrix, rank = 2."))
+ (:documentation "Basic matrix class."))
-#+(and (or cmu sbcl) gerds-pcl)
-(declaim (ext:slots (slot-boundp real-matrix complex-matrix)
- (inline standard-matrix real-matrix complex-matrix)))
+(defmethod print-object ((tensor standard-matrix) stream)
+ (print-unreadable-object (tensor stream :type t)
+ (format stream "~A x ~A~%" (nrows tensor) (ncols tensor))
+ (print-tensor tensor stream)))
-(defclass standard-matrix ()
- ((number-of-rows
- :initarg :nrows
- :initform 0
- :accessor nrows
- :type fixnum
- :documentation "Number of rows in the matrix")
- (number-of-cols
- :initarg :ncols
- :initform 0
- :accessor ncols
- :type fixnum
- :documentation "Number of columns in the matrix")
- (number-of-elements
- :initarg :nels
- :initform 0
- :accessor number-of-elements
- :type fixnum
- :documentation "Total number of elements in the matrix (nrows * ncols)")
- (store-size
- :initarg :store-size
- :initform 0
- :accessor store-size
- :type fixnum
- :documentation "Total number of elements needed to store the matrix. (Usually
-the same as nels, but not necessarily so!")
- (store
- :initarg :store
- :accessor store
- :documentation "The actual storage for the matrix. It is typically a one dimensional
-array but not necessarily so. The float and complex matrices do use
-1-D arrays. The complex matrix actually stores the real and imaginary
-parts in successive elements of the matrix because Fortran stores them
-that way."))
- (:documentation "Basic matrix class."))
+(definline nrows (matrix)
+ (declare (type standard-matrix matrix))
+ (aref (dimensions matrix) 0))
+(definline ncols (matrix)
+ (declare (type standard-matrix matrix))
+ (aref (dimensions matrix) 1))
-#+(and nil :allegro)
-(defclass standard-matrix ()
- ((number-of-rows
- :initarg :nrows
- :initform 0
- :accessor nrows
- :documentation "Number of rows in the matrix")
- (number-of-cols
- :initarg :ncols
- :initform 0
- :accessor ncols
- :documentation "Number of columns in the matrix")
- (number-of-elements
- :initarg :nels
- :initform 0
- :accessor number-of-elements
- :documentation "Total number of elements in the matrix (nrows * ncols)")
- (store-size
- :initarg :store-size
- :initform 0
- :accessor store-size
- :documentation "Total number of elements needed to store the matrix. (Usually
-the same as nels, but not necessarily so!")
- (store
- :initarg :store
- :accessor store
- :documentation "The actual storage for the matrix. It is typically a one dimensional
-array but not necessarily so. The float and complex matrices do use
-1-D arrays. The complex matrix actually stores the real and imaginary
-parts in successive elements of the matrix because Fortran stores them
-that way."))
- (:documentation "Basic matrix class."))
+(definline row-stride (matrix)
+ (declare (type standard-matrix matrix))
+ (aref (strides matrix) 0))
-(defclass real-matrix (standard-matrix)
- ((store
- :type (simple-array real-matrix-element-type (*))))
- (:documentation "A class of matrices with real elements."))
+(definline col-stride (matrix)
+ (declare (type standard-matrix matrix))
+ (aref (strides matrix) 1))
-(defclass complex-matrix (standard-matrix)
- ((store
- :type (simple-array complex-matrix-element-type (*))))
- (:documentation "A class of matrices with complex elements."))
+(definline size (matrix)
+ (declare (type standard-matrix matrix))
+ (let ((dims (dimensions matrix)))
+ (declare (type (index-array 2) dims))
+ (list (aref dims 0) (aref dims 1))))
+;;
(defmethod initialize-instance :after ((matrix standard-matrix) &rest initargs)
(declare (ignore initargs))
- (let* ((n (nrows matrix))
- (m (ncols matrix))
- (nxm (* n m)))
- (declare (type fixnum n m nxm))
- (setf (number-of-elements matrix) nxm)
- (setf (store-size matrix) nxm)))
-
-(defmethod make-load-form ((matrix standard-matrix) &optional env)
- "MAKE-LOAD-FORM allows us to determine a load time value for
- matrices, for example #.(make-matrix ...)"
- (make-load-form-saving-slots matrix :environment env))
+ (mlet*
+ ((rank (rank matrix) :type index-type))
+ (unless (= rank 2)
+ (error 'tensor-not-matrix :rank rank :tensor matrix))))
-(defgeneric row-vector-p (matrix)
- (:documentation "
+;;
+(definline row-matrix-p (matrix)
+ "
Syntax
======
- (ROW-VECTOR-P x)
+ (ROW-MATRIX-P x)
Purpose
=======
- Return T if X is a row vector (number of columns is 1)"))
+ Return T if X is a row matrix (number of columns is 1)"
+ (tensor-type-p matrix '(1 *)))
-(defgeneric col-vector-p (matrix)
- (:documentation "
+(definline col-matrix-p (matrix)
+ "
Syntax
======
- (COL-VECTOR-P x)
+ (COL-MATRIX-P x)
Purpose
=======
- Return T if X is a column vector (number of rows is 1)"))
+ Return T if X is a column matrix (number of rows is 1)"
+ (tensor-type-p matrix '(* 1)))
-(defgeneric row-or-col-vector-p (matrix)
- (:documentation "
- Syntax
- ======
- (ROW-OR-COL-VECTOR-P x)
-
- Purpose
- =======
- Return T if X is either a row or a column vector"))
-
-(defgeneric square-matrix-p (matrix)
- (:documentation "
- Syntax
- ======
- (SQUARE-MATRIX-P x)
-
- Purpose
- =======
- Return T if X is square matrix"))
-
-(defgeneric size (matrix)
- (:documentation "
+(definline row-or-col-matrix-p (matrix)
+"
Syntax
======
- (SIZE x)
+ (ROW-OR-COL-matrix-P x)
Purpose
=======
- Return the number of rows and columns of the matrix X as a list"))
-
-(declaim (inline row-vector-p))
-(defmethod row-vector-p ((matrix standard-matrix))
- (= (nrows matrix) 1))
-
-(declaim (inline col-vector-p))
-(defmethod col-vector-p ((matrix standard-matrix))
- (= (ncols matrix) 1))
-
-(declaim (inline row-or-col-vector-p))
-(defmethod row-or-col-vector-p ((matrix standard-matrix))
+ Return T if X is either a row or a column matrix."
(or (row-vector-p matrix) (col-vector-p matrix)))
-(declaim (inline square-matrix-p))
-(defmethod square-matrix-p ((matrix standard-matrix))
- (= (nrows matrix) (ncols matrix)))
+(defun square-matrix-p (matrix)
+ (and (square-p matrix) (matrix-p matrix)))
-(defmethod size ((matrix standard-matrix))
- (list (nrows matrix) (ncols matrix)))
-
-;; For compatibility with Fortran, matrices are stored in column major
-;; order instead of row major order. Also, we store the matrix as a
-;; one-dimensional array instead of a two-dimensional array. This
-;; makes it easy to interface to LAPACK routines.
;;
-;; furthermore, this next macro should really be left as a macro
-;; to avoid integer to pointer coercions, since FORTRAN-MATRIX-INDEXING
-;; will be called too many times.
-
-#+nil
-(defmacro fortran-matrix-indexing (i j l)
- `(let ((i ,i)
- (j ,j)
- (l ,l))
- (declare (optimize (speed 3) (safety 0))
- (type fixnum i j l))
- (let* ((q (* j l))
- (p (+ i q)))
- (declare (type fixnum q p))
- p)))
-
-(declaim (inline fortran-matrix-indexing))
-(defun fortran-matrix-indexing (row col nrows)
- (declare (type (and fixnum (integer 0)) row col nrows))
- (the fixnum (+ row (the fixnum (* col nrows)))))
-
-;; For matrices with complex-valued elements, we store the array as a
-;; double-length double-precision floating-point vector, as Fortran
-;; does too. The first element is the real part; the second, the
-;; imaginary part.
-
-#+nil
-(defmacro fortran-complex-matrix-indexing (i j l)
- `(let ((i ,i)
- (j ,j)
- (l ,l))
- (declare (optimize (speed 3) (safety 0))
- (type fixnum i j l))
- (let* ((q (* j l))
- (p (+ i q))
- (r (* 2 p)))
- (declare (type fixnum q p r))
- r)))
-
-(declaim (inline fortran-complex-matrix-indexing))
-(defun fortran-complex-matrix-indexing (row col nrows)
- (declare (type (and fixnum (integer 0)) row col nrows))
- (the fixnum (* 2 (the fixnum (+ row (the fixnum (* col nrows)))))))
-
-
-
-;;; coerce is broken in CMUCL. Here is a function
-;;; that implements coerce correctly for what we want.
-
-(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))))
-
(defgeneric fill-matrix (matrix fill-element)
(:documentation
"
@@ -518,390 +94,43 @@ that way."))
Fill MATRIX with FILL-ELEMENT.
"))
-(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"))
-
-(defmethod fill-matrix ((matrix complex-matrix) (fill number))
- (copy! fill matrix))
-
(defmethod fill-matrix ((matrix t) (fill t))
(error "arguments MATRIX and FILL to FILL-MATRIX must be a
matrix and a number"))
-;; Allocate an array suitable for the store part of a real matrix.
-
-(declaim (inline allocate-integer4-store))
-(defun allocate-integer4-store (size &optional (initial-element 0))
- "(ALLOCATE-INTEGER-STORE SIZE [INITIAL-ELEMENT]). Allocates
-integer storage. Default INITIAL-ELEMENT = 0."
- (make-array size
- :element-type 'integer4-matrix-element-type
- :initial-element initial-element))
-
-(declaim (inline allocate-real-store))
-(defun allocate-real-store (size &optional (initial-element 0))
- (make-array size :element-type 'real-matrix-element-type
- :initial-element (coerce initial-element 'real-matrix-element-type)))
-
-(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)))
-
-(defun make-real-matrix-dim (n m &optional (fill 0.0d0))
- "
- 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))
- (make-instance 'real-matrix :nrows n :ncols m
- :store (allocate-real-store (* n m) casted-fill))))
-
-
-;;; Make a matrix from a 2-D Lisp array
-(defun make-real-matrix-array (array)
- "
- 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))
- (dotimes (i n)
- (declare (type fixnum i))
- (dotimes (j m)
- (declare (type fixnum j))
- (setf (aref store (fortran-matrix-indexing i j n))
- (coerce (aref array i j) 'real-matrix-element-type))))
- (make-instance 'real-matrix :nrows n :ncols m :store store)))
-
-(defun make-real-matrix-seq-of-seq (seq)
- (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))
- (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 (fortran-matrix-indexing i j n))
- (coerce (elt this-row j) 'real-matrix-element-type)))))
- (make-instance 'real-matrix :nrows n :ncols m :store store)))
-
-(defun make-real-matrix-seq (seq)
- (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)))
- (make-instance 'real-matrix :nrows n :ncols 1 :store store)))
-
-(defun make-real-matrix-sequence (seq)
- (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))
- (t
- ;; Assume a simple sequence
- (make-real-matrix-seq seq)))))
- ((arrayp seq)
- (make-real-matrix-array seq))))
-
-(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")))))
-
-
-
-(defun make-complex-matrix-dim (n m &optional (fill #c(0.0d0 0.0d0)))
- "
- Syntax
- ======
- (MAKE-COMPLEX-MATRIX-DIM n m [fill-element])
-
- Purpose
- =======
- Creates an NxM COMPLEX-MATRIX with initial contents FILL-ELEMENT,
- the default #c(0.0d0 0.0d0)
-
- See MAKE-COMPLEX-MATRIX.
-"
- (declare (type fixnum n m))
- (let* ((size (* n m))
- (store (allocate-complex-store size))
- (matrix (make-instance 'complex-matrix :nrows n :ncols m :store store)))
-
- (fill-matrix matrix fill)
- matrix))
-
-(defun make-complex-matrix-array (array)
- "
- Syntax
- ======
- (MAKE-COMPLEX-MATRIX-ARRAY array)
-
- Purpose
- =======
- Creates a COMPLEX-MATRIX with the same contents as ARRAY.
-"
- (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))
- (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 (fortran-complex-matrix-indexing i j n)))
- (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 :store store)))
-
-
-(defun make-complex-matrix-seq-of-seq (seq)
- (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))
-
- (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 (fortran-complex-matrix-indexing i j n)))
- (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 :store store)))
-
-
-(defun make-complex-matrix-seq (seq)
- (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)))
-
- (make-instance 'complex-matrix :nrows n :ncols 1 :store store)))
-
-
-(defun make-complex-matrix-sequence (seq)
- (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))
- (t
- ;; Assume a simple sequence
- (make-complex-matrix-seq seq)))))
- ((arrayp seq)
- (make-complex-matrix-array seq))))
-
-
-(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:
+;;
+(defclass real-matrix (standard-matrix real-tensor)
+ ()
+ (:documentation "A class of matrices with real elements."))
- 1 2 3
- 4 5 6
+(defclass real-sub-matrix (real-matrix standard-sub-tensor)
+ ()
+ (:documentation "Sub-matrix class with real elements."))
- (make-complex-matrix #((1 2 3) (4 5 6)))
- 2x3 matrix:
+(setf (gethash 'real-matrix *sub-tensor-counterclass*) 'real-sub-matrix
+ (gethash 'real-sub-matrix *sub-tensor-counterclass*) 'real-sub-matrix
+ ;;
+ (gethash 'real-matrix *tensor-class-optimizations*) 'real-tensor
+ (gethash 'real-sub-matrix *tensor-class-optimizations*) 'real-tensor)
+;;
- 1 2 3
- 4 5 6
+(defclass complex-matrix (standard-matrix complex-tensor)
+ ()
+ (:documentation "A class of matrices with complex elements."))
- (make-complex-matrix #((1 2 3) #(4 5 6)))
- 2x3 matrix:
+(defclass complex-sub-matrix (complex-matrix standard-sub-tensor)
+ ()
+ (:documentation "Sub-matrix class with complex elements."))
- 1 2 3
- 4 5 6
+(setf (gethash 'complex-matrix *sub-tensor-counterclass*) 'complex-sub-matrix
+ (gethash 'complex-sub-matrix *sub-tensor-counterclass*) 'complex-sub-matrix
+ ;;
+ (gethash 'complex-matrix *tensor-class-optimizations*) 'complex-tensor
+ (gethash 'complex-sub-matrix *tensor-class-optimizations*) 'complex-tensor)
- (make-complex-matrix #2a((1 2 3) (4 5 6)))
- 2x3 matrix:
+;;
- 1 2 3
- 4 5 6
+(definline matrix-ref (matrix row &optional col)
+ (declare (type standard-matrix matrix))
+ (tensor-ref matrix `(,row ,col)))
-"
- (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/print.lisp b/src/print.lisp
index 6b555a7..c05478b 100644
--- a/src/print.lisp
+++ b/src/print.lisp
@@ -163,8 +163,3 @@ of a matrix (default 0)
(print-unreadable-object (tensor stream :type t)
(format stream "~A~%" (dimensions tensor))
(print-tensor tensor stream)))
-
-(defmethod print-object ((tensor standard-matrix) stream)
- (print-unreadable-object (tensor stream :type t)
- (format stream "~A x ~A~%" (nrows tensor) (ncols tensor))
- (print-tensor tensor stream)))
diff --git a/src/real-tensor.lisp b/src/real-tensor.lisp
index 8016330..c6c252c 100644
--- a/src/real-tensor.lisp
+++ b/src/real-tensor.lisp
@@ -1,14 +1,12 @@
(in-package :matlisp)
-(eval-when (load eval compile)
- (deftype real-type ()
- "The type of the elements stored in a REAL-MATRIX"
- 'double-float)
-
- (deftype real-array (size)
- "The type of the storage structure for a REAL-MATRIX"
- `(simple-array real-type (,size)))
- )
+(deftype real-type ()
+ "The type of the elements stored in a REAL-MATRIX"
+ 'double-float)
+
+(deftype real-array (size)
+ "The type of the storage structure for a REAL-MATRIX"
+ `(simple-array real-type (,size)))
;;
(make-array-allocator allocate-real-store 'real-type 0d0
@@ -64,14 +62,4 @@ Allocates real storage. Default initial-element = 0d0.")
element stream)
(format stream "~11,5,,,,,'Eg" element))
-;;
-
-(defun make-real-tensor-dims (&rest subs)
- (let* ((dims (make-index-store subs))
- (ss (reduce #'* dims))
- (store (allocate-real-store ss)))
- (make-instance 'real-tensor :store store :dimensions dims)))
-#+nil(defun make-real-tensor-array (arr)
- (let* ((dims (array-dimensions arr))
- (ret (apply #'make-real-tensor-dims dims)))))
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
deleted file mode 100644
index 192a23a..0000000
--- a/src/standard-matrix.lisp
+++ /dev/null
@@ -1,131 +0,0 @@
-(in-package :matlisp)
-
-;;
-(defclass standard-matrix (standard-tensor)
- ((rank
- :accessor rank
- :type index-type
- :initform 2
- :documentation "For a matrix, rank = 2."))
- (:documentation "Basic matrix class."))
-
-(definline nrows (matrix)
- (declare (type standard-matrix matrix))
- (aref (dimensions matrix) 0))
-
-(definline ncols (matrix)
- (declare (type standard-matrix matrix))
- (aref (dimensions matrix) 1))
-
-(definline row-stride (matrix)
- (declare (type standard-matrix matrix))
- (aref (strides matrix) 0))
-
-(definline col-stride (matrix)
- (declare (type standard-matrix matrix))
- (aref (strides matrix) 1))
-
-(definline size (matrix)
- (declare (type standard-matrix matrix))
- (let ((dims (dimensions matrix)))
- (declare (type (index-array 2) dims))
- (list (aref dims 0) (aref dims 1))))
-
-;;
-(defmethod initialize-instance :after ((matrix standard-matrix) &rest initargs)
- (declare (ignore initargs))
- (mlet*
- ((rank (rank matrix) :type index-type))
- (unless (= rank 2)
- (error 'tensor-not-matrix :rank rank :tensor matrix))))
-
-;;
-(definline row-matrix-p (matrix)
- "
- Syntax
- ======
- (ROW-MATRIX-P x)
-
- Purpose
- =======
- Return T if X is a row matrix (number of columns is 1)"
- (tensor-type-p matrix '(1 *)))
-
-(definline col-matrix-p (matrix)
- "
- Syntax
- ======
- (COL-MATRIX-P x)
-
- Purpose
- =======
- Return T if X is a column matrix (number of rows is 1)"
- (tensor-type-p matrix '(* 1)))
-
-(definline row-or-col-matrix-p (matrix)
-"
- Syntax
- ======
- (ROW-OR-COL-matrix-P x)
-
- Purpose
- =======
- Return T if X is either a row or a column matrix."
- (or (row-vector-p matrix) (col-vector-p matrix)))
-
-(defun square-matrix-p (matrix)
- (and (square-p matrix) (matrix-p matrix)))
-
-;;
-(defgeneric fill-matrix (matrix fill-element)
- (:documentation
- "
- Syntax
- ======
- (FILL-MATRIX matrix fill-element)
-
- Purpose
- =======
- Fill MATRIX with FILL-ELEMENT.
-"))
-
-(defmethod fill-matrix ((matrix t) (fill t))
- (error "arguments MATRIX and FILL to FILL-MATRIX must be a
-matrix and a number"))
-
-;;
-(defclass real-matrix (standard-matrix real-tensor)
- ()
- (:documentation "A class of matrices with real elements."))
-
-(defclass real-sub-matrix (real-matrix standard-sub-tensor)
- ()
- (:documentation "Sub-matrix class with real elements."))
-
-(setf (gethash 'real-matrix *sub-tensor-counterclass*) 'real-sub-matrix
- (gethash 'real-sub-matrix *sub-tensor-counterclass*) 'real-sub-matrix
- ;;
- (gethash 'real-matrix *tensor-class-optimizations*) 'real-tensor
- (gethash 'real-sub-matrix *tensor-class-optimizations*) 'real-tensor)
-;;
-
-(defclass complex-matrix (standard-matrix complex-tensor)
- ()
- (:documentation "A class of matrices with complex elements."))
-
-(defclass complex-sub-matrix (complex-matrix standard-sub-tensor)
- ()
- (:documentation "Sub-matrix class with complex elements."))
-
-(setf (gethash 'complex-matrix *sub-tensor-counterclass*) 'complex-sub-matrix
- (gethash 'complex-sub-matrix *sub-tensor-counterclass*) 'complex-sub-matrix
- ;;
- (gethash 'complex-matrix *tensor-class-optimizations*) 'complex-tensor
- (gethash 'complex-sub-matrix *tensor-class-optimizations*) 'complex-tensor)
-
-;;
-
-(definline matrix-ref (matrix row &optional col)
- (declare (type standard-matrix matrix))
- (tensor-ref matrix `(,row ,col)))
-
diff --git a/src/standard-tensor.lisp b/src/standard-tensor.lisp
index a780560..94db9a7 100644
--- a/src/standard-tensor.lisp
+++ b/src/standard-tensor.lisp
@@ -89,6 +89,8 @@
"
Contains a either:
o A property list containing:
+ :store-allocator (n) -> Allocates a store of size n
+ :coercer (ele) -> Coerced to store-type
:element-type
:store-type
:reader (store idx) => result
@@ -101,9 +103,9 @@
(declare (type symbol clname))
(let ((opt (gethash clname *tensor-class-optimizations*)))
(cond
+ ((null opt) nil)
((symbolp opt)
(get-tensor-class-optimization opt))
- ((null opt) nil)
(t (values opt clname)))))
;; Akshay: I have no idea what this does, or why we want it
@@ -471,3 +473,4 @@
(error 'tensor-cannot-find-sub-class :tensor-class (class-of tensor)))
:parent-tensor tensor :store (store tensor) :head nhd
:dimensions (make-index-store ndim) :strides (make-index-store nstd)))))))
+
diff --git a/src/tensor-maker.lisp b/src/tensor-maker.lisp
new file mode 100644
index 0000000..da9fdd8
--- /dev/null
+++ b/src/tensor-maker.lisp
@@ -0,0 +1,48 @@
+(in-package :matlisp)
+
+(defmacro make-tensor-maker (func-name (tensor-class))
+ (let ((opt (get-tensor-class-optimization tensor-class)))
+ (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
+ `(defun ,func-name (&rest args)
+ (labels ((make-dims (dims)
+ (declare (type cons dims))
+ (let* ((vdim (make-index-store dims))
+ (ss (reduce #'* vdim))
+ (store (,(getf opt :store-allocator) ss)))
+ (make-instance ',tensor-class :store store :dimensions vdim)))
+ (make-from-array (arr)
+ (declare (type (array * *) arr))
+ (let* ((ret (make-dims (array-dimensions arr)))
+ (st-r (store ret)))
+ (declare (type ,tensor-class ret)
+ (type ,(linear-array-type (getf opt :store-type)) st-r))
+ (very-quickly
+ (mod-dotimes (idx (dimensions ret))
+ with (linear-sums
+ (of-r (strides ret) (head ret)))
+ do ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) (apply #'aref arr (idx->list idx))) 'st-r 'of-r)))
+ ret))
+ (make-from-list (lst)
+ (let* ((ret (make-dims (list-dimensions lst)))
+ (st-r (store ret)))
+ (declare (type ,tensor-class ret)
+ (type ,(linear-array-type (getf opt :store-type)) st-r))
+ (very-quickly
+ (list-loop (idx ele lst)
+ with (linear-sums
+ (of-r (strides ret)))
+ do ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) ele) 'st-r 'of-r)))
+ ret)))
+ (let ((largs (length args)))
+ (if (= largs 1)
+ (etypecase (first args)
+ (array
+ (make-from-array (first args)))
+ (cons
+ (make-from-list (first args)))
+ (integer
+ (make-dims (list (first args)))))
+ (make-dims args)))))))
+
+(make-tensor-maker make-real-tensor (real-tensor))
+(make-tensor-maker make-complex-tensor (complex-tensor))
diff --git a/src/utilities.lisp b/src/utilities.lisp
index fec194e..a35196c 100644
--- a/src/utilities.lisp
+++ b/src/utilities.lisp
@@ -231,6 +231,17 @@
(with-output-to-string (ostr ret)
(apply #'format (append `(,ostr ,fmt) args)))
ret))
+
+(defun list-dimensions (lst)
+ (declare (type list lst))
+ (labels ((lst-tread (idx lst)
+ (if (null lst) (reverse idx)
+ (progn
+ (setf (car idx) (length lst))
+ (if (consp (car lst))
+ (lst-tread (cons 0 idx) (car lst))
+ (reverse idx))))))
+ (lst-tread (list 0) lst)))
;;---------------------------------------------------------------;;
(defstruct (foreign-vector
(:conc-name fv-)
commit 174d27300595c21a466a330fa34ab66fa7131bdf
Author: Akshay Srinivasan <aks...@gm...>
Date: Wed Jul 4 22:19:51 2012 +0530
More changes to README
diff --git a/README.org b/README.org
index 46ed60f..071d4fc 100644
--- a/README.org
+++ b/README.org
@@ -1,20 +1,20 @@
MatLisp - a base for scientific computation in Lisp.
-
This is the development branch of Matlisp.
* Progress Tracker
** What works ?
* Basic {real, complex} tensor structure in place.
- * Added a specialisation agnostic macro which generate...
[truncated message content] |
