Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch classy updated. 2012-02-23-246-g808353d
|
From: Akshay S. <ak...@us...> - 2013-08-15 02:20:05
|
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, classy has been updated
via 808353d428ddc07d365bf1de8abcc86f0179ee08 (commit)
from f9bf6a61b1860942b520069596563b2db546f927 (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 808353d428ddc07d365bf1de8abcc86f0179ee08
Author: Akshay Srinivasan <aks...@gm...>
Date: Wed Aug 14 19:20:10 2013 -0700
Ported gemm.
diff --git a/matlisp.asd b/matlisp.asd
index 8450621..3f30b7d 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -153,7 +153,6 @@
: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" "matlisp-level-2")
diff --git a/packages.lisp b/packages.lisp
index 5c7413f..04096db 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -68,7 +68,7 @@
(defpackage "MATLISP-UTILITIES"
(:use #:common-lisp #:matlisp-conditions)
- (:export #:ensure-list #:id
+ (:export #:ensure-list #:id #:ieql
#:vectorify #:copy-n
#:ensure-args #:repsym #:findsym #:find-tag
#:zip #:zip-eq #:zipsym
@@ -92,7 +92,6 @@
#:inlining #:definline
#:with-optimization #:quickly #:very-quickly #:slowly #:quickly-if))
-
(defpackage "MATLISP-TEMPLATE"
(:use #:common-lisp #:matlisp-utilities)
(:export #:deft/generic #:deft/method #:remt/method))
diff --git a/src/base/blas-helpers.lisp b/src/base/blas-helpers.lisp
index de623a8..dd3fdaa 100644
--- a/src/base/blas-helpers.lisp
+++ b/src/base/blas-helpers.lisp
@@ -24,12 +24,12 @@
(list csto-a? csto-b?)))))
(definline fortran-nop (op)
- (ecase op (#\t #\n) (#\n #\t)))
+ (ecase op (#\T #\N) (#\N #\T)))
(defun split-job (job)
(declare (type symbol job))
(let-typed ((name (symbol-name job) :type string))
- (loop :for x :across name :collect (char-downcase x))))
+ (loop :for x :across name :collect (char-upcase x))))
(definline flip-major (job)
(declare (type symbol job))
@@ -46,7 +46,7 @@
(cs (aref stds 1) :type index-type))
;;Note that it is not required that (rs = nc * cs) or (cs = nr * rs)
(cond
- ((= cs 1) (values rs (fortran-nop op) :row-major))
+ ((and (char/= op #\C) (= cs 1)) (values rs (fortran-nop op) :row-major))
((= rs 1) (values cs op :col-major)))))
;;Stride makers.
@@ -77,6 +77,29 @@
(defun make-stride (dims)
(ecase *default-stride-ordering* (:row-major (make-stride-rmj dims)) (:col-major (make-stride-cmj dims))))
-(defun call-fortran? (x lb)
+(defun call-fortran? ( x lb)
(declare (type standard-tensor x))
(> (size x) lb))
+
+(defmacro with-columnification ((type (&rest input) (&rest output)) &rest body)
+ (with-gensyms (cfunc)
+ (let ((input-syms (mapcar #'(lambda (x)
+ (assert (or (symbolp (second x)) (characterp (second x))) nil "Given a non-symbolic input.")
+ (gensym (symbol-name (car x)))) input))
+ (output-syms (mapcar #'(lambda (mat) (gensym (symbol-name mat))) output)))
+ `(labels ((,cfunc (a &optional b)
+ (declare (type ,type a))
+ (let ((ret (or b (let ((*default-stride-ordering* :col-major)) (t/zeros ,type (the index-store-vector (dimensions a)))))))
+ (declare (type ,type a ret))
+ (t/copy! (,type ,type) a ret))))
+ (let (,@(mapcar #'(lambda (x sym) (let ((mat (first x)) (job (second x)))
+ `(,sym (if (blas-matrix-compatiblep ,mat ,job) ,mat
+ (,cfunc ,mat))))) input input-syms)
+ ,@(mapcar #'(lambda (mat sym) `(,sym (if (eql (third (multiple-value-list (blas-matrix-compatiblep ,mat #\N))) :col-major) ,mat
+ (,cfunc ,mat)))) output output-syms))
+ (declare (type ,type ,@(append input-syms output-syms)))
+ (symbol-macrolet (,@(mapcar #'(lambda (mat sym) `(,mat ,sym)) (append (mapcar #'car input) output) (append input-syms output-syms)))
+ ,@body)
+ ,@(mapcar #'(lambda (mat sym) `(unless (eql (third (multiple-value-list (blas-matrix-compatiblep ,mat #\N))) :col-major)
+ (,cfunc ,sym ,mat))) output output-syms)
+ nil)))))
diff --git a/src/base/template.lisp b/src/base/template.lisp
index f027c83..50a2d83 100644
--- a/src/base/template.lisp
+++ b/src/base/template.lisp
@@ -33,6 +33,7 @@
(declare (type ,ty ,@args))
(cl:/ ,@args))))
+;;
(deft/generic (t/fc #'subtypep) ty (num))
(deft/method t/fc (ty number) (num)
(with-gensyms (num-sym)
@@ -51,7 +52,9 @@
`(defmethod fconj ((x ,clname))
(t/fc ,clname x)))
(fc x))))
-
+(defun field-realp (fil)
+ (eql (macroexpand-1 `(t/fc ,fil phi)) 'phi))
+;;
(deft/generic (t/f= #'subtypep) ty (&rest nums))
(deft/method t/f= (ty number) (&rest nums)
(let* ((decl (zipsym nums))
diff --git a/src/classes/numeric.lisp b/src/classes/numeric.lisp
index 8605abb..93893f5 100644
--- a/src/classes/numeric.lisp
+++ b/src/classes/numeric.lisp
@@ -3,6 +3,12 @@
(defclass numeric-tensor (standard-tensor) ())
(deft/method t/field-type (sym numeric-tensor) ()
'number)
+
+;;
+(defleaf integer-tensor (numeric-tensor) ())
+(deft/method t/field-type (sym integer-tensor) ()
+ 'integer)
+
;;
(defclass blas-numeric-tensor (numeric-tensor) ())
(deft/generic (t/l1-lb #'subtypep) sym ())
diff --git a/src/level-1/maker.lisp b/src/level-1/maker.lisp
index 46a3f76..cf41512 100644
--- a/src/level-1/maker.lisp
+++ b/src/level-1/maker.lisp
@@ -4,7 +4,9 @@
(deft/method t/zeros (class standard-tensor) (dims &optional initial-element)
(with-gensyms (astrs adims sizs)
`(let* ((,adims (make-index-store ,dims)))
+ (declare (type index-store-vector ,adims))
(multiple-value-bind (,astrs ,sizs) (make-stride ,adims)
+ (declare (type index-store-vector ,astrs))
(make-instance ',class
:dimensions ,adims
:head 0
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index 177d7c8..f793ccd 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -39,12 +39,14 @@
(type character ,transp))
(unless (t/f= ,(field-type sym) ,beta (t/fid* ,(field-type sym)))
(t/scdi! ,sym ,beta ,y :scal? t :numx? t))
+ ,@(when (field-realp (field-type sym))
+ `((when (char= ,transp #\C) (setq ,transp #\T))))
;;These loops are optimized for column major matrices
- (ecase (char-upcase ,transp)
+ (ecase ,transp
(#\N (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (ref ,A i j) (ref ,x j)) nil))
(#\T (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (ref ,A j i) (ref ,x j)) nil))
- (#\C (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A i j)) (ref ,x j)) nil))
- (#\H (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A j i)) (ref ,x j)) nil)))
+ ,@(unless (field-realp (field-type sym))
+ `((#\C (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A j i)) (ref ,x j)) nil)))))
,y)))
;;---------------------------------------------------------------;;
@@ -72,48 +74,29 @@
JOB Operation
---------------------------------------------------
:N (default) alpha * A * x + beta * y
- :T alpha * transpose(A)* x + beta * y
- :C alpha * conjugate(A) * x + beta * y
- :H alpha * transpose o conjugate(A) + beta * y
+ :T alpha * A^T * x + beta * y
+ :C alpha * A^H * x + beta * y
")
(:method :before (alpha (A standard-tensor) (x standard-tensor)
beta (y standard-tensor)
&optional (job :n))
- (assert (member job '(:n :t :c :h)) nil 'invalid-value
- :given job :expected `(member job '(:n :t :c :h))
+ (assert (member job '(:n :t :c)) nil 'invalid-value
+ :given job :expected `(member job '(:n :t :c))
:message "GEMV!: Given an unknown job.")
(assert (not (eq x y)) nil 'invalid-arguments
:message "GEMV!: x and y cannot be the same vector")
(assert (and
(tensor-vectorp x) (tensor-vectorp y) (tensor-matrixp A)
(= (aref (the index-store-vector (dimensions x)) 0)
- (aref (the index-store-vector (dimensions A)) (if (member job '(:t :h)) 0 1)))
+ (aref (the index-store-vector (dimensions A)) (if (member job '(:t :c)) 0 1)))
(= (aref (the index-store-vector (dimensions y)) 0)
- (aref (the index-store-vector (dimensions A)) (if (member job '(:t :h)) 1 0))))
+ (aref (the index-store-vector (dimensions A)) (if (member job '(:t :c)) 1 0))))
nil 'tensor-dimension-mismatch)))
-(defun ieql (&rest args)
- (loop :for ele :in (cdr args)
- :do (unless (eql (car args) ele)
- (return nil))
- :finally (return t)))
-
-(defun >class (cls)
- (let ((clist (reverse (sort cls #'coerceable?))))
- (loop :for ele :in (cdr clist)
- :do (unless (coerceable? ele (car clist))
- (return nil))
- :finally (return (car clist)))))
-
-(defun tensor-coerce (ten cls &optional (duplicate? t))
- (let ((clname (if (typep cls 'standard-class) (class-name cls) cls)))
- (if (and (not duplicate?) (typep ten clname)) ten
- (copy! ten (zeros (dimensions ten) clname)))))
-
(defmethod gemv! (alpha (A standard-tensor) (x standard-tensor) beta (y standard-tensor) &optional (job :n))
(let ((clx (class-name (class-of x)))
(cly (class-name (class-of y)))
- (cla (class-name (class-of y))))
+ (cla (class-name (class-of A))))
(assert (and (member cla *tensor-type-leaves*)
(member clx *tensor-type-leaves*)
(member cly *tensor-type-leaves*))
@@ -126,7 +109,7 @@
(beta (t/coerce ,(field-type clx) beta))
(cjob (aref (symbol-name job) 0)))
(declare (type ,(field-type clx) alpha beta)
- (type character trans-op))
+ (type character cjob))
,(recursive-append
(when (subtypep clx 'blas-numeric-tensor)
`(if (call-fortran? A (t/l2-lb ,cla))
@@ -144,7 +127,7 @@
y))
(gemv! alpha A x beta y job))
(t
- (error "Don't know how to apply axpy! to classes ~a, ~a." clx cly)))))
+ (error "Don't know how to apply gemv! to classes ~a." (list cla clx cly))))))
;;---------------------------------------------------------------;;
(defgeneric gemv (alpha A x beta y &optional job)
(:documentation
@@ -167,9 +150,8 @@
JOB Operation
---------------------------------------------------
:N (default) alpha * A * x + beta * y
- :T alpha * A'* x + beta * y
- :C alpha * conjugate(A) * x + beta * y
- :H alpha * transpose o conjugate(A) + beta * y
+ :T alpha * A^T * x + beta * y
+ :C alpha * A^H * x + beta * y
"))
(defmethod gemv (alpha (A standard-tensor) (x standard-tensor)
diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp
index 776406f..7596787 100644
--- a/src/level-3/gemm.lisp
+++ b/src/level-3/gemm.lisp
@@ -14,7 +14,7 @@
`(let* (,@decl
(,m (aref (the index-store-vector (dimensions ,C)) 0))
(,n (aref (the index-store-vector (dimensions ,C)) 1))
- (,k (aref (the index-store-vector (dimensions ,A)) (ecase (char-upcase ,transa) ((#\N #\C) 1) ((#\T #\H) 0)))))
+ (,k (aref (the index-store-vector (dimensions ,A)) (ecase (char-upcase ,transa) (#\N 1) ((#\T #\C) 0)))))
(declare (type ,sym ,A ,B ,C)
(type ,(field-type sym) ,alpha ,beta)
(type index-type ,lda ,ldb ,ldc ,m ,n ,k)
@@ -33,69 +33,45 @@
;;
(deft/generic (t/gemm! #'subtypep) sym (alpha A B beta C transa transb))
-;;Witness the power of macros, muggles! :)
(deft/method t/gemm! (sym standard-tensor) (alpha A B beta C transa transb)
- (using-gensyms (decl (alpha A x beta y transp))
+ (using-gensyms (decl (alpha A B beta C transa transb))
`(let (,@decl)
- (declare (type ,sym ,A ,x ,y)
+ (declare (type ,sym ,A ,B ,C)
(type ,(field-type sym) ,alpha ,beta)
- (type character ,transp))
+ (type character ,transa ,transb))
(unless (t/f= ,(field-type sym) ,beta (t/fid* ,(field-type sym)))
- (t/scdi! ,sym ,beta ,y :scal? t :numx? t))
+ (t/scdi! ,sym ,beta ,C :scal? t :numx? t))
+ ,@(when (field-realp (field-type sym))
+ `((when (char= ,transa #\C) (setq ,transa #\T))
+ (when (char= ,transb #\C) (setq ,transb #\T))))
;;These loops are optimized for column major matrices
- (ecase (char-upcase ,transp)
- (#\N (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (ref ,A i j) (ref ,x j)) nil))
- (#\T (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (ref ,A j i) (ref ,x j)) nil))
- (#\C (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A i j)) (ref ,x j)) nil))
- (#\H (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A j i)) (ref ,x j)) nil)))
- ,y)))
-
-
-;;Real
-(generate-typed-gemm! real-base-typed-gemm!
- (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
- (apply #'combine-jobs
- (mapcar #'(lambda (x)
- (ecase x ((:n :t) x) (:h :t) (:c :n)))
- (multiple-value-list (split-job job))))))
-
-;;Complex
-(generate-typed-gemm! complex-base-typed-gemm!
- (complex-tensor zgemm *complex-l3-fcall-lb*))
-
-(definline complex-typed-gemm! (alpha A B beta C job)
- (declare (type complex-matrix A B C)
- (type complex-type alpha beta)
- (type symbol job))
- (multiple-value-bind (job-A job-B) (split-job job)
- (if (and (member job-A '(:n :t))
- (member job-B '(:n :t)))
- (complex-base-typed-gemm! alpha A B beta C job)
- (let ((A (ecase job-A
- ((:n :t) A)
- ((:h :c) (let ((ret (complex-typed-copy! A (complex-typed-zeros (dimensions A)))))
- (real-typed-num-scal! -1d0 (tensor-imagpart~ ret))
- ret))))
- (B (ecase job-B
- ((:n :t) B)
- ((:h :c) (let ((ret (complex-typed-copy! A (complex-typed-zeros (dimensions A)))))
- (real-typed-num-scal! -1d0 (tensor-imagpart~ ret))
- ret))))
- (tjob (combine-jobs (ecase job-A ((:n :t) job-A) (:h :t) (:c :n))
- (ecase job-B ((:n :t) job-B) (:h :t) (:c :n)))))
- (complex-base-typed-gemm! alpha A B
- beta C tjob)))))
-
-;;Symbolic
-#+maxima
-(generate-typed-gemm! symbolic-base-typed-gemm!
- (symbolic-tensor nil 0))
-
+ ,(labels ((transpose-ref (mat)
+ `(ref ,(cadr mat) ,@(reverse (cddr mat))))
+ (conjugate-ref (mat)
+ `(t/fc ,(field-type sym) ,mat))
+ (generate-mm-code (transa transb)
+ (destructuring-bind (A-ref B-ref) (mapcar #'(lambda (mat trans) (ecase trans
+ ((#\N #\T) mat)
+ ((#\C) (conjugate-ref mat))))
+ (mapcar #'(lambda (mat trans) (ecase trans
+ ((#\N) mat)
+ ((#\T #\C) (transpose-ref mat))))
+ (list `(ref ,A i k) `(ref ,B k j)) (list transa transb))
+ (list transa transb))
+ (let ((loopo (let ((ta (member transa '(#\T #\C)))
+ (tb (member transb '(#\T #\C))))
+ (cond
+ ((and (not ta) (not tb)) `(j k i))
+ ((and (not ta) tb) `(k j i))
+ (t`(i j k))))))
+ `(einstein-sum ,sym ,loopo (ref ,C i j) (* ,alpha ,A-ref ,B-ref) nil)))))
+ `(ecase ,transa
+ ,@(loop :for ta :across (if (field-realp (field-type sym)) "NT" "NTC")
+ :collect `(,ta (ecase ,transb
+ ,@(loop :for tb :across (if (field-realp (field-type sym)) "NT" "NTC")
+ :collect `(,tb ,(generate-mm-code ta tb))))))))
+ ,C)))
;;---------------------------------------------------------------;;
-
(defgeneric gemm! (alpha A B beta C &optional job)
(:documentation
"
@@ -118,20 +94,10 @@
JOB must be a keyword with two of these alphabets
N Identity
T Transpose
- C Complex conjugate
- H Hermitian transpose {conjugate transpose}
-
- so that (there are 4x4 operations in total).
-
- JOB Operation
- ---------------------------------------------------
- :NN (default) alpha * A * B + beta * C
- :TN alpha * transpose(A) * B + beta * C
- :NH alpha * A * transpose o conjugate(B) + beta * C
- :HC alpha * transpose o conjugate(A) * conjugate(B) + beta * C
+ C Hermitian transpose {conjugate transpose}
")
- (:method :before ((alpha number) (A standard-matrix) (B standard-matrix)
- (beta number) (C standard-matrix)
+ (:method :before (alpha (A standard-tensor) (B standard-tensor)
+ beta (C standard-tensor)
&optional (job :nn))
(let ((nr-a (nrows A))
(nc-a (ncols A))
@@ -140,64 +106,46 @@
(nr-c (nrows C))
(nc-c (ncols C)))
(declare (type index-type nr-a nc-a nr-b nc-b nr-c nc-c))
- (let ((sjobs (multiple-value-list (split-job job))))
+ (let ((sjobs (split-job job)))
(assert (= (length sjobs) 2) nil 'invalid-arguments :message "Ill formed job")
- (ecase (first sjobs) ((:n :c) t) ((:t :h) (rotatef nr-a nc-a)))
- (ecase (second sjobs) ((:n :c) t) ((:t :h) (rotatef nr-b nc-b))))
+ (ecase (first sjobs) (#\N t) ((#\T #\C) (rotatef nr-a nc-a)))
+ (ecase (second sjobs) ((#\N) t) ((#\T #\C) (rotatef nr-b nc-b))))
(assert (not (or (eq A C) (eq B C))) nil 'invalid-arguments
:message "GEMM!: C = {A or B} is not allowed.")
- (assert (and (= nr-c nr-a)
+ (assert (and (tensor-matrixp A) (tensor-matrixp B) (tensor-matrixp C)
+ (= nr-c nr-a)
(= nc-a nr-b)
(= nc-b nc-c)) nil 'tensor-dimension-mismatch))))
-
-(defmethod gemm! ((alpha number) (a real-matrix) (b real-matrix)
- (beta number) (c real-matrix)
- &optional (job :nn))
- (real-typed-gemm! (coerce-real alpha) a b
- (coerce-real beta) c job))
-
-(defmethod gemm! ((alpha number) (a complex-matrix) (b complex-matrix)
- (beta number) (c complex-matrix)
- &optional (job :nn))
- (complex-typed-gemm! (coerce-complex alpha) a b
- (coerce-complex beta) c job))
-
-(defmethod gemm! ((alpha number) (a real-matrix) (b real-matrix)
- (beta number) (c complex-matrix)
- &optional (job :nn))
- (unless (= beta 1)
- (scal! beta c))
- (unless (= alpha 0)
- (if (complexp alpha)
- (let ((A.x (make-real-tensor (nrows c) (ncols c)))
- (vw.c (tensor-realpart~ c)))
- (real-typed-gemm! (coerce-real 1) A B (coerce-real 0) A.x job)
- ;;Re
- (axpy! (realpart alpha) A.x vw.c)
- ;;Im
- (incf (head vw.c))
- (axpy! (imagpart alpha) A.x vw.c))
- (let ((vw.c (tensor-realpart~ c)))
- (real-typed-gemm! (coerce-real alpha) A B
- (coerce-real 1) vw.c job))))
- C)
-
-(defmethod gemm! ((alpha number) (a real-matrix) (b complex-matrix)
- (beta number) (c complex-matrix)
- &optional (job :nn))
- (let ((A.cplx (copy! A (make-complex-tensor (nrows a) (ncols a)))))
- (complex-typed-gemm! (coerce-complex alpha) A.cplx B
- (coerce-complex beta) C job))
- C)
-
-(defmethod gemm! ((alpha number) (a complex-matrix) (b real-matrix)
- (beta number) (c complex-matrix)
- &optional (job :nn))
- (let ((B.cplx (copy! B (make-complex-tensor (nrows B) (ncols B)))))
- (complex-typed-gemm! (coerce-complex alpha) A B.cplx
- (coerce-complex beta) C job))
- C)
-
+
+(defmethod gemm! (alpha (A standard-tensor) (B standard-tensor) beta (C standard-tensor) &optional (job :nn))
+ (let ((cla (class-name (class-of A)))
+ (clb (class-name (class-of B)))
+ (clc (class-name (class-of C))))
+ (assert (and (member cla *tensor-type-leaves*)
+ (member clb *tensor-type-leaves*)
+ (member clc *tensor-type-leaves*))
+ nil 'tensor-abstract-class :tensor-class (list cla clb clc))
+ (cond
+ ((ieql cla clb clc)
+ (compile-and-eval
+ `(defmethod gemm! (alpha (A ,cla) (B ,clb) beta (C ,clc) &optional (job :nn))
+ (let ((alpha (t/coerce ,(field-type cla) alpha))
+ (beta (t/coerce ,(field-type cla) beta)))
+ (declare (type ,(field-type cla) alpha beta))
+ (destructuring-bind (joba jobb) (split-job job)
+ (declare (type character joba jobb))
+ ,(recursive-append
+ (when (subtypep clc 'blas-numeric-tensor)
+ `(if (call-fortran? C (t/l3-lb ,clc))
+ (with-columnification (,cla ((a joba) (b jobb)) (c))
+ (multiple-value-bind (lda opa) (blas-matrix-compatiblep a joba)
+ (multiple-value-bind (ldb opb) (blas-matrix-compatiblep b jobb)
+ (t/blas-gemm! ,cla alpha A lda B ldb beta C ldb opa opb))))))
+ `(t/gemm! ,cla alpha A B beta C joba jobb))))
+ C))
+ (gemm! alpha A B beta C job))
+ (t
+ (error "Don't know how to apply gemm! to classes ~a." (list cla clb clc))))))
;;---------------------------------------------------------------;;
(defgeneric gemm (alpha a b beta c &optional job)
(:documentation
@@ -221,46 +169,9 @@
JOB must be a keyword with two of these alphabets
N Identity
T Transpose
- C Complex conjugate
- H Hermitian transpose {conjugate transpose}
-
- so that (there are 4x4 operations in total).
-
- JOB Operation
- ---------------------------------------------------
- :NN (default) alpha * A * B + beta * C
- :TN alpha * transpose(A) * B + beta * C
- :NH alpha * A * transpose o conjugate(B) + beta * C
- :HC alpha * transpose o conjugate(A) * conjugate(B) + beta * C
+ C Hermitian conjugate
"))
-(defmethod gemm ((alpha number) (a standard-matrix) (b standard-matrix)
- (beta number) (c complex-matrix)
- &optional (job :nn))
- (let ((result (copy C)))
- (gemm! alpha A B beta result job)))
-
-;; if all args are not real then at least one of them
-;; is complex, so we need to call GEMM! with a complex C
-(defmethod gemm ((alpha number) (a standard-matrix) (b standard-matrix)
- (beta number) (c real-matrix)
- &optional (job :nn))
- (let ((result (funcall (if (or (complexp alpha) (complexp beta)
- (typep a 'complex-matrix) (typep b 'complex-matrix))
- #'complex-typed-zeros
- #'real-typed-zeros)
- (dimensions C))))
- (copy! C result)
- (gemm! alpha A B beta result job)))
-
-(defmethod gemm ((alpha number) (a standard-matrix) (b standard-matrix)
- (beta (eql nil)) (c (eql nil))
- &optional (job :nn))
- (multiple-value-bind (job-A job-B) (split-job job)
- (let ((result (funcall (if (or (complexp alpha) (complexp beta)
- (typep a 'complex-matrix) (typep b 'complex-matrix))
- #'complex-typed-zeros
- #'real-typed-zeros)
- (make-index-store (list (if (member job-A '(:n :c)) (nrows A) (ncols A))
- (if (member job-B '(:n :c)) (ncols B) (nrows B)))))))
- (gemm! alpha A B 0 result job))))
+(defmethod gemm (alpha (A standard-tensor) (B standard-tensor)
+ beta (C standard-tensor) &optional (job :nn))
+ (gemm! alpha A B beta (copy C) job))
diff --git a/src/utilities/functions.lisp b/src/utilities/functions.lisp
index a224404..f790208 100644
--- a/src/utilities/functions.lisp
+++ b/src/utilities/functions.lisp
@@ -6,6 +6,12 @@
(declaim (inline id))
(defun id (x) x)
+(defun ieql (&rest args)
+ (loop :for ele :in (cdr args)
+ :do (unless (eql (car args) ele)
+ (return nil))
+ :finally (return t)))
+
(declaim (inline vectorify))
(defun vectorify (seq n &optional (element-type t))
(declare (type (or vector list) seq))
-----------------------------------------------------------------------
Summary of changes:
matlisp.asd | 1 -
packages.lisp | 3 +-
src/base/blas-helpers.lisp | 31 +++++-
src/base/template.lisp | 5 +-
src/classes/numeric.lisp | 6 +
src/level-1/maker.lisp | 2 +
src/level-2/gemv.lisp | 50 +++------
src/level-3/gemm.lisp | 241 +++++++++++++-----------------------------
src/utilities/functions.lisp | 6 +
9 files changed, 138 insertions(+), 207 deletions(-)
hooks/post-receive
--
matlisp
|
