[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-141-g9a94775

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 422-6466

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  9a94775cd4eb5593fea88f5cf665bcadc198fb6f (commit)
      from  cd30ca81e687388cf532e30e08f79b68cf56c325 (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 9a94775cd4eb5593fea88f5cf665bcadc198fb6f
Author: Akshay Srinivasan <aks...@gm...>
Date:   Tue Jul 24 13:56:44 2012 +0530

    o Added (:h, :c) job handling ability to gemm!

diff --git a/src/base/blas-helpers.lisp b/src/base/blas-helpers.lisp
index 339f417..1e3fcb6 100644
--- a/src/base/blas-helpers.lisp
+++ b/src/base/blas-helpers.lisp
@@ -57,3 +57,12 @@
     ((string= sop "N") "T")
     ((string= sop "T") "N")
     (t (error "Unrecognised fortran-op."))))
+
+(defun split-job (job)
+  (values-list
+   (map 'list #'(lambda (x) (intern (string x) "KEYWORD")) (symbol-name job))))
+
+(defun combine-jobs (&rest jobs)
+  (let ((job (intern (apply #'concatenate 'string (mapcar #'symbol-name jobs)) "KEYWORD")))
+    job))
+  
diff --git a/src/level-1/trans.lisp b/src/level-1/trans.lisp
index 56bad3e..d5c0087 100644
--- a/src/level-1/trans.lisp
+++ b/src/level-1/trans.lisp
@@ -153,7 +153,7 @@
   =======
   Like mconjugate!, but non-destructive."
   (etypecase A
-    (standard-tensor (copy (mconjugate! A)))
+    (standard-tensor (mconjugate! (copy A)))
     (number (conjugate A))))
 
 ;;
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index d4dfadd..78004a5 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -72,7 +72,9 @@
 
 (defun complex-typed-gemv! (alpha A x beta y job)
   (declare (type complex-matrix A)
-	   (type complex-vector x y))
+	   (type complex-vector x y)
+	   (type complex-type alpha beta)
+	   (type symbol job))
   (if (member job '(:n :t))
       (complex-base-typed-gemv! alpha A x beta y job)
       ;;The CBLAS way.
@@ -83,8 +85,6 @@
 
 ;;---------------------------------------------------------------;;
 
-;;Can't support "C" because its dual (complex-conjugate without transpose)
-;;isn't supported by BLAS (which'd be needed for row-major matrices).
 (defgeneric gemv! (alpha A x beta y &optional job)
   (:documentation
 "
diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp
index f5172b7..325378c 100644
--- a/src/level-3/gemm.lisp
+++ b/src/level-3/gemm.lisp
@@ -152,31 +152,36 @@
 							  finally ,(funcall (getf opt :value-writer) '(+ (* alpha sum) val) 'sto-C 'of-C)))))))))))
        C)))
 
-;;Tweakable
-(defparameter *real-gemm-fortran-call-lower-bound* 100
-  "
-  If the maximum dimension in the MM is lower than this
-  parameter, then the lisp code is used by default, instead of
-  calling BLAS. Used to avoid the FFI overhead when calling
-  MM with small matrices.
-  Default set with SBCL on x86-64 linux. A reasonable value
-  is something between 20 and 200.")
-(generate-typed-gemm! real-typed-gemm! (real-matrix
-					dgemm dgemv
-					*real-gemm-fortran-call-lower-bound*))
+;;Real
+(generate-typed-gemm! real-base-typed-gemm!
+  (real-matrix dgemm dgemv *real-l3-fcall-lb*))
+
+(definline real-typed-gemm! (alpha A B beta C job)
+  (real-base-typed-gemv! 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-matrix zgemm zgemv *complex-l3-fcall-lb*))
+
+(defun 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 ((:h :c) (mconjugate A)) ((:n :t) A)))
+	      (B (ecase job-B ((:h :c) (mconjugate B)) ((:n :t) B)))
+	      (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)))))
 
-;;Tweakable
-(defparameter *complex-gemm-fortran-call-lower-bound* 60
-  "
-  If the maximum dimension in the MM is lower than this
-  parameter, then the lisp code is used by default, instead of
-  calling BLAS. Used to avoid the FFI overhead when calling
-  MM with small matrices.
-  Default set with SBCL on x86-64 linux. A reasonable value
-  is something between 20 and 200.")
-(generate-typed-gemm! complex-typed-gemm! (complex-matrix
-					   zgemm zgemv
-					   *complex-gemm-fortran-call-lower-bound*))
 ;;---------------------------------------------------------------;;
 
 (defgeneric gemm! (alpha A B beta C &optional job)
@@ -198,12 +203,20 @@
   alpha,beta are scalars and A,B,C are matrices.
   op(A) means either A or A'.
 
+  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 * A'* B + beta * C
-     :NT                alpha * A * B'+ beta * C
-     :TT                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
 ")
   (:method :before ((alpha number) (A standard-matrix) (B standard-matrix)
 		    (beta number) (C standard-matrix)
@@ -215,12 +228,10 @@
 	(nr-c (nrows C))
 	(nc-c (ncols C)))
     (declare (type index-type nr-a nc-a nr-b nc-b nr-c nc-c))
-    (case job
-      (:nn t)
-      (:tn (rotatef nr-a nc-a))
-      (:nt (rotatef nr-b nc-b))
-      (:tt (rotatef nr-a nc-a) (rotatef nr-b nc-b))
-      (t (error 'invalid-value :given job :expected '(member job '(:nn :tn :nt :tt)))))
+    (let ((sjobs (multiple-value-list (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))))
     (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)
@@ -295,12 +306,20 @@
   alpha,beta are scalars and A,B,C are matrices.
   op(A) means either A or A'.
 
+  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 * A'* B + beta * C
-     :NT                alpha * A * B'+ beta * C
-     :TT                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
 "))
 
 (defmethod gemm ((alpha number) (a standard-matrix) (b standard-matrix)
@@ -324,16 +343,12 @@
 (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) (ecase job
-				       (:nn (values :n :n))
-				       (:nt (values :n :t))
-				       (:tn (values :t :n))
-				       (:tt (values :t :t)))
+  (multiple-value-bind (job-A job-B) (split-job job)
     (let ((result (apply
 		   (if (or (complexp alpha) (complexp beta)
 			   (typep a 'complex-matrix) (typep b 'complex-matrix))
 		       #'make-complex-tensor
 		       #'make-real-tensor)
-		   (list (if (eq job-A :n) (nrows A) (ncols A))
-			 (if (eq job-B :n) (ncols B) (nrows B))))))
+		   (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))))

-----------------------------------------------------------------------

Summary of changes:
 src/base/blas-helpers.lisp |    9 ++++
 src/level-1/trans.lisp     |    2 +-
 src/level-2/gemv.lisp      |    6 +-
 src/level-3/gemm.lisp      |  101 +++++++++++++++++++++++++-------------------
 4 files changed, 71 insertions(+), 47 deletions(-)


hooks/post-receive
-- 
matlisp


