From: Akshay S. <ak...@us...> - 2012-12-28 16:47:27
|
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 e916823ab6bd97795ad7eaea63ad778a423b0919 (commit) via 856c60140465482aaa72f021360c4c795073ad6f (commit) via 0166ce8014b662aca4a91484eb2458e7e87be8ac (commit) via 1c74913ff22ddc869220e6ee124bcf272b188d12 (commit) via 376399e23fbbb868c8eb3ef80ee8bc9c65c5d98e (commit) via ff3082257b6f984b30131dd170f011eacd78f7e6 (commit) from c0ba7e46b390f6e744e6865192b6eb57eb95c585 (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 e916823ab6bd97795ad7eaea63ad778a423b0919 Author: Akshay Srinivasan <aks...@gm...> Date: Fri Dec 28 10:41:02 2012 -0600 o Ported GEMM. Should get to work on automatic method generation. diff --git a/matlisp.asd b/matlisp.asd index 17f0086..cbb1979 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -145,12 +145,10 @@ :depends-on ("copy" "scal")) (:file "trans" :depends-on ("scal" "copy")))) - #+nil (:module "matlisp-level-2" :pathname "level-2" :depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1") :components ((:file "gemv"))) - #+nil (:module "matlisp-level-3" :pathname "level-3" :depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1") diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp index e6ba5ab..ecf1756 100644 --- a/src/level-3/gemm.lisp +++ b/src/level-3/gemm.lisp @@ -28,133 +28,144 @@ (in-package #:matlisp) -(defmacro generate-typed-gemm! (func (matrix-class blas-gemm-func blas-gemv-func fortran-lb-parameter)) - (let* ((opt (get-tensor-class-optimization matrix-class))) - (assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class) +(defmacro generate-typed-gemm! (func (tensor-class blas-gemm-func blas-gemv-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))) `(definline ,func (alpha A B beta C job) (declare (type ,(getf opt :element-type) alpha beta) (type ,matrix-class A B C) (type symbol job)) - (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)) - (call-fortran? (> (max (nrows C) (ncols C) (if (eq job-A :n) (ncols A) (nrows A))) - ,fortran-lb-parameter)) - ((maj-A ld-A fop-A) (if call-fortran? (blas-matrix-compatible-p A job-A) (values nil 0 "?")) :type (symbol index-type (string 1))) - ((maj-B ld-B fop-B) (if call-fortran? (blas-matrix-compatible-p B job-B) (values nil 0 "?")) :type (symbol index-type (string 1))) - ((maj-C ld-C fop-C) (if call-fortran? (blas-matrix-compatible-p C :n) (values nil 0 "?")) :type (symbol index-type nil))) - (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))) + ,(let + ((lisp-routine + `(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) ; - (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))) + (rstp-A (row-stride A) :type index-type) + (cstp-A (col-stride A) :type index-type) + (hd-A (head 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))))) - (t - (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) + (rstp-B (row-stride B) :type index-type) + (cstp-B (col-stride B) :type index-type) + (hd-B (head B) :type index-type) + (sto-B (store B) :type ,(linear-array-type (getf opt :store-type))) ; - (rstp-A (row-stride A) :type index-type) - (cstp-A (col-stride A) :type index-type) - (hd-A (head A) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) + (rstp-C (row-stride C) :type index-type) + (cstp-C (col-stride C) :type index-type) + (hd-C (head C) :type index-type) + (sto-C (store C) :type ,(linear-array-type (getf opt :store-type)))) + (when (eq job-A :t) + (rotatef rstp-A cstp-A)) + (when (eq job-B :t) + (rotatef rstp-B cstp-B)) + (very-quickly + (loop :repeat nr-C + :for rof-A :of-type index-type := hd-A :then (+ rof-A rstp-A) + :for rof-C :of-type index-type := hd-C :then (+ rof-C rstp-C) + :do (loop :repeat nc-C + :for cof-B :of-type index-type := hd-B :then (+ cof-B cstp-B) + :for of-C :of-type index-type := rof-C :then (+ of-C cstp-C) + :do (let-typed ((val (,(getf opt :f*) beta (,(getf opt :reader) sto-C of-C)) :type ,(getf opt :element-type))) + (loop :repeat dotl + :for of-A :of-type index-type := rof-A :then (+ of-A cstp-A) + :for of-B :of-type index-type := cof-B :then (+ of-B rstp-B) + :with sum :of-type ,(getf opt :element-type) := (,(getf opt :fid+)) + :do (let-typed ((A-val (,(getf opt :reader) sto-A of-A) :type ,(getf opt :element-type)) + (B-val (,(getf opt :reader) sto-B of-B) :type ,(getf opt :element-type))) + (setf sum (,(getf opt :f+) sum (,(getf opt :f*) A-val B-val)))) + :finally (,(getf opt :value-writer) (,(getf opt :f+) (,(getf opt :f*) alpha sum) val) sto-C of-C))))))))) + `(mlet* ,(recursive-append + '(((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))) + (when blas? + `((call-fortran? (> (max (nrows C) (ncols C) (if (eq job-A :n) (ncols A) (nrows A))) + ,fortran-lb-parameter)) + ((maj-A ld-A fop-A) (if call-fortran? (blas-matrix-compatible-p A job-A) (values nil 0 "?")) :type (symbol index-type (string 1))) + ((maj-B ld-B fop-B) (if call-fortran? (blas-matrix-compatible-p B job-B) (values nil 0 "?")) :type (symbol index-type (string 1))) + ((maj-C ld-C fop-C) (if call-fortran? (blas-matrix-compatible-p C :n) (values nil 0 "?")) :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))) ; - (rstp-B (row-stride B) :type index-type) - (cstp-B (col-stride B) :type index-type) - (hd-B (head B) :type index-type) - (sto-B (store B) :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) ; - (rstp-C (row-stride C) :type index-type) - (cstp-C (col-stride C) :type index-type) - (hd-C (head C) :type index-type) - (sto-C (store C) :type ,(linear-array-type (getf opt :store-type)))) - (when (eq job-A :t) - (rotatef rstp-A cstp-A)) - (when (eq job-B :t) - (rotatef rstp-B cstp-B)) - (very-quickly - (loop repeat nr-C - for rof-A of-type index-type = hd-A then (+ rof-A rstp-A) - for rof-C of-type index-type = hd-C then (+ rof-C rstp-C) - do (loop repeat nc-C - for cof-B of-type index-type = hd-B then (+ cof-B cstp-B) - for of-C of-type index-type = rof-C then (+ of-C cstp-C) - do (let-typed ((val (* beta ,(funcall (getf opt :reader) 'sto-C 'of-C)) :type ,(getf opt :element-type))) - (loop repeat dotl - for of-A of-type index-type = rof-A then (+ of-A cstp-A) - for of-B of-type index-type = cof-B then (+ of-B rstp-B) - summing (* ,(funcall (getf opt :reader) 'sto-A 'of-A) - ,(funcall (getf opt :reader) 'sto-B 'of-B)) into sum of-type ,(getf opt :element-type) - finally ,(funcall (getf opt :value-writer) '(+ (* alpha sum) val) 'sto-C 'of-C)))))))))) + (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))))) + (t + ,lisp-routine)) + lisp-routine))) C))) ;;Real (generate-typed-gemm! real-base-typed-gemm! - (real-matrix dgemm dgemv *real-l3-fcall-lb*)) + (real-tensor dgemm dgemv *real-l3-fcall-lb*)) (definline real-typed-gemm! (alpha A B beta C job) (real-base-typed-gemm! alpha A B beta C @@ -165,7 +176,7 @@ ;;Complex (generate-typed-gemm! complex-base-typed-gemm! - (complex-matrix zgemm zgemv *complex-l3-fcall-lb*)) + (complex-tensor zgemm zgemv *complex-l3-fcall-lb*)) (definline complex-typed-gemm! (alpha A B beta C job) (declare (type complex-matrix A B C) @@ -190,6 +201,11 @@ (complex-base-typed-gemm! alpha A B beta C tjob))))) +;;Symbolic +#+maxima +(generate-typed-gemm! symbolic-base-typed-gemm! + (symbolic-tensor nil nil 0)) + ;;---------------------------------------------------------------;; (defgeneric gemm! (alpha A B beta C &optional job) ----------------------------------------------------------------------- Summary of changes: matlisp.asd | 8 +- src/base/blas-helpers.lisp | 3 +- src/base/standard-tensor.lisp | 174 +++++++++++++++------------ src/classes/complex-tensor.lisp | 126 +++++++++++++------- src/classes/matrix.lisp | 4 +- src/classes/real-tensor.lisp | 113 ++++++++++++------ src/classes/symbolic-tensor.lisp | 126 +++++++++++++++++++ src/conditions.lisp | 7 - src/level-1/axpy.lisp | 82 ++++++++----- src/level-1/copy.lisp | 93 +++++++++------ src/level-1/dot.lisp | 128 ++++++++++---------- src/level-1/scal.lisp | 141 ++++++++++++++-------- src/level-1/swap.lisp | 40 ++++--- src/level-1/tensor-maker.lisp | 12 +- src/level-2/gemv.lisp | 106 +++++++++------- src/level-3/gemm.lisp | 246 ++++++++++++++++++++------------------ 16 files changed, 868 insertions(+), 541 deletions(-) create mode 100644 src/classes/symbolic-tensor.lisp hooks/post-receive -- matlisp |