Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-124-g83d4618
|
From: Akshay S. <ak...@us...> - 2012-07-10 07:38: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 83d461878939dc99bbe82269d113041fdbdc9e52 (commit)
from 6770dbf44302c7d981ea50386827106748b8f3cc (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 83d461878939dc99bbe82269d113041fdbdc9e52
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Jul 10 13:04:32 2012 +0530
o Added gemm generating macro.
diff --git a/README b/README
index 029550d..3bec6c9 100644
--- a/README
+++ b/README
@@ -29,6 +29,7 @@ This is the development branch of Matlisp.
* QUADPACK: Move from f2cl-ed version to the Fortran one.
* MINPACK: Move from f2cl-ed version to the Fortran one.
* ODEPACK: Add abstraction for DLSODE, and DLSODAR may others too.
+ * Add a Lisp generic wrapper for every BLAS func {low priority}.
*** Syntactic sugar
* Add array slicing macros
diff --git a/packages.lisp b/packages.lisp
index b90e858..190bca9 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -74,7 +74,7 @@
#:list-dimensions
;;Macros
#:when-let #:if-let #:if-ret #:with-gensyms #:let-rec
- #:mlet* #:make-array-allocator
+ #:mlet* #:make-array-allocator #:let-typed
#:nconsc #:define-constant
#:macrofy
;;
diff --git a/src/classes/matrix.lisp b/src/classes/matrix.lisp
index ad059aa..df19960 100644
--- a/src/classes/matrix.lisp
+++ b/src/classes/matrix.lisp
@@ -1,8 +1,5 @@
(in-package #:matlisp)
-;;
-
-
(definline nrows (matrix)
(declare (type standard-matrix matrix))
(aref (dimensions matrix) 0))
@@ -26,9 +23,6 @@
(list (aref dims 0) (aref dims 1))))
;;
-
-
-;;
(definline row-matrix-p (matrix)
"
Syntax
@@ -62,32 +56,32 @@
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)
+(definline square-matrix-p (matrix)
(and (square-p matrix) (matrix-p matrix)))
-;;
-(defgeneric fill-matrix (matrix fill-element)
- (:documentation
- "
- Syntax
- ======
- (FILL-MATRIX matrix fill-element)
+;; ;;
+;; (defgeneric fill-matrix (matrix fill-element)
+;; (:documentation
+;; "
+;; Syntax
+;; ======
+;; (FILL-MATRIX matrix fill-element)
- Purpose
- =======
- Fill MATRIX with 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"))
+;; (defmethod fill-matrix ((matrix t) (fill t))
+;; (error "arguments MATRIX and FILL to FILL-MATRIX must be a
+;; matrix and a number"))
-;;
-;;
+;; ;;
+;; ;;
-;;
+;; ;;
-(definline matrix-ref (matrix row &optional col)
- (declare (type standard-matrix matrix))
- (tensor-ref matrix `(,row ,col)))
+;; (definline matrix-ref (matrix row &optional col)
+;; (declare (type standard-matrix matrix))
+;; (tensor-ref matrix `(,row ,col)))
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index fb708f0..29a4d64 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -80,6 +80,8 @@
(zgemv zaxpy zdotu zscal))
;;---------------------------------------------------------------;;
+
+;;Can't support "C" because the dual isn't supported by BLAS.
(defgeneric gemv! (alpha A x beta y &optional job)
(:documentation
"
diff --git a/src/old/gemm.lisp b/src/level-3/gemm.lisp
similarity index 62%
rename from src/old/gemm.lisp
rename to src/level-3/gemm.lisp
index 5ca6f21..c7d2e00 100644
--- a/src/old/gemm.lisp
+++ b/src/level-3/gemm.lisp
@@ -25,138 +25,130 @@
;;; ENHANCEMENTS, OR MODIFICATIONS.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Originally written by Raymond Toy.
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; $Id: gemm.lisp,v 1.8 2011/01/25 18:36:56 rtoy Exp $
-;;;
-;;; $Log: gemm.lisp,v $
-;;; Revision 1.8 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.7.2.1 2011/01/25 18:16:53 rtoy
-;;; Use cl:real instead of real.
-;;;
-;;; Revision 1.7 2004/05/24 16:34:22 rtoy
-;;; More SBCL support from Robert Sedgewick. The previous SBCL support
-;;; was incomplete.
-;;;
-;;; Revision 1.6 2001/06/22 12:52:41 rtoy
-;;; Use ALLOCATE-REAL-STORE and ALLOCATE-COMPLEX-STORE to allocate space
-;;; instead of using the error-prone make-array.
-;;;
-;;; Revision 1.5 2001/02/26 17:44:54 rtoy
-;;; Remove the complex-alpha,beta special variables. (Make a closure out
-;;; of them.)
-;;;
-;;; Revision 1.4 2000/07/11 18:02:03 simsek
-;;; o Added credits
-;;;
-;;; Revision 1.3 2000/07/11 02:11:56 simsek
-;;; o Added support for Allegro CL
-;;;
-;;; Revision 1.2 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.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")
-
-(defmacro generate-typed-gemm!-func (func element-type store-type matrix-type blas-gemm-func lisp-gemv-func)
- `(defun ,func (alpha a b beta c job)
- (declare (optimize (safety 0) (speed 3))
- (type ,element-type alpha beta)
- (type ,matrix-type a b c)
- (type symbol job))
- (mlet* ((job-a (ecase job ((:nn :nt) :n) ((:tn :tt) :t)) :type symbol)
- (job-b (ecase job ((:nn :tn) :n) ((:nt :tt) :t)) :type symbol)
- ((hd-c nr-c nc-c st-c) (slot-values c '(head number-of-rows number-of-cols store))
- :type (fixnum fixnum fixnum (,store-type *)))
- ((hd-a st-a) (slot-values a '(head store))
- :type (fixnum (,store-type *)))
- ((hd-b st-b) (slot-values b '(head store))
- :type (fixnum (,store-type *)))
- (k (if (eq job-a :n)
- (ncols a)
- (nrows a))
- :type fixnum)
- ((order-a lda fort-job-a) (blas-matrix-compatible-p a job-a)
- :type (symbol fixnum (string 1)))
- ((order-b ldb fort-job-b) (blas-matrix-compatible-p b job-b)
- :type (symbol fixnum (string 1)))
- ((order-c ldc fort-job-c) (blas-matrix-compatible-p c :n)
- :type (nil fixnum (string 1))))
- ;;
- (if (and (> lda 0) (> ldb 0) (> ldc 0))
- (progn
- (when (string= fort-job-c "T")
- (rotatef a b)
- (rotatef lda ldb)
- (rotatef fort-job-a fort-job-b)
- (rotatef hd-a hd-b)
- (rotatef st-a st-b)
- (rotatef nr-c nc-c)
- ;;
- (setf fort-job-a (fortran-snop fort-job-a))
- (setf fort-job-b (fortran-snop fort-job-b)))
- (,blas-gemm-func fort-job-a fort-job-b
- nr-c nc-c k
- alpha
- st-a lda
- st-b ldb
- beta
- st-c ldc
- :head-a hd-a :head-b hd-b :head-c hd-c))
- (progn
- (when (eq job-a :t) (transpose! a))
- (when (eq job-b :t) (transpose! b))
- ;;
- (symbol-macrolet
- ((loop-col
- (mlet* ((cs-b (col-stride b) :type fixnum)
- (cs-c (col-stride c) :type fixnum)
- (col-b (col~ b 0) :type ,matrix-type)
- (col-c (col~ c 0) :type ,matrix-type))
- (dotimes (j nc-c)
- (when (> j 0)
- (setf (head col-b) (+ (head col-b) cs-b))
- (setf (head col-c) (+ (head col-c) cs-c)))
- (,lisp-gemv-func alpha a col-b beta col-c :n))))
- (loop-row
- (mlet* ((rs-a (row-stride a) :type fixnum)
- (rs-c (row-stride c) :type fixnum)
- (row-a (transpose! (row~ a 0)) :type ,matrix-type)
- (row-c (transpose! (row~ c 0)) :type ,matrix-type))
- (dotimes (i nr-c)
- (when (> i 0)
- (setf (head row-a) (+ (head row-a) rs-a))
- (setf (head row-c) (+ (head row-c) rs-c)))
- (,lisp-gemv-func alpha b row-a beta row-c :t)))))
- (cond
- (order-a loop-col)
- (order-b loop-row)
- ((< nr-c nc-c) loop-row)
- (t loop-col)))
- ;;
- (when (eq job-a :t) (transpose! a))
- (when (eq job-b :t) (transpose! b))
- )))
- c))
+(in-package #:matlisp)
+
+(defmacro generate-typed-gemm! (func (matrix-class) (blas-gemm-func blas-gemv-func))
+ (let* ((opt (get-tensor-class-optimization matrix-class)))
+ (assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class)
+ `(defun ,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))
+ ((maj-a ld-a fop-a) (blas-matrix-compatible-p A job-a) :type (symbol index-type (string 1)))
+ ((maj-b ld-b fop-b) (blas-matrix-compatible-p B job-b) :type (symbol index-type (string 1)))
+ ((maj-c ld-c fop-c) (blas-matrix-compatible-p C :n) :type (symbol index-type nil)))
+ (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)))
+ (cond
+ (maj-a
+ (let-typed ((nc-c (ncols C) :type index-type)
+ (sto-c (store C) :type ,(linear-array-type (getf opt :store-type)))
+ (stp-c (row-stride C) :type index-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)
+ (strd-c (col-stride C) :type index-type))
+ (when (eq maj-a :row-major)
+ (rotatef nr-a nc-a))
+ (very-quickly
+ (mod-dotimes (idx (idxv nc-c))
+ with (linear-sums
+ (of-b (idxv strd-b) (head B))
+ (of-c (idxv strd-c) (head 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)))))
+ (maj-b
+ (let-typed ((nr-c (nrows C) :type index-type)
+ (stp-c (col-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 B) (row-stride B)) :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)))
+ (strd-a (if (eq job-A :n) (row-stride A) (col-stride A)) :type index-type)
+ (strd-c (row-stride C) :type index-type))
+ (when (eq maj-b :row-major)
+ (rotatef nr-b nc-b))
+ (very-quickly
+ (mod-dotimes (idx (idxv nr-c))
+ with (linear-sums
+ (of-A (idxv strd-a) (head A))
+ (of-c (idxv strd-c) (head 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 ((dotl (ecase job-a (:n (ncols A)) (:t (nrows A))) :type index-type)
+ (rstp-a (row-stride A) :type index-type)
+ (cstp-a (col-stride A) :type index-type)
+ (rstp-b (row-stride A) :type index-type)
+ (cstp-b (col-stride A) :type index-type)
+ (sto-a (store A) :type ,(linear-array-type (getf opt :store-type)))
+ (sto-b (store B) :type ,(linear-array-type (getf opt :store-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
+ (mod-dotimes (idx (dimensions C))
+ with (loop-order :row-major)
+ with (linear-sums
+ (of-a (idxv rstp-a 0) (head A)) ; cstp-a))
+ (of-b (idxv 0 cstp-b) (head B)) ; rstp-b))
+ (of-c (strides C) (head C))) ; 0)))
+ do (let-typed ((tmp (,(getf opt :coercer) 0) :type ,(getf opt :element-type))
+ (val (* beta ,(funcall (getf opt :reader) 'sto-c 'of-c)) :type ,(getf opt :element-type)))
+ (loop repeat dotl
+ for dof-a of-type index-type = of-a then (+ dof-a cstp-a)
+ for dof-b of-type index-type = of-b then (+ dof-b rstp-b)
+ do (incf tmp (* ,(funcall (getf opt :reader) 'sto-a 'dof-a)
+ ,(funcall (getf opt :reader) 'sto-b 'dof-b))))
+ ,(funcall (getf opt :value-writer) '(+ (* alpha tmp) (* beta val)) 'sto-c 'of-c)))))))))
+ C)))
+
+(generate-typed-gemm! real-typed-gemm! (real-matrix) (dgemm dgemv))
+(generate-typed-gemm! complex-typed-gemm! (complex-matrix) (zgemm zgemv))
+
+(let ((A (tensor-realpart~
+ (make-complex-tensor '((1 2 3)
+ (4 5 6)
+ (7 8 9)))))
+ (C (make-real-tensor 3 3)))
+ (real-typed-gemm! 1d0 A A 0d0 C :nn))
+
;;;;
(defgeneric gemm! (alpha a b beta c &optional job)
(:documentation
@@ -406,4 +398,4 @@
(complex-coerce beta)
beta)
c)))
- (gemm! alpha a b 1d0 result job)))
\ No newline at end of file
+ (gemm! alpha a b 1d0 result job)))
diff --git a/src/utilities.lisp b/src/utilities.lisp
index 4fbdeab..e2e74c4 100644
--- a/src/utilities.lisp
+++ b/src/utilities.lisp
@@ -1,5 +1,6 @@
(in-package #:matlisp-utilities)
+;;TODO: cleanup!
(defmacro mlet* (decls &rest body)
"
mlet* ({ {(var*) | var} values-form &keyform declare type}*) form*
@@ -55,6 +56,49 @@
`(progn
,@body))))
+(defmacro let-typed (bindings &rest body)
+"
+ let-typed ({var form &key type}*) form*
+
+ This macro also handles type declarations.
+
+ Example:
+ > (let-typed ((x 1 :type fixnum))
+ (+ 1 x))
+
+ expands into:
+
+ > (let ((x 1))
+ (declare (type fixnum x))
+ (+ 1 x))
+"
+ (labels ((parse-bindings (bdng let-decl type-decl)
+ (if (null bdng) (values (reverse let-decl) (reverse type-decl))
+ ;;Unless the user gives a initialisation form, no point declaring type
+ ;; {var is bound to nil}.
+ (destructuring-bind (var &optional form &key (type nil)) (ensure-list (car bdng))
+ (parse-bindings (cdr bdng)
+ (cons (if form `(,var ,form) var) let-decl)
+ (if type
+ (cons `(type ,type ,var) type-decl)
+ type-decl))))))
+ (multiple-value-bind (let-bdng type-decl) (parse-bindings bindings nil nil)
+ (let ((decl-code (recursive-append
+ (cond
+ ((and (consp (first body))
+ (eq (caar body) 'declare))
+ (first body))
+ ((consp type-decl)
+ '(declare ))
+ (t nil))
+ type-decl)))
+ `(let (,@let-bdng)
+ ,@(if (null decl-code) nil `(,decl-code))
+ ,@(if (and (consp (first body))
+ (eq (caar body) 'declare))
+ (cdr body)
+ body))))))
+
(defmacro let-rec (name arglist &rest code)
"
(let-rec name ({var [init-form]}*) declaration* form*) => result*
diff --git a/tests/loopy-tests.lisp b/tests/loopy-tests.lisp
index 9f63e62..de334ec 100644
--- a/tests/loopy-tests.lisp
+++ b/tests/loopy-tests.lisp
@@ -58,7 +58,36 @@
(of-a (idxv n 1 0))
(of-b (idxv 0 n 1))
(of-c (idxv n 0 1)))
- do (incf (aref st-c of-c) (* (aref st-a of-a) (aref st-b of-b)))))))))
+ do (incf (aref st-c of-c) (* (aref st-a of-a) (aref st-b of-b)))))))))
+
+(defun test-mm-daxpy (n)
+ (let* ((t-a (make-real-tensor n n))
+ (t-b (make-real-tensor n n))
+ (t-c (make-real-tensor n n))
+ (st-a (store t-a))
+ (st-b (store t-b))
+ (st-c (store t-c)))
+ (declare (type real-tensor t-a t-b t-c)
+ (type (real-array *) st-a st-b st-c))
+ (mod-dotimes (idx (dimensions t-a))
+ with (linear-sums
+ (of-a (strides t-a))
+ (of-b (strides t-b))
+ (of-c (strides t-c)))
+ do (setf (aref st-a of-a) (random 1d0)
+ (aref st-b of-b) (random 1d0)
+ (aref st-c of-c) 0d0))
+ (time
+ (very-quickly
+ (mod-dotimes (idx (idxv n n))
+ with (loop-order :row-major)
+ with (linear-sums
+ (of-a (idxv 1 0))
+ (of-b (idxv n 1))
+ (of-c (idxv 1 0)))
+ do (daxpy n (aref st-b of-b) st-a n st-c n
+ of-a of-c))))))
+
(defun test-mm-ddot (n)
(let* ((t-a (make-real-tensor n n))
-----------------------------------------------------------------------
Summary of changes:
README | 1 +
packages.lisp | 2 +-
src/classes/matrix.lisp | 48 ++++----
src/level-2/gemv.lisp | 2 +
src/{old => level-3}/gemm.lisp | 256 +++++++++++++++++++---------------------
src/utilities.lisp | 44 +++++++
tests/loopy-tests.lisp | 31 +++++-
7 files changed, 223 insertions(+), 161 deletions(-)
rename src/{old => level-3}/gemm.lisp (62%)
hooks/post-receive
--
matlisp
|
