From: Akshay S. <ak...@us...> - 2012-06-30 19:05:56
|
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 695636685fd91ce1602b135d0c0e782ca06d47e7 (commit) via adf78b01d61996d75fda7ce045e3f3f11aa3f9ed (commit) from 1231d97cfa4e89109805a7a5284d939bbd65f5f9 (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 695636685fd91ce1602b135d0c0e782ca06d47e7 Author: Akshay Srinivasan <aks...@gm...> Date: Sat Jun 30 23:35:35 2012 +0530 o Optimised mod-dotimes (finally!). A 1000x1000 matrix multiplication is now a mere 2 times slower than fully optimized(-Ofast) C code using incremental index-offsets (like mod-dotimes), but about 5 times slower than naive-3-loops in Fortran. OpenBLAS is still much much faster. This more than justifies jumping from Python to Lisp (as if the language itself was not enough!). Timings ======= SBCL: 12.5 s C: 6.0 s Fortran: 2.5 s OpenBLAS: 0.3 s o realimag, copy and scal are now in a usable form. o Tweaks to permutations. Added sort-permute function. diff --git a/README b/README.old similarity index 100% rename from README rename to README.old diff --git a/matlisp.asd b/matlisp.asd index 787e4ca..e300c13 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -85,23 +85,35 @@ :depends-on ("foreign-interface" "foreign-functions") :components ((:file "conditions") - (:file "standard-tensor") + ;; + (:file "standard-tensor" + :depends-on ("conditions")) + ;; (:file "loopy" :depends-on ("standard-tensor")) (:file "blas-helpers" :depends-on ("standard-tensor")) + ;; (:file "real-tensor" :depends-on ("standard-tensor")) (:file "complex-tensor" :depends-on ("standard-tensor")) (:file "standard-matrix" - :depends-on ("standard-tensor")) + :depends-on ("standard-tensor" "real-tensor" "complex-tensor")) ;; (:file "real-matrix" ;; :depends-on ("standard-matrix")) ;; (:file "complex-matrix" ;; :depends-on ("standard-matrix")) (:file "print" - :depends-on ("standard-tensor" "standard-matrix")))))) + :depends-on ("standard-tensor" "standard-matrix")) + ;;Copy, Scal + (:file "copy" + :depends-on ("real-tensor" "complex-tensor" "loopy")) + (:file "scal" + :depends-on ("copy" "loopy")) + (:file "realimag" + :depends-on ("real-tensor" "complex-tensor" "copy")) + )))) ;; (defclass f2cl-cl-source-file (asdf:cl-source-file) diff --git a/src/blas-helpers.lisp b/src/blas-helpers.lisp index 10dc90b..3817137 100644 --- a/src/blas-helpers.lisp +++ b/src/blas-helpers.lisp @@ -1,9 +1,11 @@ (in-package :matlisp) (definline idx-max (seq) + (declare (type (index-array *) seq)) (reduce #'max seq)) (definline idx-min (seq) + (declare (type (index-array *) seq)) (reduce #'min seq)) (defun idx= (a b) @@ -22,12 +24,23 @@ (loop for ele across a collect ele)) -(defun blas-copyable-p (&rest tensors) - (let ((stdi-list (very-quickly +(defun blas-copyable-p (ten-a ten-b) + ;; (declare (type standard-tensor ten-a ten-b)) + ;; (let ((stdi-a (very-quickly + ;; (sort (apply #'vector + ;; (loop + ;; for std across (strides ten-a) + ;; and dim across (dimensions ten-a) + ;; collect `(,std ,dim))) + ;; #'< :key #'first)))) + ;; t)) + + + (let ((stdi-list (very-quickly (loop - for ten in tensors + for ten of-type standard-tensor in tensors and pten = nil then ten - for i = 0 then (1+ i) + for i of-type index-type = 0 then (1+ i) when (> i 0) do (unless (idx= (dimensions ten) (dimensions pten)) (return nil)) @@ -37,8 +50,8 @@ (very-quickly (sort (apply #'vector (loop - for std across (strides ten) - and dim across (dimensions ten) + for std of-type index-type across (strides ten) + and dim of-type index-type across (dimensions ten) collect `(,std ,dim))) #'< :key #'car))))))) (if (null stdi-list) (values nil nil) @@ -46,7 +59,7 @@ (loop for stdi in stdi-list and p-stdi = (first stdi-list) then stdi - for i = 0 then (1+ i) + for i of-type index-type = 0 then (1+ i) when (> i 0) do (unless (loop for a-stdi across stdi diff --git a/src/copy.lisp b/src/copy.lisp index e620330..95bcaf2 100644 --- a/src/copy.lisp +++ b/src/copy.lisp @@ -143,6 +143,7 @@ (generate-typed-copy! complex-typed-copy! (complex-tensor zcopy)) (generate-typed-num-copy! complex-typed-num-copy! (complex-tensor zcopy)) ;;---------------------------------------------------------------;; + (defgeneric copy! (from-tensor to-tensor) (:documentation " @@ -182,18 +183,27 @@ (defmethod copy! ((x number) (y real-tensor)) (real-typed-num-copy! (coerce-real x) y)) -(defmethod copy! ((x complex-matrix) (y complex-tensor)) +(defmethod copy! ((x complex-tensor) (y complex-tensor)) (complex-typed-copy! x y)) -(defmethod copy! ((x real-matrix) (y complex-tensor)) - (real-double-copy!-typed x (mrealpart~ y)) - (scal! 0d0 (mimagpart~ y)) +(defmethod copy! ((x real-tensor) (y complex-tensor)) + ;;Borrowed from realimag.lisp + (let ((tmp (make-instance 'real-sub-tensor + :parent-tensor y :store (store y) + :dimensions (dimensions y) + :strides (map '(index-array *) #'(lambda (n) (* 2 n)) (strides y)) + :head (the index-type (* 2 (head y)))))) + (declare (type real-sub-tensor tmp)) + (real-typed-copy! x tmp) + ;;Increasing the head by 1 points us to the imaginary part. + (incf (head tmp)) + (real-typed-num-copy! 0d0 tmp)) y) (defmethod copy! ((x number) (y complex-tensor)) (complex-typed-num-copy! (coerce-complex x) y)) -;;;; +;; (defgeneric copy (tensor) (:documentation " @@ -206,16 +216,12 @@ Return a copy of the tensor X")) (defmethod copy ((tensor real-tensor)) - (let* ((ret (apply #'make-real-tensor-dims - (loop for dim across (dimensions tensor) - collect dim)))) + (let* ((ret (apply #'make-real-tensor-dims (idx->list (dimensions tensor))))) (declare (type real-tensor ret)) (copy! tensor ret))) (defmethod copy ((tensor complex-tensor)) - (let* ((ret (apply #'make-complex-tensor-dims - (loop for dim across (dimensions tensor) - collect dim)))) + (let* ((ret (apply #'make-complex-tensor-dims (idx->list (dimensions tensor))))) (declare (type complex-tensor ret)) (copy! tensor ret))) diff --git a/src/loopy.lisp b/src/loopy.lisp index 9b27ad1..32037bc 100644 --- a/src/loopy.lisp +++ b/src/loopy.lisp @@ -76,18 +76,19 @@ (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 (dims-sym rank-sym count-sym) - `(let* ((,dims-sym ,dims) - (,rank-sym (length ,dims-sym)) - (,idx (allocate-index-store ,rank-sym)) - ,@(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))) - ,@(let ((decl `(,@(when (getf sdecl :linear-sums) - `((type (index-array *) ,@(mapcar #'(lambda (x) (getf x :stride-sym)) (getf sdecl :linear-sums))))) - ,@(loop for x in (getf sdecl :variables) + `(let ((,dims-sym ,dims)) + (declare (type (index-array *) ,dims-sym)) + (let ((,rank-sym (length ,dims-sym))) + (declare (type index-type ,rank-sym)) + (let ((,idx (allocate-index-store ,rank-sym)) + ,@(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))))) + ,@(loop for x in (getf sdecl :variables) unless (null (getf x :type)) - collect `(type ,(getf x :type) ,(getf x :variable)))))) - (unless (null decl) - `((declare ,@decl)))) + collect `(type ,(getf x :type) ,(getf x :variable)))) (loop ,@(loop for decl in (getf sdecl :linear-sums) append `(with ,(getf decl :offset-sym) of-type index-type = ,(getf decl :offset-init))) ,@(unless (null code) @@ -107,7 +108,7 @@ `(let ((,cstrd (aref ,(getf decl :stride-sym) ,count-sym))) (declare (type index-type ,cstrd)) (unless (= ,cstrd 0) - (decf ,(getf decl :offset-sym) (* ,cstrd (1- (aref ,dims-sym ,count-sym))))))))) + (decf ,(getf decl :offset-sym) (the index-type (* ,cstrd (1- (aref ,dims-sym ,count-sym)))))))))) (progn (incf (aref ,idx ,count-sym)) ,@(loop @@ -118,4 +119,4 @@ (unless (= ,cstrd 0) (incf ,(getf decl :offset-sym) ,cstrd))))) (return t))) - finally (return nil)))))))))) + finally (return nil)))))))))))) diff --git a/src/permutations.lisp b/src/permutations.lisp index 797e198..c68d97d 100644 --- a/src/permutations.lisp +++ b/src/permutations.lisp @@ -4,7 +4,40 @@ ((message :reader message :initform "Object is not a permutation.")) (:documentation "Object is not a permutation.")) -;;---------------------------------------------------------------;; +;;Class definitions----------------------------------------------;; +(defclass permutation () + ((representation :accessor repr + :initarg :repr) + (group-rank :accessor group-rank + :type index-type))) +;; +(defclass permutation-cycle (permutation) + ((representation :type cons))) + +(defmethod initialize-instance :after ((per permutation-cycle) &rest initargs) + (declare (ignore initargs)) + (let ((cls 0)) + (declare (type index-type cls)) + (unless (very-quickly + (dolist (cyc (r-value per) t) + (unless (cycle-p cyc) + (return nil)) + (setf cls (max cls (idx-max cyc))))) + (error 'permutation-error)) + (setf (group-rank per) (the index-type (1+ cls))))) +;; +(defclass permutation-action (permutation) + ((representation :type (index-array *)))) + +(defmethod initialize-instance :after ((per permutation-action) &rest initargs) + (declare (ignore initargs)) + (let ((act (r-value per))) + (declare (type (index-array *) act)) + (unless (action-p act) + (error 'permutation-error)) + (setf (group-rank per) (idx-max act)))) + +;;Conversions and validation-------------------------------------;; (defun insert-element (x sort l-b u-b) "Does a binary-esque sort to keep track of elements in a permutation, in descending order. If there are duplicates @@ -17,47 +50,61 @@ (declare (type index-type l-b u-b)) (let* ((midx (+ l-b (floor (- u-b l-b) 2))) (mid (aref sort midx))) + (declare (type index-type midx mid)) (cond ((or (< x 0) (member x `(,(aref sort u-b) ,(aref sort l-b) ,mid))) (error 'permutation-error)) ((= midx l-b) (when (> x (aref sort u-b)) - (loop - with sidx = (+ midx (if (> x mid) 0 1)) - for i downfrom (- len 1) to sidx - do (setf (aref sort (+ i 1)) (aref sort i)) - finally (setf (aref sort sidx) x)))) + (very-quickly + (loop + with sidx of-type index-type = (+ midx (if (> x mid) 0 1)) + for i of-type index-type downfrom (1- len) to sidx + do (setf (aref sort (+ i 1)) (aref sort i)) + finally (setf (aref sort sidx) x))))) ((< x mid) (insert-ele midx u-b)) ((> x mid) (insert-ele l-b midx))) sort))) (insert-ele l-b u-b)))) -(defun cycle-p (perm) +(defun cycle-new-p (perm) "Does a sorting operation to check for duplicate elements in the cycle representation of a permutation." + (declare (type (index-array *) perm)) (let* ((len (length perm)) - (sort (allocate-index-store len -1))) - (dotimes (i len t) - (handler-case (insert-element (aref perm i) sort 0 i) - (permutation-error () (return nil)))))) - -(defun action-p (arr) + (sort (very-quickly (sort (copy-seq perm) #'<)))) + (declare (type (index-array *) sort) + (type index-type len)) + (very-quickly + (loop for i of-type index-type from 1 below len + when (= (aref sort i) (aref sort (1- i))) + do (return nil) + finally (return t))))) + +(defun action-p (act) "Checks if ARR is a possible permutation vector. A permutation pi is characterized by a vector containing the indices from 0,..., @function{length}(@arg{perm})-1 in some order." - (declare (type (index-array *) arr)) - (let ((s-arr (sort (copy-seq arr) #'<))) - (dotimes (i (length s-arr) t) - (unless (= i (aref s-arr i)) - (return nil))))) - -(defun action->cycle (per) + (declare (type (index-array *) act)) + (let* ((len (length act)) + (sort (very-quickly (sort (copy-seq act) #'<)))) + (declare (type (index-array *) sort) + (type index-type len)) + (very-quickly + (loop for i of-type index-type from 0 below len + unless (= (aref sort i) i) + do (return nil) + finally (return t))))) + +(defun action->cycle (act) ;;Caution: will go into an infinite loop if object is not proper. - "This function obtains the canonical cycle representation + " + This function obtains the canonical cycle representation of a permutation. The first argument is the action of the permutation on the array #(0 1 2 3 ..). \"Canonical\" may be a bit of an overstatement; this is the way - S_n was presented by Van der Waerden." + S_n was presented by Van der Waerden. +" (declare (type permutation-action per)) (mlet* ((arr (r-value per) :type (index-array *))) @@ -85,37 +132,9 @@ (cycle-walk nil nil)))) ;;---------------------------------------------------------------;; -(defclass permutation () - ((representation :accessor r-value - :initarg :r-value) - (group-rank :accessor group-rank - :type index-type))) - -(defclass permutation-cycle (permutation) - ((representation :type cons))) - -(defmethod initialize-instance :after ((per permutation-cycle) &rest initargs) - (declare (ignore initargs)) - (let ((cls 0)) - (unless (dolist (cyc (r-value per) t) - (unless (cycle-p cyc) - (return nil)) - (setf cls (max cls (reduce #'max cyc)))) - (error 'permutation-error)) - (setf (group-rank per) (1+ cls)))) - -(defclass permutation-action (permutation) - ((:representation :type (index-array *)))) - -(defmethod initialize-instance :after ((per permutation-action) &rest initargs) - (declare (ignore initargs)) - (unless (action-p (r-value per)) - (error 'permutation-error))) (defun cycles->action (cyc) ) - - ;; (defun apply-cycle! (seq cyc) @@ -152,4 +171,60 @@ (defun seqrnd (seq) "Randomize the elements of a sequence. Destructive on SEQ." - (sort seq #'> :key #'(lambda (x) (random 1.0)))) \ No newline at end of file + (sort seq #'> :key #'(lambda (x) (random 1.0)))) + +;; + +(defun allocate-unit-permutation (n) + (declare (type fixnum n)) + (let ((ret (allocate-index-store n))) + (declare (type (index-array *) ret)) + (very-quickly + (loop + for i of-type index-type from 0 below n + do (setf (aref ret i) i))) + ret)) + +(defun sort-permute (seq predicate) + " + (sort-permute seq predicate) + + Sorts a index-array and also returns + the permutation-action required to move + from the given sequence to the sorted form. + + Takes about 10x the running time which can be + achieved with cl:sort. + " + (declare (type (index-array *) seq) + (type function predicate)) + (let* ((len (length seq)) + (perm (allocate-unit-permutation len))) + (declare (type index-type len) + (type (index-array *) perm)) + (labels ((qsort-bounds (lb ub) + (declare (type index-type lb ub)) + #+nil(format t "~a lb:~a ub:~a ~%" seq lb ub) + (if (= ub (1+ lb)) t + (let* ((ele (aref seq lb)) + (ele-idx (very-quickly + (loop + for i of-type index-type from (1+ lb) below ub + with ele-idx of-type index-type = lb + do (unless (funcall predicate ele (aref seq i)) + (when (> i (1+ ele-idx)) + (rotatef (aref seq ele-idx) (aref seq (1+ ele-idx))) + (rotatef (aref perm ele-idx) (aref perm (1+ ele-idx)))) + (rotatef (aref seq ele-idx) (aref seq i)) + (rotatef (aref perm ele-idx) (aref perm i)) + (incf ele-idx) + #+nil(format t " ~a ~%" seq)) + finally (return ele-idx))))) + (when (> (- ub ele-idx) 2) + (qsort-bounds (1+ ele-idx) ub)) + (when (> (- ele-idx lb) 1) + (qsort-bounds lb ele-idx)))))) + (qsort-bounds 0 len) + (values seq perm)))) + +(quicksort-with-action (idxv 10 9 8 7 6 5 4 3 2 1) #'<) diff --git a/src/realimag.lisp b/src/realimag.lisp index 015ae12..4dbc7c3 100644 --- a/src/realimag.lisp +++ b/src/realimag.lisp @@ -84,7 +84,7 @@ (complex-tensor (make-instance 'real-sub-tensor :parent-tensor tensor :store (store tensor) :dimensions (dimensions tensor) - :strides (map '(index-array *) #'(lambda (x) (* 2 x)) (strides xten)) + :strides (map '(index-array *) #'(lambda (x) (* 2 x)) (strides tensor)) :head (the index-type (* 2 (head tensor))))) (number (realpart tensor)))) @@ -106,38 +106,31 @@ (complex-tensor (make-instance 'real-sub-tensor :parent-tensor tensor :store (store tensor) :dimensions (dimensions tensor) - :strides (map '(index-array *) #'(lambda (x) (* 2 x)) (strides xten)) + :strides (map '(index-array *) #'(lambda (x) (* 2 x)) (strides tensor)) :head (the index-type (+ 1 (* 2 (head tensor)))))) (number (imagpart tensor)))) -(defun tensor-realpart (mat) +(definline tensor-realpart (tensor) " Syntax ====== - (MREALPART matrix) + (tensor-realpart tensor) Purpose ======= - Returns a copy of the real part of \"matrix\". + Returns a copy of the real part of tensor. - If \"matrix\" is a scalar, returns its real part. + If \"tensor\" is a scalar, returns its real part. See IMAG, REALPART, IMAGPART " - (typecase mat - (real-matrix (copy mat)) - (complex-matrix (copy (make-instance 'sub-real-matrix - :parent mat :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))))) - (number (cl:realpart mat)))) + (copy (tensor-realpart~ tensor))) -(defun tensor-imagpart (mat) +(definline tensor-imagpart (tensor) " Syntax ====== - (MIMAGPART~ matrix) + (tensor-imagpart matrix) Purpose ======= @@ -147,45 +140,4 @@ See IMAG, REALPART, IMAGPART " - - (typecase mat - (real-matrix (make-real-matrix-dim (nrows mat) (ncols mat))) - (complex-matrix (copy (make-instance 'sub-real-matrix - :parent mat :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)))))) - (number (cl:imagpart mat)))) - - -(declaim (inline real)) -(defun real (matrix) -" - Syntax - ====== - (REAL matrix) - - Purpose - ======= - Returns a new REAL-MATRIX which is the real part of MATRIX. - If MATRIX is a scalar, returns its real part. - - See IMAG, REALPART, IMAGPART -" - (mrealpart matrix)) - - -(defun imag (matrix) -" - Syntax - ====== - (IMAG matrix) - - Purpose - ======= - Returns a new REAL-MATRIX which is the imaginary part of MATRIX. - If MATRIX is a scalar, returns its imaginary part. - - See REAL, REALPART, IMAGPART -" - (mimagpart matrix)) + (copy (tensor-imagpart tensor))) diff --git a/src/scal.lisp b/src/scal.lisp index 15c6ea6..134d168 100644 --- a/src/scal.lisp +++ b/src/scal.lisp @@ -75,22 +75,17 @@ `(defun ,func (alpha to) (declare (type ,tensor-class to) (type ,(getf opt :element-type) alpha)) - (let ((t-dims (dimensions to)) - (t-stds (strides to)) - (t-sto (store to)) - (t-hd (head to))) - (declare (type (index-array *) t-dims t-stds) - (type index-type t-hd) - (type ,(linear-array-type (getf opt :store-type)) t-sto)) - (if-let (min-stride (consecutive-store-p t-stds t-dims)) - (,blas-func (number-of-elements to) alpha t-sto min-stride t-hd) + (if-let (min-stride (consecutive-store-p to)) + (,blas-func (number-of-elements to) alpha (store to) min-stride (head to)) + (let ((t-sto (store to))) + (declare (type ,(linear-array-type (getf opt :store-type)) t-sto)) (very-quickly ;;Can possibly make this faster (x2) by using ,blas-func in one of ;;the inner loops, but this is to me messy and as of now unnecessary. ;;SBCL can already achieve Fortran-ish speed inside this loop. - (mod-dotimes (idx t-dims) + (mod-dotimes (idx (dimensions to)) with (linear-sums - (t-of t-stds t-hd)) + (t-of (strides to) (head to))) do (let ((scal-val (* ,(funcall (getf opt :reader) 't-sto 't-of) alpha))) ,(funcall (getf opt :value-writer) 'scal-val 't-sto 't-of)))))) to))) @@ -99,8 +94,8 @@ ;; zdscal and zscal is significant, except for very large arrays. (generate-typed-scal! real-typed-scal! (real-tensor dscal)) (generate-typed-scal! complex-typed-scal! (complex-tensor zscal)) - ;;---------------------------------------------------------------;; + (defgeneric scal! (alpha x) (:documentation " @@ -144,14 +139,12 @@ (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)) +(defmethod scal ((alpha complex) (x real-tensor)) + (let* ((result (apply #'make-complex-tensor-dims (idx->list (dimensions x))))) + (declare (type complex-tensor result)) (copy! x result) (scal! alpha result))) -(defmethod scal ((alpha number) (x complex-matrix)) +(defmethod scal ((alpha number) (x complex-tensor)) (let ((result (copy x))) (scal! alpha result))) diff --git a/src/standard-tensor.lisp b/src/standard-tensor.lisp index ae18f0b..905fe8e 100644 --- a/src/standard-tensor.lisp +++ b/src/standard-tensor.lisp @@ -1,30 +1,10 @@ (in-package :matlisp) -;; -(eval-when (load eval compile) - (deftype integer4-type () - '(signed-byte 32)) - (deftype integer4-array (size) - `(simple-array integer4-type (,size))) - - ;; - (deftype index-type () - #+cmu '(signed-byte 32) - #-cmu '(signed-byte 64)) - (deftype index-array (size) - `(simple-array index-type (,size))) - ) - -(declaim (inline allocate-integer4-store)) -(make-array-allocator allocate-integer4-store 'integer4-type 0 -" - Syntax - ====== - (ALLOCATE-INT32-STORE SIZE [INITIAL-ELEMENT 0]) +(deftype index-type () + 'fixnum) - Purpose - ======= - Allocates integer-32 storage.") +(deftype index-array (size) + `(simple-array index-type (,size))) (make-array-allocator allocate-index-store 'index-type 0 " @@ -51,8 +31,8 @@ (definline idxv (&rest contents) (make-index-store contents)) - ;; + (defclass standard-tensor () ((rank :accessor rank @@ -161,12 +141,13 @@ (very-quickly (loop for i of-type index-type from 0 below rank - and sto-idx of-type index-type = hd then (+ sto-idx (* cidx (aref strides i))) - for cidx of-type index-type = (aref idx i) - do (unless (< -1 cidx (aref dims i)) - (error 'tensor-index-out-of-bounds :argument i :index cidx :dimension (aref dims i))) + for cidx across idx + with sto-idx of-type index-type = hd + do (if (< -1 cidx (aref dims i)) + (incf sto-idx (the index-type (* (aref strides i) cidx))) + (error 'tensor-index-out-of-bounds :argument i :index cidx :dimension (aref dims i))) finally (return sto-idx)))))) - +x (defun store-indexing-lst (idx hd strides dims) " Syntax diff --git a/src/utilities.lisp b/src/utilities.lisp index 6441673..63fef05 100644 --- a/src/utilities.lisp +++ b/src/utilities.lisp @@ -237,6 +237,7 @@ (defstruct (foreign-vector (:conc-name fv-) (:print-function (lambda (obj stream depth) + (declare (ignore depth)) (format stream "#F(") (let ((sz (fv-size obj))) (dotimes (i sz) commit adf78b01d61996d75fda7ce045e3f3f11aa3f9ed Author: Akshay Srinivasan <aks...@gm...> Date: Fri Jun 29 21:34:48 2012 +0530 More tweaks to "copy" diff --git a/matlisp.asd b/matlisp.asd index 3ff85aa..787e4ca 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -88,6 +88,8 @@ (:file "standard-tensor") (:file "loopy" :depends-on ("standard-tensor")) + (:file "blas-helpers" + :depends-on ("standard-tensor")) (:file "real-tensor" :depends-on ("standard-tensor")) (:file "complex-tensor" diff --git a/src/blas-helpers.lisp b/src/blas-helpers.lisp index 78c33cb..10dc90b 100644 --- a/src/blas-helpers.lisp +++ b/src/blas-helpers.lisp @@ -1,89 +1,108 @@ (in-package :matlisp) -(definline fortran-op (op) - (ecase op (:n "N") (:t "T"))) +(definline idx-max (seq) + (reduce #'max seq)) -(definline fortran-nop (op) - (ecase op (:t "N") (:n "T"))) +(definline idx-min (seq) + (reduce #'min seq)) -(defun fortran-snop (sop) - (cond - ((string= sop "N") "T") - ((string= sop "T") "N") - (t (error "Unrecognised fortran-op.")))) +(defun idx= (a b) + (declare (type (index-array *) a b)) + (when (= (length a) (length b)) + (very-quickly + (loop + for ele-a across a + for ele-b across b + unless (= ele-a ele-b) + do (return nil) + finally (return t))))) -(defun blas-copyable-p (matrix) - (declare (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)) - (very-quickly - (cond - ((or (= nc 1) (= cs (* nr rs))) (values t rs ne)) - ((or (= nr 1) (= rs (* nc cs))) (values t cs ne)) - (t (values nil -1 -1)))))) +(definline idx->list (a) + (declare (type (index-array *) a)) + (loop for ele across a + collect ele)) -(defun blas-matrix-compatible-p (matrix &optional (op :n)) - (declare (optimize (safety 0) (speed 3)) - (type (or real-matrix complex-matrix) matrix)) - (mlet* (((rs cs) (slot-values matrix '(row-stride col-stride)) - :type (fixnum fixnum))) - (cond - ((= cs 1) (values :row-major rs (fortran-nop op))) - ((= rs 1) (values :col-major cs (fortran-op op))) - ;;Lets not confound lisp's type declaration. - (t (values nil -1 "?"))))) +(defun blas-copyable-p (&rest tensors) + (let ((stdi-list (very-quickly + (loop + for ten in tensors + and pten = nil then ten + for i = 0 then (1+ i) + when (> i 0) + do (unless (idx= (dimensions ten) (dimensions pten)) + (return nil)) + collect (progn + (assert (typep ten 'standard-tensor) nil + 'invalid-type :given (type-of ten) :expected 'standard-tensor) + (very-quickly + (sort (apply #'vector + (loop + for std across (strides ten) + and dim across (dimensions ten) + collect `(,std ,dim))) + #'< :key #'car))))))) + (if (null stdi-list) (values nil nil) + (very-quickly + (loop + for stdi in stdi-list + and p-stdi = (first stdi-list) then stdi + for i = 0 then (1+ i) + when (> i 0) + do (unless (loop + for a-stdi across stdi + and a-aoff = (first (aref stdi 0)) then (* a-aoff (second a-stdi)) + for b-stdi across p-stdi + and b-aoff = (first (aref p-stdi 0)) then (* b-aoff (second b-stdi)) + do (unless (and (= (first a-stdi) a-aoff) + (= (first b-stdi) b-aoff) + (= (second a-stdi) (second b-stdi))) + (return nil)) + finally (return t)) + (return (values t nil))) + finally (return (values t (mapcar #'(lambda (x) (first (aref x 0))) stdi-list)))))))) +(defun consecutive-store-p (tensor) + (declare (type standard-tensor tensor)) + (let ((strides (strides tensor)) + (dims (dimensions tensor))) + (declare (type (index-array *) strides dims)) + (let* ((stride-dims (very-quickly + (sort (apply #'vector + (loop + for std across strides + and dim across dims + collect `(,std ,dim))) + #'< :key #'car))) + (stride-min (first (aref stride-dims 0)))) + (declare (type index-type stride-min) + (type (simple-vector *) stride-dims)) + (very-quickly + (loop + for st-di across stride-dims + and accumulated-off = stride-min then (* accumulated-off (second st-di)) + unless (= (first st-di) accumulated-off) do (return nil) + finally (return stride-min)))))) -(defun col-major-p (strides dims) - (declare (type (index-array *) strides dims)) - (very-quickly - (loop - for off across strides - and dim across dims - and accumulated-off = 1 then (* accumulated-off dim) - unless (= off accumulated-off) do (return nil) - finally (return t)))) -(defun row-major-p (strides dims) - (declare (type (index-array *) strides dims)) - (very-quickly - (loop - for idx of-type index-type from (1- (length dims)) downto 0 - for dim of-type index-type = (aref dims idx) - for off of-type index-type = (aref strides idx) - and accumulated-off of-type index-type = 1 then (* accumulated-off dim) - unless (= off accumulated-off) do (return nil) - finally (return t)))) +;; (defun blas-matrix-compatible-p (matrix &optional (op :n)) +;; (declare (optimize (safety 0) (speed 3)) +;; (type (or real-matrix complex-matrix) matrix)) +;; (mlet* (((rs cs) (slot-values matrix '(row-stride col-stride)) +;; :type (fixnum fixnum))) +;; (cond +;; ((= cs 1) (values :row-major rs (fortran-nop op))) +;; ((= rs 1) (values :col-major cs (fortran-op op))) +;; ;;Lets not confound lisp's type declaration. +;; (t (values nil -1 "?"))))) -(defun same-dimension-p (a b) - (declare (type (index-array *) a b)) - (let ((l-a (length a))) - (when (= l-a (length b)) - (very-quickly - (loop - for i from 0 below l-a - unless (= (aref a i) (aref b i)) - do (return nil) - finally (return t)))))) +;; (definline fortran-op (op) +;; (ecase op (:n "N") (:t "T"))) -(defun consecutive-store-p (strides dims) - (declare (type (index-array *) strides dims)) - (let* ((stride-dims (very-quickly - (sort (apply #'vector - (loop - for std across strides - and dim across dims - collect `(,std ,dim))) - #'< :key #'car))) - (stride-min (first (aref stride-dims 0)))) - (declare (type index-type stride-min) - (type (simple-vector *) stride-dims)) - (very-quickly - (loop - for st-di across stride-dims - and accumulated-off = stride-min then (* accumulated-off (second st-di)) - unless (= (first st-di) accumulated-off) do (return nil) - finally (return stride-min))))) +;; (definline fortran-nop (op) +;; (ecase op (:t "N") (:n "T"))) + +;; (defun fortran-snop (sop) +;; (cond +;; ((string= sop "N") "T") +;; ((string= sop "T") "N") +;; (t (error "Unrecognised fortran-op.")))) diff --git a/src/conditions.lisp b/src/conditions.lisp index 781f435..ad16f1e 100644 --- a/src/conditions.lisp +++ b/src/conditions.lisp @@ -55,11 +55,10 @@ (to :reader to :initarg :to)) (:documentation "Cannot coerce one type into another.")) -(defmethod print-object ((c coercion) stream) +(defmethod print-object ((c coercion-error) stream) (format stream "Cannot coerce ~a into ~a." (from c) (to c)) (call-next-method)) - ;;---------------------------------------------------------------;; (define-condition matlisp-error (error) ;;Optional argument for error-handling. diff --git a/src/copy.lisp b/src/copy.lisp index 4c21264..e620330 100644 --- a/src/copy.lisp +++ b/src/copy.lisp @@ -86,30 +86,24 @@ (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) `(defun ,func (from to) (declare (type ,tensor-class from to)) - (let ((f-dims (dimensions from)) - (f-stds (strides from)) - (f-sto (store from)) - (f-hd (head from)) - (t-dims (dimensions to)) - (t-stds (strides to)) - (t-sto (store to)) - (t-hd (head to))) - (declare (type (index-array *) f-dims f-stds t-dims t-stds) - (type index-type f-hd t-hd) - (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto)) - (if (or (and (row-major-p t-stds t-dims) (row-major-p f-stds f-dims)) - (and (col-major-p t-stds t-dims) (col-major-p f-stds f-dims))) - (,blas-func (number-of-elements from) f-sto 1 t-sto 1 f-hd t-hd) - (very-quickly - ;;Can possibly make this faster (x2) by using ,blas-func in one of - ;;the inner loops, but this is to me messy and as of now unnecessary. - ;;SBCL can already achieve Fortran-ish speed inside this loop. - (mod-dotimes (idx f-dims) - with (linear-sums - (f-of f-stds f-hd) - (t-of t-stds t-hd)) - do ,(funcall (getf opt :reader-writer) 'f-sto 'f-of 't-sto 't-of)))) - to)))) + (multiple-value-bind (dims-p strd-p) (blas-copyable-p from to) + (unless dims-p + (error 'tensor-dimension-mismatch)) + (if strd-p + (,blas-func (number-of-elements from) (store from) (first strd-p) (store to) (second strd-p) (head from) (head to)) + (let ((f-sto (store from)) + (t-sto (store to))) + (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto)) + (very-quickly + ;;Can possibly make this faster (x2) by using ,blas-func in one of + ;;the inner loops, but this is to me messy and as of now unnecessary. + ;;SBCL can already achieve Fortran-ish speed inside this loop. + (mod-dotimes (idx (dimensions from)) + with (linear-sums + (f-of (strides from) (head from)) + (t-of (strides to) (head to))) + do ,(funcall (getf opt :reader-writer) 'f-sto 'f-of 't-sto 't-of)))))) + to))) (defmacro generate-typed-num-copy! (func (tensor-class blas-func)) ;;Be very careful when using functions generated by this macro. @@ -128,19 +122,19 @@ (declare (type (index-array *) t-dims t-stds) (type index-type t-hd) (type ,(linear-array-type (getf opt :store-type)) t-sto)) - (if (consecutive-p t-stds t-dims) - (let ((num-array (,(getf opt :store-allocator) 1))) - (declare (type ,(linear-array-type (getf opt :store-type)) num-array)) - ,(funcall (getf opt :value-writer) 'num-from 'num-array 0) - (,blas-func (number-of-elements to) num-array 0 t-sto 1 0 t-hd)) - (very-quickly - ;;Can possibly make this faster (x2) by using ,blas-func in one of - ;;the inner loops, but this is to me messy and as of now unnecessary. - ;;SBCL can already achieve Fortran-ish speed inside this loop. - (mod-dotimes (idx t-dims) - with (linear-sums - (t-of t-stds t-hd)) - do ,(funcall (getf opt :value-writer) 'num-from 't-sto 't-of)))) + (if-let (min-stride (consecutive-store-p to)) + (let ((num-array (,(getf opt :store-allocator) 1))) + (declare (type ,(linear-array-type (getf opt :store-type)) num-array)) + ,(funcall (getf opt :value-writer) 'num-from 'num-array 0) + (,blas-func (number-of-elements to) num-array 0 t-sto min-stride 0 t-hd)) + (very-quickly + ;;Can possibly make this faster (x2) by using ,blas-func in one of + ;;the inner loops, but this is to me messy and as of now unnecessary. + ;;SBCL can already achieve Fortran-ish speed inside this loop. + (mod-dotimes (idx t-dims) + with (linear-sums + (t-of t-stds t-hd)) + do ,(funcall (getf opt :value-writer) 'num-from 't-sto 't-of)))) to)))) (generate-typed-copy! real-typed-copy! (real-tensor dcopy)) @@ -173,7 +167,7 @@ REAL-MATRIX but the converse is possible. ") (:method :before ((x standard-tensor) (y standard-tensor)) - (unless (same-dimension-p (dimensions x) (dimensions y)) + (unless (idx= (dimensions x) (dimensions y)) (error 'tensor-dimension-mismatch))) (:method ((x standard-tensor) (y standard-tensor)) (mod-dotimes (idx (dimensions x)) @@ -191,16 +185,16 @@ (defmethod copy! ((x complex-matrix) (y complex-tensor)) (complex-typed-copy! x y)) -;; (defmethod copy! ((x real-matrix) (y complex-tensor)) -;; (real-double-copy!-typed x (mrealpart~ y)) -;; (scal! 0d0 (mimagpart~ y)) -;; y) +(defmethod copy! ((x real-matrix) (y complex-tensor)) + (real-double-copy!-typed x (mrealpart~ y)) + (scal! 0d0 (mimagpart~ y)) + y) (defmethod copy! ((x number) (y complex-tensor)) (complex-typed-num-copy! (coerce-complex x) y)) ;;;; -(defgeneric copy (matrix) +(defgeneric copy (tensor) (:documentation " Syntax @@ -209,72 +203,44 @@ Purpose ======= - Return a copy of the matrix X")) + Return a copy of the tensor X")) -(defmethod copy ((matrix real-matrix)) - (let* ((n (nrows matrix)) - (m (ncols matrix)) - (result (make-real-matrix-dim n m))) - (declare (type fixnum n m)) - (copy! matrix result))) +(defmethod copy ((tensor real-tensor)) + (let* ((ret (apply #'make-real-tensor-dims + (loop for dim across (dimensions tensor) + collect dim)))) + (declare (type real-tensor ret)) + (copy! tensor ret))) -(defmethod copy ((matrix complex-matrix)) - (let* ((n (nrows matrix)) - (m (ncols matrix)) - (result (make-complex-matrix-dim n m))) - (declare (type fixnum n m)) - (copy! matrix result))) +(defmethod copy ((tensor complex-tensor)) + (let* ((ret (apply #'make-complex-tensor-dims + (loop for dim across (dimensions tensor) + collect dim)))) + (declare (type complex-tensor ret)) + (copy! tensor ret))) -(defmethod copy ((matrix number)) - matrix) +(defmethod copy ((tensor number)) + tensor) ;; -(defgeneric convert-to-lisp-array (matrix) - (:documentation - " +(defun convert-to-lisp-array (tensor) +" Syntax ====== - (CONVERT-TO-LISP-ARRAY matrix) + (convert-to-lisp-array tensor) Purpose ======= - Create a new Lisp array with the same dimensions as the matrix and - with the same elements. This is a copy of the matrix. - - Row and column vectors are converted to a 1D lisp vector. Other - matrices are converted a 2D lisp array. -")) - -(defun convert-1d-array (m eltype) - (let ((array (make-array (* (number-of-rows m) - (number-of-cols m)) - :element-type eltype))) - ;; We could do this faster by accessing the storage directly, but - ;; this is easy. - (dotimes (k (length array)) - (setf (aref array k) (matrix-ref m k))) - array)) - -(defun convert-2d-array (m eltype) - (let* ((nrows (number-of-rows m)) - (ncols (number-of-cols m)) - (array (make-array (list (number-of-rows m) - (number-of-cols m)) - :element-type eltype))) - ;; We could do this faster by accessing the storage directly, but - ;; this is easy. - (dotimes (r nrows) - (dotimes (c ncols) - (setf (aref array r c) - (matrix-ref m r c)))) - array)) - -(defmethod convert-to-lisp-array ((m real-matrix)) - (if (or (row-vector-p m) (col-vector-p m)) - (convert-1d-array m 'double-float) - (convert-2d-array m 'double-float))) - -(defmethod convert-to-lisp-array ((m complex-matrix)) - (if (or (row-vector-p m) (col-vector-p m)) - (convert-1d-array m '(complex double-float)) - (convert-2d-array m '(complex double-float)))) + Create a new Lisp array with the same dimensions as the tensor and + with the same elements. This is a copy of the tensor. +" + (declare (type standard-tensor tensor)) + (let* ((dims (dimensions tensor)) + (ret (make-array (idx->list dims) + :element-type (if-ret (getf (get-tensor-class-optimization (class-name (class-of tensor))) :element-type) + (error 'tensor-cannot-find-optimization :tensor-class (class-name (class-of tensor))))))) + (declare (type (index-array *) dims)) + (very-quickly + (mod-dotimes (idx dims) + do (setf (apply #'aref ret (idx->list idx)) (tensor-ref tensor idx)))) + ret)) diff --git a/src/realimag.lisp b/src/realimag.lisp index 2036121..015ae12 100644 --- a/src/realimag.lisp +++ b/src/realimag.lisp @@ -133,7 +133,7 @@ :head (* 2 (head mat))))) (number (cl:realpart mat)))) -(defun mimagpart (mat) +(defun tensor-imagpart (mat) " Syntax ====== diff --git a/src/scal.lisp b/src/scal.lisp index e0de476..15c6ea6 100644 --- a/src/scal.lisp +++ b/src/scal.lisp @@ -140,8 +140,7 @@ (defmethod scal ((alpha number) (x number)) (* alpha x)) -;; -(defmethod scal ((alpha cl:real) (x real-matrix)) +(defmethod scal ((alpha number) (x real-tensor)) (let ((result (copy x))) (scal! alpha result))) @@ -153,7 +152,6 @@ (copy! x result) (scal! alpha result))) -;; (defmethod scal ((alpha number) (x complex-matrix)) (let ((result (copy x))) (scal! alpha result))) diff --git a/src/standard-tensor.lisp b/src/standard-tensor.lisp index 9e9bcd5..ae18f0b 100644 --- a/src/standard-tensor.lisp +++ b/src/standard-tensor.lisp @@ -8,9 +8,9 @@ `(simple-array integer4-type (,size))) ;; - (deftype index-type () + (deftype index-type () #+cmu '(signed-byte 32) - #-cmu '(signed-byte 64)) + #-cmu '(signed-byte 64)) (deftype index-array (size) `(simple-array index-type (,size))) ) @@ -189,7 +189,7 @@ (type cons idx)) (let ((rank (length strides))) (declare (type index-type rank)) - (labels ((rec-sum (sum i lst) + (labels ((rec-sum (sum i lst) (cond ((consp lst) (let ((cidx (car lst))) @@ -219,7 +219,7 @@ HD + \ STRIDES * IDX /_ i i i = 0 -" +" (declare (type standard-tensor tensor) (type (or (index-array *) cons) idx)) (typecase idx @@ -308,7 +308,7 @@ (let ((,tstore (store ,tensym))) (declare (type ,(linear-array-type store-element-type) ,tstore)) ,@body)))) - (let ((hst (list + (let ((hst (list :reader (macrofy ,reader) :value-writer (macrofy ,value-writer) :reader-writer (macrofy ,reader-writer) @@ -408,7 +408,7 @@ nil))))))) (parse-sub subscripts 0))))) -(definline vector-p (tensor) +(definline vector-p (tensor) (declare (type standard-tensor tensor)) (tensor-type-p tensor '(*))) @@ -425,7 +425,7 @@ ;;---------------------------------------------------------------;; (define-constant +array-slicing-symbols+ '(\:) -" +" Symbols which are used to refer to slicing operations.") (defun sub-tensor~ (tensor subscripts) @@ -453,7 +453,7 @@ ;; Get [:, :, 0:10:2] (0:10:2 = [i : 0 <= i < 10, i % 2 = 0]) > (sub-tensor~ X '(\: \: ((\: 2) 0 *))) -" +" (declare (type standard-tensor tensor)) (let ((rank (rank tensor)) (dims (dimensions tensor)) diff --git a/src/utilities.lisp b/src/utilities.lisp index 3a26d44..6441673 100644 --- a/src/utilities.lisp +++ b/src/utilities.lisp @@ -233,15 +233,6 @@ (with-output-to-string (ostr ret) (apply #'format (append `(,ostr ,fmt) args))) ret)) - -(declaim (inline seq-max)) -(defun seq-max (seq) - (reduce #'max seq)) - -(declaim (inline seq-max)) -(defun seq-min (seq) - (reduce #'min seq)) - ;;---------------------------------------------------------------;; (defstruct (foreign-vector (:conc-name fv-) ----------------------------------------------------------------------- Summary of changes: README => README.old | 0 matlisp.asd | 20 ++++- src/blas-helpers.lisp | 186 +++++++++++++++++++++++++++------------------- src/conditions.lisp | 3 +- src/copy.lisp | 180 +++++++++++++++++++-------------------------- src/loopy.lisp | 27 ++++--- src/permutations.lisp | 177 +++++++++++++++++++++++++++++++------------- src/realimag.lisp | 68 +++-------------- src/scal.lisp | 33 +++----- src/standard-tensor.lisp | 53 ++++--------- src/utilities.lisp | 10 +-- 11 files changed, 383 insertions(+), 374 deletions(-) rename README => README.old (100%) hooks/post-receive -- matlisp |