Menu
▾
▴
matlisp-commit
|
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
|
Thanks for helping keep SourceForge clean.
X