Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-123-g6770dbf
|
From: Akshay S. <ak...@us...> - 2012-07-09 15:15:25
|
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 6770dbf44302c7d981ea50386827106748b8f3cc (commit)
via 38da10cc73eaa514e7bcacc1f996eeefda503078 (commit)
from cfd1ff7fa12112dcb0df038f9ecd60a5d637aa18 (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 6770dbf44302c7d981ea50386827106748b8f3cc
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Jul 9 20:40:26 2012 +0530
o Added gemv!
diff --git a/matlisp.asd b/matlisp.asd
index 0644b0e..26df9c3 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -117,7 +117,11 @@
(:file "realimag"
:depends-on ("copy"))
(:file "axpy"
- :depends-on ("copy"))))))
+ :depends-on ("copy"))))
+ (:module "matlisp-level-2"
+ :pathname "level-2"
+ :depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1")
+ :components ((:file "gemv")))))
;; (defclass f2cl-cl-source-file (asdf:cl-source-file)
diff --git a/packages.lisp b/packages.lisp
index 4ce172e..b90e858 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -34,6 +34,7 @@
;;Generic errors
#:generic-error #:message
#:invalid-type #:given #:expected
+ #:invalid-arguments
#:invalid-value #:given #:expected
#:unknown-token #:token
#:parser-error
diff --git a/src/base/blas-helpers.lisp b/src/base/blas-helpers.lisp
index 39de052..8e5ed08 100644
--- a/src/base/blas-helpers.lisp
+++ b/src/base/blas-helpers.lisp
@@ -36,18 +36,15 @@
unless (= so-st accumulated-off) do (return nil)
finally (return (aref sort-std 0))))))
-
-(defun blas-matrix-compatible-p (matrix &optional (op :n))
- (declare (type standard-tensor matrix))
- (let ((stds (strides matrix)))
- (declare (type (index-array *) stds))
- (if (not (= (array-dimension stds 0) 2)) nil
- (let ((rs (aref stds 0))
- (cs (aref stds 1)))
- (declare (type index-type rs cs))
- (cond
- ((= cs 1) (values :row-major rs (fortran-nop op)))
- ((= rs 1) (values :col-major cs (fortran-op op))))))))
+(defun blas-matrix-compatible-p (matrix op)
+ (declare (type standard-matrix matrix))
+ (let ((rs (aref (strides matrix) 0))
+ (cs (aref (strides matrix) 1)))
+ (declare (type index-type rs cs))
+ (cond
+ ((= cs 1) (values :row-major rs (fortran-nop op)))
+ ((= rs 1) (values :col-major cs (fortran-op op)))
+ (t (values nil 0 "?")))))
(definline fortran-op (op)
(ecase op (:n "N") (:t "T")))
diff --git a/src/conditions.lisp b/src/conditions.lisp
index a4e17a3..88c991f 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -32,6 +32,10 @@
(format stream "Given object of type ~A, expected ~A.~%" (given c) (expected c))
(call-next-method)))
+(defcondition invalid-arguments (generic-error)
+ ()
+ (:documentation "Given invalid arguments to the function."))
+
(defcondition invalid-value (generic-error)
((given-value :reader given :initarg :given)
(expected-value :reader expected :initarg :expected))
diff --git a/src/ffi/ffi-cffi.lisp b/src/ffi/ffi-cffi.lisp
index 344c8ec..8f57fa9 100644
--- a/src/ffi/ffi-cffi.lisp
+++ b/src/ffi/ffi-cffi.lisp
@@ -6,6 +6,8 @@
;; Callbacks : (:function <output-type> {(params)})
+;;TODO add declarations to generated wrappers.
+
(in-package #:matlisp-ffi)
(define-constant +ffi-types+ '(:single-float :double-float
diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp
index 53be2cd..edef2d3 100644
--- a/src/level-1/axpy.lisp
+++ b/src/level-1/axpy.lisp
@@ -77,8 +77,8 @@
is stored in Y and Y is returned.
")
(:method :before ((alpha number) (x standard-tensor) (y standard-tensor))
- (unless (idx= (dimensions x) (dimensions y))
- (error 'tensor-dimension-mismatch)))
+ (assert (idx= (dimensions x) (dimensions y)) nil
+ 'tensor-dimension-mismatch))
(:method ((alpha number) (x complex-tensor) (y real-tensor))
(error 'coercion-error :from 'complex-tensor :to 'real-tensor)))
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index 8651f95..fb708f0 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -1,80 +1,85 @@
-(in-package :matlisp)
+(in-package #:matlisp)
-(defmacro generate-typed-copy! (func (tensor-class blas-func))
+(defmacro generate-typed-gemv! (func
+ (matrix-class vector-class)
+ (blas-gemv-func blas-axpy-func blas-dot-func blas-scal-func))
;;Be very careful when using functions generated by this macro.
- ;;Indexes can be tricky and this has no safety net
+ ;;Indexes can be tricky and this has no safety net.
;;Use only after checking the arguments for compatibility.
- (let* ((opt (get-tensor-class-optimization tensor-class)))
- (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (from to)
- (declare (type ,tensor-class from to))
- (if-let (strd-p (blas-copyable-p from to))
- (,blas-func (number-of-elements from) (store from) (first strd-p) (store to) (second strd-p) (head from) (head to))
- (let ((f-sto (store from))
- (t-sto (store to)))
- (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto))
- (very-quickly
- ;;Can possibly make this faster (x2) by using ,blas-func in one of
- ;;the inner loops, but this is to me messy and as of now unnecessary.
- ;;SBCL can already achieve Fortran-ish speed inside this loop.
- (mod-dotimes (idx (dimensions from))
- with (linear-sums
- (f-of (strides from) (head from))
- (t-of (strides to) (head to)))
- do ,(funcall (getf opt :reader-writer) 'f-sto 'f-of 't-sto 't-of)))))
- to)))
-
-
-(defmacro generate-typed-gemv! (func (tensor-class blas-func))
- (let* ((opt (get-tensor-class-optimization tensor-class)))
- (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
+ (let* ((opt (get-tensor-class-optimization matrix-class)))
+ (assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class)
`(defun ,func (alpha A x beta y job)
- (declare (type (getf opt :element-type) alpha beta)
- (type ,tensor-class A x y)
- (type boolean job))
- (tensor-t
-
-
-;;There's no support for ":c", because there is no
-;;equivalent of ":n" with complex conjugation.
-(defmacro generate-typed-gemv!-func (func element-type store-type matrix-type blas-gemv-func blas-axpy-func blas-dot-func)
- ;;Be very careful when you use functions generated by this macro!
- ;;Indexes can be tricky and this has no safety net
- ;;Use only after checking the arguments for compatibility.
- `(defun ,func (alpha A x beta y job)
- (declare (tyep
- (declare (type ,element-type alpha beta)
- (type ,matrix-type A x y)
- (type symbol job))
- (mlet* (((st-a hd-a nr-a nc-a rs-a cs-a) (slot-values A '(store head number-of-rows number-of-cols row-stride col-stride))
- :type ((,store-type *) fixnum fixnum fixnum fixnum fixnum))
- ((st-x hd-x rs-x) (slot-values x '(store head row-stride))
- :type ((,store-type *) fixnum fixnum))
- ((st-y hd-y rs-y) (slot-values y '(store head row-stride))
- :type ((,store-type *) fixnum fixnum))
- ((sym lda tf-op) (blas-matrix-compatible-p A job) :type (symbol fixnum (string 1))))
- (if (not (string= tf-op "?"))
- (progn
- (when (eq sym :row-major)
- (rotatef nr-a nc-a)
- (rotatef rs-a cs-a))
- (,blas-gemv-func tf-op nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y))
- (progn
- (when (eq job :t)
- (rotatef nr-a nc-a)
- (rotatef rs-a cs-a))
- ;;Use the smaller of the loops.
- (if (> nr-a nc-a)
- (progn
- (scal! beta y)
- (dotimes (i nc-a)
- (,blas-axpy-func nr-a (* alpha (matrix-ref-2d x i 0)) st-a rs-a st-y rs-y :head-x (+ hd-a (* i cs-a)) :head-y hd-y)))
- (dotimes (i nr-a)
- (setf (matrix-ref-2d y i 0) (+ (* alpha (,blas-dot-func nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
- (* beta (matrix-ref-2d y i 0)))))))))
- y))
-
-;;
+ (declare (type ,(getf opt :element-type) alpha beta)
+ (type ,matrix-class A)
+ (type ,vector-class x y)
+ (type symbol job))
+ (mlet*
+ (((maj-a ld-a fop-a) (blas-matrix-compatible-p A job) :type (symbol index-type (string 1))))
+ (if maj-a
+ (let ((nr-a (aref (dimensions A) 0))
+ (nc-a (aref (dimensions A) 1)))
+ (declare (type index-type nr-a nc-a))
+ (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)))
+ (let ((nr-a (aref (dimensions A) 0))
+ (nc-a (aref (dimensions A) 1))
+ (rs-a (aref (strides A) 0))
+ (cs-a (aref (strides A) 1)))
+ (declare (type index-type nr-a nc-a rs-a cs-a))
+ (when (eq job :t)
+ (rotatef nr-a nc-a)
+ (rotatef rs-a cs-a))
+ (let ((sto-a (store a))
+ (sto-x (store x))
+ (std-x (aref (strides x) 0))
+ (hd-x (head x))
+ (sto-y (store y))
+ (std-y (aref (strides y) 0))
+ (hd-y (head y)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) sto-a sto-x sto-y)
+ (type index-type std-y std-x hd-x hd-y))
+ (if (> nr-a nc-a)
+ (progn
+ (unless (= beta 1d0)
+ (,blas-scal-func nr-a beta
+ sto-y std-y hd-y))
+ (very-quickly
+ (mod-dotimes (idx (idxv nc-a))
+ with (linear-sums
+ (of-x (strides x) (head x))
+ (of-a (idxv cs-a) (head A)))
+ do (,blas-axpy-func nr-a (* alpha ,(funcall (getf opt :reader) 'sto-x 'of-x))
+ sto-a rs-a sto-y std-y
+ of-a hd-y))))
+ (very-quickly
+ (mod-dotimes (idx (idxv nr-a))
+ with (linear-sums
+ (of-y (strides y) (head y))
+ (of-a (idxv rs-a) (head A)))
+ do (let ((val (* beta ,(funcall (getf opt :reader) 'sto-y 'of-y)))
+ (dotp (,blas-dot-func nc-a
+ sto-a cs-a sto-x std-x
+ of-a hd-x)))
+ (declare (type ,(getf opt :element-type) val dotp))
+ ,(funcall (getf opt :value-writer)
+ `(+ val (* alpha dotp)) 'sto-y 'of-y)))))))))
+ y)))
+
+(generate-typed-gemv! real-typed-gemv!
+ (real-matrix real-vector)
+ (dgemv daxpy ddot dscal))
+
+(generate-typed-gemv! complex-typed-gemv!
+ (complex-matrix complex-vector)
+ (zgemv zaxpy zdotu zscal))
+
+;;---------------------------------------------------------------;;
(defgeneric gemv! (alpha A x beta y &optional job)
(:documentation
"
@@ -100,184 +105,76 @@
---------------------------------------------------
:N (default) alpha * A * x + beta * y
:T alpha * A'* x + beta * y
-
- Note
- ====
- Take caution when using GEMM! as follows:
-
- (GEMV! alpha a x beta x)
-
- The results may be unpredictable depending
- on the underlying DGEMM, ZGEMM routines
- from BLAS, ATLAS or LIBCRUFT.
-"))
-
-(defmethod gemv! :before ((alpha number) (A standard-matrix) (x standard-matrix)
- (beta number) (y standard-matrix)
- &optional (job :n))
- (mlet* (((nr-a nc-a) (slot-values A '(number-of-rows number-of-cols)) :type (fixnum fixnum))
- ((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols)) :type (fixnum fixnum))
- ((nr-y nc-y) (slot-values y '(number-of-rows number-of-cols)) :type (fixnum fixnum)))
- (unless (member job '(:n :t))
- (error "Argument JOB to GEMV! is not recognized"))
- (when (eq job :t)
- (rotatef nr-a nc-a))
- (unless (and (= nc-x 1) (= nc-y 1)
- (= nc-a nr-x) (= nr-a nr-y))
- (error "Dimensions of A,x,y given to GEMV! do not match"))))
-
-;;
-(generate-typed-gemv!-func real-double-gemv!-typed
- double-float real-matrix-store-type real-matrix
- blas:dgemv blas:daxpy blas:ddot)
-
-(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
- (beta cl:real) (y real-matrix) &optional (job :n))
- ;; y <- \beta . y + \alpha . A o x
- (real-double-gemv!-typed (coerce alpha 'double-float) A x
- (coerce beta 'double-float) y job))
-
-;;
-(generate-typed-gemv!-func complex-double-gemv!-typed
- complex-double-float complex-matrix-store-type complex-matrix
- blas:zgemv blas:zaxpy blas:zdotu)
-
-(defmethod gemv! ((alpha number) (A complex-matrix) (x complex-matrix)
- (beta number) (y complex-matrix) &optional (job :n))
- ;; y <- \beta . y + \alpha . A o x
- (complex-double-gemv!-typed (complex-coerce alpha) A x
- (complex-coerce beta) y job))
-
-;
-(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
- (beta complex) (y complex-matrix) &optional (job :n))
- (let ((r-y (mrealpart~ y)))
- (declare (type real-matrix r-y))
- ;; y <- \beta * y
- (scal! (complex-coerce beta) y)
- ;; y <- y + \alpha * A o x
- (real-double-gemv!-typed (coerce alpha 'double-float) A x 1d0 r-y job)))
-
-(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
- (beta complex) (y complex-matrix) &optional (job :n))
- ;; y <- \beta * y
- (scal! (complex-coerce beta) y)
- ;; y <- y + \alpha * A o x
- (gemv! alpha A x 1d0 y job))
-
-(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
- (beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-be (coerce beta 'double-float))
- (r-al (coerce alpha 'double-float))
- (r-y (mrealpart~ y)))
- (declare (type double-float r-be r-al)
- (type real-matrix r-y))
- ;; y <- \beta * y
- (scal! r-be y)
- ;; (mrealpart~ y) <- (mrealpart~ y) + \alpha * A o x
- (real-double-gemv!-typed r-al A x 1d0 r-y job))
+")
+ (:method :before ((alpha number) (A standard-matrix) (x standard-vector)
+ (beta number) (y standard-vector)
+ &optional (job :n))
+ (assert (member job '(:n :t)) nil 'invalid-value
+ :given job :expected `(member job '(:n :t))
+ :message "Inside gemv!")
+ (assert (not (eq x y)) nil 'invalid-arguments
+ :message "GEMV!: x and y cannot be the same vector")
+ (assert (and
+ (= (aref (dimensions x) 0)
+ (aref (dimensions A) (if (eq job :t) 0 1)))
+ (= (aref (dimensions y) 0)
+ (aref (dimensions A) (if (eq job :t) 1 0))))
+ nil 'tensor-dimension-mismatch)))
+
+(defmethod gemv! ((alpha number) (A real-matrix) (x real-vector)
+ (beta number) (y real-vector) &optional (job :n))
+ (real-typed-gemv! (coerce-real alpha) A x
+ (coerce-real beta) y job))
+
+(defmethod gemv! ((alpha number) (A complex-matrix) (x complex-vector)
+ (beta number) (y complex-vector) &optional (job :n))
+ (complex-typed-gemv! (coerce-complex alpha) A x
+ (coerce-complex beta) y job))
+
+(defmethod gemv! ((alpha number) (A real-matrix) (x real-vector)
+ (beta number) (y complex-vector) &optional (job :n))
+ (unless (= beta 1)
+ (complex-typed-scal! (coerce-complex beta) y))
+ (unless (= alpha 0)
+ (if (complexp alpha)
+ (let ((A.x (make-real-tensor (aref (dimensions y) 0)))
+ (vw-y (tensor-realpart~ y)))
+ (real-typed-gemv! (coerce-real 1) A x (coerce-real 0) A.x job)
+ ;;
+ (real-typed-axpy! (coerce-real (realpart alpha)) A.x vw-y)
+ ;;Move view to the imaginary part
+ (incf (head vw-y))
+ (real-typed-axpy! (coerce-real (imagpart alpha)) A.x vw-y))
+ (real-typed-gemv! (coerce-real alpha) A x
+ (coerce-real 1) (tensor-realpart~ y) job)))
y)
-(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
- (beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-al (coerce (realpart alpha) 'double-float))
- (i-al (coerce (imagpart alpha) 'double-float))
- (r-be (coerce beta 'double-float))
- (r-y (mrealpart~ y))
- (i-y (mimagpart~ y)))
- (declare (type double-float r-al i-al r-be)
- (type real-matrix r-y i-y))
- ;; (mrealpart~ y) <- \beta * (mrealpart~ y) + (realpart \alpha) . A o x
- (real-double-gemv!-typed r-al A x r-be r-y job)
- ;; (mimagpart~ y) <- \beta * (mimagpart~ y) + (imagpart \alpha) . A o x
- (real-double-gemv!-typed i-al A x r-be i-y job))
- y)
-
-;
-(defmethod gemv! ((alpha number) (A real-matrix) (x complex-matrix)
- (beta complex) (y complex-matrix) &optional (job :n))
- ;; y <- \beta y
- (scal! beta y)
- ;; y <- y + \alpha . A o x
- (gemv! alpha A x 1d0 y job))
-
-(defmethod gemv! ((alpha cl:real) (A real-matrix) (x complex-matrix)
- (beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-x (mrealpart~ x))
- (i-x (mimagpart~ x))
- (r-y (mrealpart~ y))
- (i-y (mimagpart~ y))
- (r-al (coerce (realpart alpha) 'double-float))
- (r-be (coerce beta 'double-float)))
- (declare (type double-float r-al r-be)
- (type real-matrix r-x i-x r-y i-y))
- ;; (mrealpart~ y) <- \beta * (mrealpart~ y) + \alpha . A o (mrealpart~ x)
- (real-double-gemv!-typed r-al A r-x r-be r-y job)
- ;; (mimagpart~ y) <- \beta * (mimagpart~ y) + \alpha . A o (mrealpart~ x)
- (real-double-gemv!-typed r-al A i-x r-be i-y job))
- y)
-
-(defmethod gemv! ((alpha complex) (A real-matrix) (x complex-matrix)
- (beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-x (mrealpart~ x))
- (i-x (mimagpart~ x))
- (r-y (mrealpart~ y))
- (i-y (mimagpart~ y))
- (r-al (coerce (realpart alpha) 'double-float))
- (i-al (coerce (imagpart alpha) 'double-float))
- (r-be (coerce beta 'double-float)))
- (declare (type double-float r-al r-be i-al)
- (type real-matrix r-x i-x r-y i-y))
- (real-double-gemv!-typed r-al A r-x r-be r-y job)
- (real-double-gemv!-typed (- i-al) A i-x 1d0 r-y job)
- ;;
- (real-double-gemv!-typed i-al A r-x r-be i-y job)
- (real-double-gemv!-typed r-al A i-x 1d0 i-y job))
- y)
-
-;
-(defmethod gemv! ((alpha number) (A complex-matrix) (x real-matrix)
- (beta complex) (y complex-matrix) &optional (job :n))
- ;; y <- \beta y
- (scal! beta y)
- ;; y <- y + \alpha . A o x
- (gemv! alpha A x 1d0 y job))
-
-(defmethod gemv! ((alpha cl:real) (A complex-matrix) (x real-matrix)
- (beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-A (mrealpart~ A))
- (i-A (mimagpart~ A))
- (r-y (mrealpart~ y))
- (i-y (mimagpart~ y))
- (r-al (coerce (realpart alpha) 'double-float))
- (r-be (coerce beta 'double-float)))
- (declare (type double-float r-al r-be)
- (type real-matrix r-A i-A r-y i-y))
- ;; (mrealpart~ y) <- \beta * (mrealpart~ y) + \alpha . A o (mrealpart~ x)
- (real-double-gemv!-typed r-al r-A x r-be r-y job)
- ;; (mimagpart~ y) <- \beta * (mimagpart~ y) + \alpha . A o (mrealpart~ x)
- (real-double-gemv!-typed r-al i-A x r-be i-y job))
+(defmethod gemv! ((alpha number) (A real-matrix) (x complex-vector)
+ (beta number) (y complex-matrix) &optional (job :n))
+ (unless (= beta 1)
+ (complex-typed-scal! (coerce-complex beta) y))
+ (unless (= alpha 0)
+ (let ((A.x (make-complex-tensor (aref (dimensions y) 0))))
+ (let ((vw-x (tensor-realpart~ x))
+ (vw-A.x (tensor-realpart~ x)))
+ ;;Re
+ (real-typed-gemv! (coerce-real 1) A vw-x (coerce-real 0) vw-A.x job)
+ ;;Im
+ (incf (head vw-x))
+ (incf (head vw-A.x))
+ (real-typed-gemv! (coerce-real 1) A vw-x (coerce-real 0) vw-A.x job))
+ (complex-typed-axpy! (coerce-complex alpha) A.x y)))
y)
-(defmethod gemv! ((alpha complex) (A complex-matrix) (x real-matrix)
- (beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-A (mrealpart~ A))
- (i-A (mimagpart~ A))
- (r-y (mrealpart~ y))
- (i-y (mimagpart~ y))
- (r-al (coerce (realpart alpha) 'double-float))
- (i-al (coerce (imagpart alpha) 'double-float))
- (r-be (coerce beta 'double-float)))
- (declare (type double-float r-al r-be i-al)
- (type real-matrix r-A i-A r-y i-y))
- (real-double-gemv!-typed r-al r-A x r-be r-y job)
- (real-double-gemv!-typed (- i-al) i-A x 1d0 r-y job)
- ;;
- (real-double-gemv!-typed i-al r-A x r-be i-y job)
- (real-double-gemv!-typed r-al i-A x 1d0 i-y job))
+(defmethod gemv! ((alpha number) (A complex-matrix) (x real-vector)
+ (beta number) (y complex-vector) &optional (job :n))
+ (let ((cplx-x (make-complex-tensor (aref (dimensions x) 0))))
+ (real-typed-copy! x (tensor-realpart~ cplx-x))
+ (complex-typed-gemv! (coerce-complex alpha) A cplx-x
+ (coerce-complex beta) y job))
y)
-;;;;
+;;---------------------------------------------------------------;;
(defgeneric gemv (alpha A x beta y &optional job)
(:documentation
"
@@ -302,12 +199,17 @@
:T alpha * A'* x + beta * y
"))
-(defmethod gemv ((alpha number) (A standard-matrix) (x standard-matrix)
- (beta number) (y standard-matrix) &optional (job :n))
- (let ((result (scal (if (or (typep alpha 'complex) (typep beta 'complex)
- (typep A 'complex-matrix) (typep x 'complex-matrix))
- (complex-coerce beta)
- beta)
- y)))
- (declare (type standard-matrix y))
+(defmethod gemv ((alpha number) (A standard-matrix) (x standard-vector)
+ (beta number) (y complex-vector) &optional (job :n))
+ (let ((result (copy y)))
(gemv! alpha A x 1d0 result job)))
+
+(defmethod gemv ((alpha number) (A standard-matrix) (x standard-vector)
+ (beta number) (y real-vector) &optional (job :n))
+ (let ((result (if (or (complexp alpha) (complexp beta)
+ (typep A 'complex-matrix) (typep x'complex-matrix))
+ (make-complex-tensor (aref (dimensions y) 0))
+ (make-real-tensor (aref (dimensions y) 0)))))
+ (copy! y result)
+ (gemv! alpha A x beta result job)))
+
commit 38da10cc73eaa514e7bcacc1f996eeefda503078
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Jul 9 12:44:34 2012 +0530
o Added optimization for matrix, vector classes {follows to <>-tensor}.
o All methods in dot are now defined for <>-vector classes.
diff --git a/src/classes/complex-tensor.lisp b/src/classes/complex-tensor.lisp
index 3ba3f1e..828af4b 100644
--- a/src/classes/complex-tensor.lisp
+++ b/src/classes/complex-tensor.lisp
@@ -80,6 +80,8 @@ Cannot hold complex numbers."))
(rotatef (aref tstore (* 2 tidx)) (aref fstore (* 2 fidx)))
(rotatef (aref tstore (1+ (* 2 tidx))) (aref fstore (1+ (* 2 fidx)))))))
+(setf (get-tensor-class-optimization 'complex-matrix) 'complex-tensor
+ (get-tensor-class-optimization 'complex-vector) 'complex-tensor)
;;
(defmethod print-element ((tensor complex-tensor)
element stream)
diff --git a/src/classes/real-tensor.lisp b/src/classes/real-tensor.lisp
index 89bfefc..d896616 100644
--- a/src/classes/real-tensor.lisp
+++ b/src/classes/real-tensor.lisp
@@ -57,6 +57,9 @@ Allocates real storage. Default initial-element = 0d0.")
(lambda (fstore fidx tstore tidx)
(rotatef (aref tstore tidx) (aref fstore fidx))))
+(setf (get-tensor-class-optimization 'real-matrix) 'real-tensor
+ (get-tensor-class-optimization 'real-vector) 'real-tensor)
+
;;
(defmethod (setf tensor-ref) ((value number) (tensor real-tensor) subscripts)
(let ((sto-idx (store-indexing subscripts tensor)))
diff --git a/src/level-1/dot.lisp b/src/level-1/dot.lisp
index 6e7a18a..59d00ab 100644
--- a/src/level-1/dot.lisp
+++ b/src/level-1/dot.lisp
@@ -58,33 +58,25 @@
If X and Y are both scalars then this is the same
as (* (CONJUGATE X) Y) if CONJUAGTE-P and (* X Y)
otherwise.
-"))
+")
+ (:method :before ((x standard-vector) (y standard-vector) &optional (conjugate-p t))
+ (declare (ignore conjugate-p))
+ (unless (idx= (dimensions x) (dimensions y))
+ (error 'tensor-dimension-mismatch))))
(defmethod dot ((x number) (y number) &optional (conjugate-p t))
(if conjugate-p
(* (conjugate x) y)
(* x y)))
-(defmethod dot :before ((x standard-tensor) (y standard-tensor) &optional (conjugate-p t))
- (declare (ignore conjugate-p))
- (unless (and (vector-p x) (vector-p y))
- (error 'tensor-not-vector
- :rank (cond
- ((not (vector-p x))
- (rank x))
- ((not (vector-p y))
- (rank y)))))
- (unless (idx= (dimensions x) (dimensions y))
- (error 'tensor-dimension-mismatch)))
-
-(defmethod dot ((x real-tensor) (y real-tensor) &optional (conjugate-p t))
+(defmethod dot ((x real-vector) (y real-vector) &optional (conjugate-p t))
(declare (ignore conjugate-p))
(ddot (number-of-elements x)
(store x) (aref (strides x) 0)
(store y) (aref (strides y) 0)
(head x) (head y)))
-(defmethod dot ((x real-tensor) (y complex-tensor) &optional (conjugate-p t))
+(defmethod dot ((x real-vector) (y complex-vector) &optional (conjugate-p t))
(declare (ignore conjugate-p))
(let ((nele (number-of-elements x))
(std-x (aref (strides x) 0))
@@ -99,13 +91,13 @@
rpart
(complex rpart ipart)))))
-(defmethod dot ((x complex-tensor) (y real-tensor) &optional (conjugate-p t))
+(defmethod dot ((x complex-vector) (y real-vector) &optional (conjugate-p t))
(let ((cres (dot y x)))
(if conjugate-p
(conjugate cres)
cres)))
-(defmethod dot ((x complex-tensor) (y complex-tensor) &optional (conjugate-p t))
+(defmethod dot ((x complex-vector) (y complex-vector) &optional (conjugate-p t))
(let ((nele (number-of-elements x))
(std-x (aref (strides x) 0))
(hd-x (head x))
-----------------------------------------------------------------------
Summary of changes:
matlisp.asd | 6 +-
packages.lisp | 1 +
src/base/blas-helpers.lisp | 21 +--
src/classes/complex-tensor.lisp | 2 +
src/classes/real-tensor.lisp | 3 +
src/conditions.lisp | 4 +
src/ffi/ffi-cffi.lisp | 2 +
src/level-1/axpy.lisp | 4 +-
src/level-1/dot.lisp | 26 +--
src/level-2/gemv.lisp | 408 +++++++++++++++------------------------
10 files changed, 192 insertions(+), 285 deletions(-)
hooks/post-receive
--
matlisp
|
