From: Akshay S. <ak...@us...> - 2013-01-26 19:07:39
|
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 a9eca9b0cc4287b81b325360972175771d7f6c71 (commit) from 4862a338530bb1b435f2d6535913abe9947931b6 (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 a9eca9b0cc4287b81b325360972175771d7f6c71 Author: Akshay Srinivasan <aks...@gm...> Date: Sat Jan 26 11:00:46 2013 -0800 o gemm! now copies matrices and then calls gemm if they have weird incompatible strides (replacing the elaborate scheme from before of using gemv). o define-tensor now takes a new field "value-incfer". diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 101c7cb..60b3d8d 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -351,9 +351,9 @@ ((tensor-class element-type store-element-type store-type &rest class-decls) &key f+ f- finv+ fid+ f* f/ finv* fid* fconj f= matrix vector - store-allocator coercer coercer-unforgiving reader value-writer reader-writer swapper) + store-allocator coercer coercer-unforgiving reader value-writer value-incfer reader-writer swapper) ;;Error checking - (assert (and f+ f- finv+ fid+ f* f/ finv* fid* f= store-allocator coercer coercer-unforgiving matrix vector reader value-writer reader-writer swapper)) + (assert (and f+ f- finv+ fid+ f* f/ finv* fid* f= store-allocator coercer coercer-unforgiving matrix vector reader value-writer value-incfer reader-writer swapper)) ;; `(eval-when (:compile-toplevel :load-toplevel :execute) ;;Class definitions @@ -399,6 +399,7 @@ :fconj ',fconj :reader ',reader :value-writer ',value-writer + :value-incfer ',value-incfer :reader-writer ',reader-writer :swapper ',swapper :store-allocator ',store-allocator diff --git a/src/classes/complex-tensor.lisp b/src/classes/complex-tensor.lisp index 276e9f5..382641a 100644 --- a/src/classes/complex-tensor.lisp +++ b/src/classes/complex-tensor.lisp @@ -89,6 +89,13 @@ (setf (aref store (* 2 idx)) (realpart value) (aref store (1+ (* 2 idx))) (imagpart value))) +(definline complex-type.value-incfer (value store idx) + (declare (type complex-store-vector store) + (type index-type idx) + (type complex-type value)) + (incf (aref store (* 2 idx)) (realpart value)) + (incf (aref store (1+ (* 2 idx))) (imagpart value))) + (definline complex-type.reader-writer (fstore fidx tstore tidx) (declare (type complex-store-vector fstore tstore) (type index-type fidx tidx)) @@ -125,6 +132,7 @@ ;; :reader complex-type.reader :value-writer complex-type.value-writer + :value-incfer complex-type.value-incfer :reader-writer complex-type.reader-writer :swapper complex-type.swapper) diff --git a/src/classes/real-tensor.lisp b/src/classes/real-tensor.lisp index 4924b69..ac5d144 100644 --- a/src/classes/real-tensor.lisp +++ b/src/classes/real-tensor.lisp @@ -55,6 +55,12 @@ (type real-type value)) (setf (aref store idx) value)) +(definline real-type.value-incfer (value store idx) + (declare (type index-type idx) + (type real-store-vector store) + (type real-type value)) + (incf (aref store idx) value)) + (definline real-type.reader-writer (fstore fidx tstore tidx) (declare (type index-type fidx tidx) (type real-store-vector fstore tstore)) @@ -99,6 +105,7 @@ Allocates real storage. Default initial-element = 0d0.") ;; :reader real-type.reader :value-writer real-type.value-writer + :value-incfer real-type.value-incfer :reader-writer real-type.reader-writer :swapper real-type.swapper) diff --git a/src/classes/symbolic-tensor.lisp b/src/classes/symbolic-tensor.lisp index d837051..0365929 100644 --- a/src/classes/symbolic-tensor.lisp +++ b/src/classes/symbolic-tensor.lisp @@ -67,6 +67,12 @@ (type symbolic-type value)) (setf (aref store idx) value)) +(definline symbolic-type.value-incfer (value store idx) + (declare (type index-type idx) + (type symbolic-store-vector store) + (type symbolic-type value)) + (setf (aref store idx) (symbolic-type.f+ (aref store idx) value))) + (definline symbolic-type.reader-writer (fstore fidx tstore tidx) (declare (type index-type fidx tidx) (type symbolic-store-vector fstore tstore)) @@ -110,6 +116,7 @@ Allocates symbolic storage. Default initial-element = 0.") ;; :reader symbolic-type.reader :value-writer symbolic-type.value-writer + :value-incfer symbolic-type.value-incfer :reader-writer symbolic-type.reader-writer :swapper symbolic-type.swapper) diff --git a/src/level-1/realimag.lisp b/src/level-1/realimag.lisp index b38569c..db9566a 100644 --- a/src/level-1/realimag.lisp +++ b/src/level-1/realimag.lisp @@ -44,7 +44,7 @@ (etypecase tensor (real-tensor tensor) (complex-tensor (make-instance (ecase (rank tensor) (2 'real-matrix) (1 'real-vector) (t 'real-tensor)) - :parent-tensor tensor :store (store tensor) :store-size (store-size tensor) + :parent-tensor tensor :store (store tensor) :store-size (length (store tensor)) :dimensions (dimensions tensor) :strides (map 'index-store-vector #'(lambda (x) (* 2 x)) (strides tensor)) :head (the index-type (* 2 (head tensor))))) @@ -66,7 +66,7 @@ (etypecase tensor (real-tensor tensor) (complex-tensor (make-instance (ecase (rank tensor) (2 'real-matrix) (1 'real-vector) (t 'real-tensor)) - :parent-tensor tensor :store (store tensor) :store-size (store-size tensor) + :parent-tensor tensor :store (store tensor) :store-size (length (store tensor)) :dimensions (dimensions tensor) :strides (map 'index-store-vector #'(lambda (x) (* 2 x)) (strides tensor)) :head (the index-type (+ 1 (* 2 (head tensor)))))) diff --git a/src/level-1/tensor-maker.lisp b/src/level-1/tensor-maker.lisp index 5961798..37c0b20 100644 --- a/src/level-1/tensor-maker.lisp +++ b/src/level-1/tensor-maker.lisp @@ -59,3 +59,24 @@ ;;Had to move it here in the wait for copy! (definline sub-tensor (tensor subscripts &optional (preserve-rank nil)) (copy (sub-tensor~ tensor subscripts preserve-rank))) + +(defmacro make-zeros-dims (func-name (tensor-class)) + (let ((opt (get-tensor-class-optimization-hashtable tensor-class))) + (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) + `(eval-when (:compile-toplevel :load-toplevel :execute) + (let ((opt (get-tensor-class-optimization-hashtable ',tensor-class))) + (assert opt nil 'tensor-cannot-find-optimization :tensor-class ',tensor-class) + (setf (getf opt :zero-maker) ',func-name + (get-tensor-class-optimization ',tensor-class) opt)) + (defun ,func-name (dims) + (declare (type index-store-vector dims)) + (let-typed ((rnk (length dims) :type index-type) + (size (very-quickly (lvec-foldl #'(lambda (a b) (declare (type index-type a b)) (the index-type (* a b))) dims)))) + (make-instance (case rnk (2 ',(getf opt :matrix)) (1 ',(getf opt :vector)) (t ',tensor-class)) + :dimensions (copy-seq dims) :store (,(getf opt :store-allocator) size) :store-size size)))))) + +(make-zeros-dims real-typed-zeros (real-tensor)) +(make-zeros-dims complex-typed-zeros (complex-tensor)) + +#+maxima +(make-zeros-dims symbolc-typed-tensor (symbolic-tensor)) diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp index 727e89b..c7ad77a 100644 --- a/src/level-2/gemv.lisp +++ b/src/level-2/gemv.lisp @@ -50,23 +50,25 @@ (Aval (,(getf opt :reader) sto-A of-A) :type ,(getf opt :element-type))) (setf dot (,(getf opt :f+) dot (,(getf opt :f*) xval Aval)))) :finally (,(getf opt :value-writer) (,(getf opt :f+) (,(getf opt :f*) alpha dot) val) sto-y of-y)))))))) - (if blas-gemv-func + (if blas-gemv-func `(mlet* ((call-fortran? (> (max (nrows A) (ncols A)) ,fortran-call-lb)) - ((maj-A ld-A fop-A) (if call-fortran? (blas-matrix-compatible-p A job) (values nil 0 "?")) :type (symbol index-type (string 1)))) + ((maj-A ld-A fop-A) (blas-matrix-compatible-p A job) :type (symbol index-type (string 1)))) (cond - ((and maj-a call-fortran?) - (let-typed ((nr-A (nrows A) :type index-type) - (nc-A (ncols A) :type index-type)) - (when (eq maj-A :row-major) - (rotatef nr-A nc-A)) - (,blas-gemv-func fop-a nr-A nc-A - alpha (store A) ld-A - (store x) (aref (strides x) 0) - beta - (store y) (aref (strides y) 0) - (head A) (head x) (head y)))) - (t + (call-fortran? + (if maj-A + (let-typed ((nr-A (nrows A) :type index-type) + (nc-A (ncols A) :type index-type)) + (when (eq maj-A :row-major) + (rotatef nr-A nc-A)) + (,blas-gemv-func fop-a nr-A nc-A + alpha (store A) ld-A + (store x) (aref (strides x) 0) + beta + (store y) (aref (strides y) 0) + (head A) (head x) (head y))) + (,func alpha (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A))) x beta y job))) + (t ,lisp-routine))) lisp-routine)) y)))) diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp index 24a7519..48831ee 100644 --- a/src/level-3/gemm.lisp +++ b/src/level-3/gemm.lisp @@ -28,11 +28,11 @@ (in-package #:matlisp) -(defmacro generate-typed-gemm! (func (tensor-class blas-gemm-func blas-gemv-func fortran-lb-parameter)) +(defmacro generate-typed-gemm! (func (tensor-class blas-gemm-func fortran-lb-parameter)) (let* ((opt (if-ret (get-tensor-class-optimization-hashtable tensor-class) (error 'tensor-cannot-find-optimization :tensor-class tensor-class))) (matrix-class (getf opt :matrix)) - (blas? (and blas-gemm-func blas-gemv-func))) + (blas? blas-gemm-func)) `(eval-when (:compile-toplevel :load-toplevel :execute) (let ((opt (get-tensor-class-optimization-hashtable ',tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class ',tensor-class) @@ -72,40 +72,66 @@ (unless (,(getf opt :f=) beta (,(getf opt :fid*))) (,(getf opt :num-scal) beta C)) ;; - (let-typed ((of-A hd-A :type index-type) - (of-B hd-B :type index-type) - (of-C hd-C :type index-type) - (r.cstp-C (* cstp-C nc-C) :type index-type) - (d.rstp-B (- rstp-B (* cstp-B nc-C)) :type index-type) - (d.rstp-A (- rstp-A (* cstp-A dotl)) :type index-type)) - (very-quickly - (loop :repeat nr-C + (cond + ((and (= cstp-C 1) (= cstp-B 1) nil) + (let-typed ((of-A hd-A :type index-type) + (of-B hd-B :type index-type) + (of-C hd-C :type index-type) + (d.rstp-B (- rstp-B nc-C) :type index-type) + (d.rstp-A (- rstp-A (* cstp-A dotl)) :type index-type)) + (very-quickly + (loop :repeat nr-C :do (progn (loop :repeat dotl - :do (let-typed - ((ele-A (,(getf opt :f*) alpha (,(getf opt :reader) sto-A of-A)) :type ,(getf opt :element-type))) - (loop :repeat nc-C - :do (progn - (,(getf opt :value-writer) - (,(getf opt :f+) - (,(getf opt :reader) sto-C of-C) - (,(getf opt :f*) ele-A (,(getf opt :reader) sto-B of-B))) - sto-C of-C) - (incf of-C cstp-C) - (incf of-B cstp-B))) - (decf of-C r.cstp-C) - (incf of-A cstp-A) - (incf of-B d.rstp-B))) + :do (let-typed + ((ele-A (,(getf opt :f*) alpha (,(getf opt :reader) sto-A of-A)) :type ,(getf opt :element-type))) + (loop :repeat nc-C + :do (progn + (,(getf opt :value-incfer) (,(getf opt :f*) ele-A (,(getf opt :reader) sto-B of-B)) + sto-C of-C) + (incf of-C) + (incf of-B))) + (decf of-C nc-C) + (incf of-A cstp-A) + (incf of-B d.rstp-B))) (incf of-C rstp-C) (incf of-A d.rstp-A) - (setf of-B hd-B)))))))) + (setf of-B hd-B)))))) + (t + (let-typed ((of-A hd-A :type index-type) + (of-B hd-B :type index-type) + (of-C hd-C :type index-type) + (r.cstp-C (* cstp-C nc-C) :type index-type) + (d.rstp-B (- rstp-B (* cstp-B nc-C)) :type index-type) + (d.rstp-A (- rstp-A (* cstp-A dotl)) :type index-type)) + (very-quickly + (loop :repeat nr-C + :do (progn + (loop :repeat dotl + :do (let-typed + ((ele-A (,(getf opt :f*) alpha (,(getf opt :reader) sto-A of-A)) :type ,(getf opt :element-type))) + (loop :repeat nc-C + :do (progn + (,(getf opt :value-writer) + (,(getf opt :f+) + (,(getf opt :reader) sto-C of-C) + (,(getf opt :f*) ele-A (,(getf opt :reader) sto-B of-B))) + sto-C of-C) + (incf of-C cstp-C) + (incf of-B cstp-B))) + (decf of-C r.cstp-C) + (incf of-A cstp-A) + (incf of-B d.rstp-B))) + (incf of-C rstp-C) + (incf of-A d.rstp-A) + (setf of-B hd-B)))))))))) ;;Tie together Fortran and lisp-routines. `(mlet* (((job-A job-B) (ecase job (:nn (values :n :n)) (:nt (values :n :t)) (:tn (values :t :n)) (:tt (values :t :t))) - :type (symbol symbol)) + :type (symbol symbol)) ,@(when blas? `((call-fortran? (> (max (nrows C) (ncols C) (if (eq job-A :n) (ncols A) (nrows A))) ,fortran-lb-parameter)) @@ -114,79 +140,37 @@ ((maj-C ld-C fop-C) (blas-matrix-compatible-p C :n) :type (symbol index-type nil))))) ,(if blas? `(cond - ((and call-fortran? maj-A maj-B maj-C) - (let-typed ((nr-C (nrows C) :type index-type) - (nc-C (ncols C) :type index-type) - (dotl (ecase job-A (:n (ncols A)) (:t (nrows A))) :type index-type)) - (when (eq maj-C :row-major) - (rotatef A B) - (rotatef ld-A ld-B) - (rotatef maj-A maj-B) - (rotatef nr-C nc-C) - (setf (values fop-A fop-B) - (values (fortran-snop fop-B) (fortran-snop fop-A)))) - (,blas-gemm-func fop-A fop-B nr-C nc-C dotl - alpha (store A) ld-A (store B) ld-B - beta (store C) ld-C - (head A) (head B) (head C)))) - ((and call-fortran? maj-A) - (let-typed ((nc-C (ncols C) :type index-type) - (strd-C (col-stride C) :type index-type) - (stp-C (row-stride C) :type index-type) - (sto-C (store C) :type ,(linear-array-type (getf opt :store-type))) - ; - (nr-A (nrows A) :type index-type) - (nc-A (ncols A) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - (hd-A (head A) :type index-type) - ; - (stp-B (if (eq job-B :n) (row-stride B) (col-stride B)) :type index-type) - (sto-B (store B) :type ,(linear-array-type (getf opt :store-type))) - (strd-B (if (eq job-B :n) (col-stride B) (row-stride B)) :type index-type)) - (when (eq maj-A :row-major) - (rotatef nr-A nc-A)) - (very-quickly - (loop repeat nc-C - for of-B of-type index-type = (head B) then (+ of-B strd-B) - for of-C of-type index-type = (head C) then (+ of-C strd-C) - do (,blas-gemv-func fop-A nr-A nc-A - alpha sto-A ld-A - sto-B stp-B - beta sto-C stp-C - hd-A of-B of-C))))) - ((and call-fortran? maj-B) - (let-typed ((nr-C (nrows C) :type index-type) - (stp-C (col-stride C) :type index-type) - (strd-C (row-stride C) :type index-type) - (sto-C (store c) :type ,(linear-array-type (getf opt :store-type))) - ; - (stp-A (if (eq job-A :n) (col-stride A) (row-stride A)) :type index-type) - (strd-A (if (eq job-A :n) (row-stride A) (col-stride A)) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - ; - (nr-B (nrows B) :type index-type) - (nc-B (ncols B) :type index-type) - (hd-B (head B) :type index-type) - (fop-B (fortran-snop fop-B) :type (string 1)) - (sto-B (store B) :type ,(linear-array-type (getf opt :store-type)))) - (when (eq maj-B :row-major) - (rotatef nr-B nc-B)) - (very-quickly - (loop repeat nr-C - for of-A of-type index-type = (head A) then (+ of-A strd-A) - for of-C of-type index-type = (head C) then (+ of-C strd-C) - do (,blas-gemv-func fop-B nr-B nc-B - alpha sto-B ld-B - sto-A stp-A - beta sto-C stp-C - hd-B of-A of-C))))) + (call-fortran? + (if (and maj-A maj-B maj-C) + (let-typed ((nr-C (nrows C) :type index-type) + (nc-C (ncols C) :type index-type) + (dotl (ecase job-A (:n (ncols A)) (:t (nrows A))) :type index-type)) + (when (eq maj-C :row-major) + (rotatef A B) + (rotatef ld-A ld-B) + (rotatef maj-A maj-B) + (rotatef nr-C nc-C) + (setf (values fop-A fop-B) + (values (fortran-snop fop-B) (fortran-snop fop-A)))) + (,blas-gemm-func fop-A fop-B nr-C nc-C dotl + alpha (store A) ld-A (store B) ld-B + beta (store C) ld-C + (head A) (head B) (head C))) + (let ((ret (,func alpha + (if maj-A A (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A)))) + (if maj-B B (,(getf opt :copy) B (,(getf opt :zero-maker) (dimensions B)))) + beta + (if maj-C C (,(getf opt :copy) C (,(getf opt :zero-maker) (dimensions C)))) + job))) + (unless maj-C + (,(getf opt :copy) ret C))))) (t ,lisp-routine)) lisp-routine))) C)))) ;;Real (generate-typed-gemm! real-base-typed-gemm! - (real-tensor dgemm dgemv *real-l3-fcall-lb*)) + (real-tensor dgemm *real-l3-fcall-lb*)) (definline real-typed-gemm! (alpha A B beta C job) (real-base-typed-gemm! alpha A B beta C @@ -197,7 +181,7 @@ ;;Complex (generate-typed-gemm! complex-base-typed-gemm! - (complex-tensor zgemm zgemv *complex-l3-fcall-lb*)) + (complex-tensor zgemm *complex-l3-fcall-lb*)) (definline complex-typed-gemm! (alpha A B beta C job) (declare (type complex-matrix A B C) @@ -229,7 +213,7 @@ ;;Symbolic #+maxima (generate-typed-gemm! symbolic-base-typed-gemm! - (symbolic-tensor nil nil 0)) + (symbolic-tensor nil 0)) ;;---------------------------------------------------------------;; ----------------------------------------------------------------------- Summary of changes: src/base/standard-tensor.lisp | 5 +- src/classes/complex-tensor.lisp | 8 ++ src/classes/real-tensor.lisp | 7 ++ src/classes/symbolic-tensor.lisp | 7 ++ src/level-1/realimag.lisp | 4 +- src/level-1/tensor-maker.lisp | 21 +++++ src/level-2/gemv.lisp | 30 ++++--- src/level-3/gemm.lisp | 174 +++++++++++++++++--------------------- 8 files changed, 143 insertions(+), 113 deletions(-) hooks/post-receive -- matlisp |