From: Akshay S. <ak...@us...> - 2012-03-14 05:58:28
|
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, matlisp-cffi has been updated via 517949e31b41b303dab91670e80b207bc45d3256 (commit) via 452bd980d066101fad780815ff5ddfccdcf5e683 (commit) via f4b94df81eab2264c071a5a6f592e6becf76f770 (commit) via 8a2722a7452cac1de3edd5fc8cc65d97c0d1a0a8 (commit) via ecab7fc4f75f615eabed198e0e7a325b2c055fd0 (commit) via 14707c10ced9ee08431313323cc1905ac77e2a21 (commit) via 7f45afca17c2501c0b2a777809e1998683e5a17f (commit) via 8708b5e4f149f2ba16b879bfa87aa3085097cb34 (commit) via 4d23b62fbced07e732614e290a08e90a731942ab (commit) via 579b17087f685d198095ef8616a6374ed59d6848 (commit) via ec16280b56855cff043ec1fc6bbddc11cf2cebca (commit) via 6695b66525322e6817de899b467fb3505ec22775 (commit) via b4da181a45007d71d24371851089844b257e2fa8 (commit) from 2cbbcb64c997e7ed95c2e9a344347ca9606fce66 (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 517949e31b41b303dab91670e80b207bc45d3256 Author: Akshay Srinivasan <aks...@gm...> Date: Wed Mar 14 11:24:44 2012 +0530 -> gemv core function works. diff --git a/src/copy.lisp b/src/copy.lisp index cf987d9..03c6e8e 100644 --- a/src/copy.lisp +++ b/src/copy.lisp @@ -80,6 +80,10 @@ ;; (defmacro generate-typed-copy!-func (func store-type matrix-type blas-func) + ;;Be very careful when using functions generated by this macro. + ;;Indexes can be tricky and this has no safety net + ;;(you don't see a matrix-ref do you ?) + ;;Use only after checking the arguments for compatibility. `(defun ,func (mat-a mat-b) (declare (type ,matrix-type mat-a mat-b) (optimize (safety 0) (speed 3))) @@ -89,24 +93,26 @@ ((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *)))) (if (and cp-a cp-b) (,blas-func sz st-a inc-a st-b inc-b :head-x hd-a :head-y hd-b) - (symbol-macrolet - ((common-code - (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride)) - :type (fixnum fixnum fixnum fixnum)) - ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride)) - :type (fixnum fixnum))) - (loop for i from 0 below nr-a - do (,blas-func nc-a st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b))))))) - ;;Choose the smaller of the loops - (if (> (nrows mat-a) (ncols mat-a)) - (with-transpose! (mat-a mat-b) - common-code) - common-code))) - mat-b))) + (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride)) + :type (fixnum fixnum fixnum fixnum)) + ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride)) + :type (fixnum fixnum))) + ;;Choose the smaller of the loops + (when (> (nrows mat-a) (ncols mat-a)) + (rotatef nr-a nc-a) + (rotatef rs-a cs-a) + (rotatef rs-b cs-b)) + (loop for i from 0 below nr-a + do (,blas-func nc-a st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b))))))) + mat-b)) ;; (defmacro generate-typed-num-copy!-func (func element-type store-type matrix-type blas-func array-decl) + ;;Be very careful when using functions generated by this macro. + ;;Indexes can be tricky and this has no safety net + ;;(you don't see a matrix-ref do you ?) + ;;Use only after checking the arguments for compatibility. (let ((num-arg (car array-decl))) (destructuring-bind (var var-maker-form var-setf-form &key type) (cadr array-decl) `(mlet* (((,var) (,@var-maker-form) ,@(if type @@ -120,18 +126,15 @@ ((hd-x st-x sz) (slot-values mat-x '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum))) (if cp-x (,blas-func sz ,var 0 st-x inc-x :head-y hd-x) - (symbol-macrolet - ((common-code - (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride)) - :type (fixnum fixnum fixnum fixnum))) - (loop for i from 0 below nr-x - do (,blas-func nc-x ,var 0 st-x cs-x :head-y (+ hd-x (* i rs-x))))))) - ;;Choose the smaller of the loops - (if (> (nrows mat-x) (ncols mat-x)) - (with-transpose! (mat-x) - common-code) - common-code))) - mat-x)))))) + (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride)) + :type (fixnum fixnum fixnum fixnum))) + ;;Choose the smaller of the loops + (when (> (nrows mat-x) (ncols mat-x)) + (rotatef nr-x nc-x) + (rotatef rs-x cs-x)) + (loop for i from 0 below nr-x + do (,blas-func nc-x ,var 0 st-x cs-x :head-y (+ hd-x (* i rs-x))))))) + mat-x))))) ;; (defgeneric copy! (matrix new-matrix) diff --git a/src/gemv.lisp b/src/gemv.lisp index 9adb6cc..cffa7d6 100644 --- a/src/gemv.lisp +++ b/src/gemv.lisp @@ -1,35 +1,32 @@ (in-package :matlisp) (defun real-double-gemv!-typed (alpha A x beta y job) - (declare ((type double-float alpha beta) - (type real-matrix A x y) - (type symbol job))) - (symbol-macrolet - ((common-block - (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 ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum)) - ((st-x hd-x rs-x) (slot-values x '(store head row-stride)) - :type ((real-matrix-store-type *) fixnum fixnum)) - ((st-y hd-y rs-y) (slot-values y '(store head row-stride)) - :type ((real-matrix-store-type *) fixnum fixnum)) - ((sym lda tf-job) (blas-matrix-compatible-p A) :type (nil fixnum (string 1)))) - (if (not (string= tf-job "?")) - (blas:dgemv tf-job 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) - (dotimes (i nr-a) - (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot 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)))) - - (mlet* ((fortran-job (ecase job (:n "N") (:t "T")) :type ((string 1))) + ;;Be very careful when you use this function! + ;;Indexes can be tricky and this has no safety net + ;;(you don't see a matrix-ref do you ?) + ;;Use only after checking the arguments for compatibility. + (declare (type double-float alpha beta) + (type real-matrix A x y) + (type symbol job)) + (mlet* ((fort-op (ecase job (:n "N") (:t "T")) :type ((string 1))) ((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 ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum)) ((st-x hd-x rs-x) (slot-values x '(store head row-stride)) :type ((real-matrix-store-type *) fixnum fixnum)) ((st-y hd-y rs-y) (slot-values y '(store head row-stride)) :type ((real-matrix-store-type *) fixnum fixnum)) - ((sym lda tf-job) (blas-matrix-compatible-p A fortran-job) :type (nil fixnum (string 1)))) - (if (not (string= tf-job "?")) - (blas:dgemv tf-job 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) - (dotimes (i nr-a) - (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot 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)))) - \ No newline at end of file + ((sym lda tf-op) (blas-matrix-compatible-p A fort-op) :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:dgemv 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 (string= fort-op "T") + (rotatef nr-a nc-a) + (rotatef rs-a cs-a)) + (dotimes (i nr-a) + (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot 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) \ No newline at end of file diff --git a/src/scal.lisp b/src/scal.lisp index 88da573..a81af20 100644 --- a/src/scal.lisp +++ b/src/scal.lisp @@ -70,6 +70,10 @@ (in-package "MATLISP") (defmacro generate-typed-scal!-func (func element-type store-type matrix-type blas-func) + ;;Be very careful when using functions generated by this macro. + ;;Indexes can be tricky and this has no safety net + ;;(you don't see a matrix-ref do you ?) + ;;Use only after checking the arguments for compatibility. `(defun ,func (alpha mat-x) (declare (type ,matrix-type mat-x) (type ,element-type alpha) @@ -80,17 +84,15 @@ :type (fixnum (,store-type *)))) (if cp-x (,blas-func sz-x alpha st-x inc-x :head-x hd-x) - (symbol-macrolet - ((common-code - (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride)) - :type (fixnum fixnum fixnum fixnum))) - (loop for i from 0 below nr-x - do (,blas-func nc-x alpha st-x cs-x :head-x (+ hd-x (* i rs-x))))))) - (if (> (nrows mat-x) (ncols mat-x)) - (with-transpose! (mat-x) - common-code) - common-code))) - mat-x))) + (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride)) + :type (fixnum fixnum fixnum fixnum))) + ;;Choose the smaller of the loops. + (when (> (nrows mat-x) (ncols mat-x)) + (rotatef nr-x nc-x) + (rotatef rs-x cs-x)) + (loop for i from 0 below nr-x + do (,blas-func nc-x alpha st-x cs-x :head-x (+ hd-x (* i rs-x))))))) + mat-x)) ;; (defgeneric scal! (alpha x) diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp index 3c2cfce..0f1ec71 100644 --- a/src/standard-matrix.lisp +++ b/src/standard-matrix.lisp @@ -318,8 +318,8 @@ matrix and a number")) ;; (defun blas-matrix-compatible-p (matrix &optional (fortran-op "N")) (declare (optimize (safety 0) (speed 3)) - (type (or real-matrix complex-matrix) mat)) - (mlet* (((rs cs) (slot-values mat '(row-stride col-stride)) + (type (or real-matrix complex-matrix) matrix)) + (mlet* (((rs cs) (slot-values matrix '(row-stride col-stride)) :type (fixnum fixnum))) (cond ((= cs 1) (values :row-major rs (cond commit 452bd980d066101fad780815ff5ddfccdcf5e683 Author: Akshay Srinivasan <aks...@gm...> Date: Tue Mar 13 22:03:18 2012 +0530 -> Ported axpy to support new classes. Started work on Level-2 stuff. Does not work yet. diff --git a/src/axpy.lisp b/src/axpy.lisp index 57104cf..d6cee37 100644 --- a/src/axpy.lisp +++ b/src/axpy.lisp @@ -74,23 +74,78 @@ (in-package "MATLISP") -#+nil (use-package "LAPACK") -#+nil (use-package "BLAS") -#+nil (use-package "FORTRAN-FFI-ACCESSORS") +(defmacro generate-typed-axpy!-func (func element-type store-type matrix-type blas-func) + `(defun ,func (alpha mat-a mat-b) + (declare (type ,element-type alpha) + (type ,matrix-type mat-a mat-b) + (optimize (safety 0) (speed 3))) + (mlet* (((cp-a inc-a sz-a) (blas-copyable-p mat-a) :type (boolean fixnum nil)) + ((cp-b inc-b sz-b) (blas-copyable-p mat-b) :type (boolean fixnum nil)) + ((hd-a st-a sz) (slot-values mat-a '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum)) + ((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *)))) + (if (and cp-a cp-b) + (,blas-func sz alpha st-a inc-a st-b inc-b :head-x hd-a :head-y hd-b) + (symbol-macrolet + ((common-code + (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride)) + :type (fixnum fixnum fixnum fixnum)) + ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride)) + :type (fixnum fixnum))) + (loop for i from 0 below nr-a + do (,blas-func nc-a alpha st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b))))))) + ;;Choose the smaller of the loops + (if (> (nrows mat-a) (ncols mat-a)) + (with-transpose! (mat-a mat-b) + common-code) + common-code))) + mat-b))) -#+nil (export '(axpy! - axpy)) +;; +(defgeneric axpy! (alpha x y) + (:documentation + " + Syntax + ====== + (AXPY! alpha x y) + + Y <- alpha * x + y + + Purpose + ======= + Same as AXPY except that the result + is stored in Y and Y is returned. +")) + +(defmethod axpy! :before ((alpha number) (x standard-matrix) (y standard-matrix)) + (mlet* (((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 (and (= nr-x nr-y) (= nc-x nc-y)) + (error "Arguments X,Y to AXPY! are of different dimensions.")))) -;; note: we should optimize the calls to axpy from other axpy's since -;; we know exactly which one we are calling. ;; -;; also, the most common type of bug is with ! operators, e.g. when -;; you say (axpy! s a a) +(generate-typed-axpy!-func real-double-axpy!-typed double-float real-matrix-store-type real-matrix blas:daxpy) + +(defmethod axpy! ((alpha number) (x complex-matrix) (y real-matrix)) + (error "cannot AXPY! a complex X to a real Y, +don't know how to coerce COMPLEX to REAL")) -(deftype complex-double-float () - '(cl:complex (double-float * *))) +(defmethod axpy! ((alpha cl:real) (x real-matrix) (y real-matrix)) + (real-double-axpy!-typed (coerce alpha 'double-float) x y)) ;; +(generate-typed-axpy!-func complex-double-axpy!-typed complex-double-float complex-matrix-store-type complex-matrix blas:zaxpy) + +(defmethod axpy! ((alpha cl:real) (x real-matrix) (y complex-matrix)) + (real-double-axpy!-typed (coerce alpha 'double-float) x (realpart! y))) + +(defmethod axpy! ((alpha complex) (x real-matrix) (y complex-matrix)) + (real-double-axpy!-typed (coerce (realpart alpha) 'double-float) x (realpart! y)) + (real-double-axpy!-typed (coerce (imagpart alpha) 'double-float) x (imagpart! y))) + +(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix)) + (complex-double-axpy!-typed (complex-coerce alpha) x y)) + +;;;; (defgeneric axpy (alpha x y) (:documentation " @@ -118,187 +173,29 @@ ")) (defmethod axpy :before ((alpha number) (x standard-matrix) (y standard-matrix)) - (let ((nxm-x (number-of-elements x)) - (nxm-y (number-of-elements y))) - (declare (type fixnum nxm-x nxm-y)) - - (if (not (= nxm-x nxm-y)) - (error "arguments X and Y to AXPY not of the same size")))) + (mlet* (((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 (and (= nr-x nr-y) (= nc-x nc-y)) + (error "Arguments X,Y to AXPY are of different dimensions.")))) ;; (defmethod axpy ((alpha cl:real) (x real-matrix) (y real-matrix)) - (axpy (coerce alpha 'real-matrix-element-type) x y)) - -(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix) (y real-matrix)) - (let* ((nxm (number-of-elements y)) - (result (copy y))) - (declare (type fixnum nxm)) + (let ((result (copy y))) + (axpy! alpha x result))) - (daxpy nxm alpha (store x) 1 (store result) 1) - result)) - -;; -(defmethod axpy ((alpha cl:real) (x complex-matrix) (y real-matrix)) - (axpy (coerce alpha 'real-matrix-element-type) x y)) +(defmethod axpy ((alpha complex) (x real-matrix) (y real-matrix)) + (let ((result (scal alpha x))) + (axpy! 1d0 y result))) -(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix) (y real-matrix)) - (let* ((nxm (number-of-elements y)) - (n (nrows y)) - (m (ncols y)) - (result (make-complex-matrix-dim n m)) - (store-x (store x)) - (store-y (store y)) - (store-result (store result))) - (declare (type fixnum n m nxm) - (type (real-matrix-store-type *) store-y) - (type (complex-matrix-store-type *) store-x store-result)) - - (zcopy nxm store-x 1 store-result 1) ;; same as (COPY! x result) - (zdscal nxm alpha store-result 1) ;; same as (SCAL! alpha result) - (daxpy nxm 1.0d0 store-y 1 store-result 2) ;; same as (AXPY! 1d0 y result) - result)) - -;; -(defmethod axpy ((alpha cl:real) (x real-matrix) (y complex-matrix)) - (axpy (coerce alpha 'complex-matrix-element-type) x y)) - -(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix) (y complex-matrix)) - (let* ((nxm (number-of-elements y)) - (result (copy y))) - (declare (type fixnum nxm)) - (daxpy nxm alpha (store x) 1 (store result) 2) - result)) - -;; -(defmethod axpy ((alpha cl:real) (x complex-matrix) (y complex-matrix)) - (axpy (coerce alpha 'complex-matrix-element-type) x y)) - -(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix) (y complex-matrix)) - (let ((nxm (number-of-elements y)) - (result (copy y))) - (declare (type fixnum nxm)) - (daxpy (* 2 nxm) alpha (store x) 1 (store result) 1) - result)) - -;; -(defmethod axpy ((alpha number) (x real-matrix) (y complex-matrix)) - (let* ((nxm (number-of-elements y)) - (n (nrows y)) - (m (ncols y)) - (result (make-complex-matrix-dim n m)) - (store-x (store x)) - (store-y (store y)) - (store-result (store result)) - (c-alpha (complex-coerce alpha))) - (declare (type complex-double-float c-alpha) - (type fixnum n m nxm) - (type (real-matrix-store-type *) store-x) - (type (complex-matrix-store-type *) store-y store-result)) - - (dcopy nxm store-x 1 store-result 2) - (zscal nxm c-alpha store-result 1) - (daxpy (* 2 nxm) 1.0d0 store-y 1 store-result 1) - - result)) - -;; (defmethod axpy ((alpha number) (x complex-matrix) (y real-matrix)) - (let* ((nxm (number-of-elements y)) - (result (copy x)) - (store-result (store result)) - (c-alpha (complex-coerce alpha))) - (declare (type complex-double-float c-alpha) - (type fixnum nxm) - (type (complex-matrix-store-type *) store-result)) - - (zscal nxm c-alpha store-result 1) - (daxpy nxm 1.0d0 (store y) 1 store-result 2) - - result)) - -;; -(defmethod axpy ((alpha number) (x complex-matrix) (y complex-matrix)) - (let ((nxm (number-of-elements y)) - (result (copy y)) - (c-alpha (complex-coerce alpha))) - (declare (type complex-double-float c-alpha) - (type fixnum nxm)) - - (zaxpy nxm c-alpha (store x) 1 (store result) 1) - - result)) + (let ((result (scal alpha x))) + (axpy! 1d0 y result))) ;; -(defgeneric axpy! (alpha x y) - (:documentation - " - Syntax - ====== - (AXPY! alpha x y) - - Purpose - ======= - Same as AXPY except that the result - is stored in Y and Y is returned. -")) - -(defmethod axpy! :before ((alpha number) (x standard-matrix) (y standard-matrix)) - (let ((nxm-x (number-of-elements x)) - (nxm-y (number-of-elements y))) - (declare (type fixnum nxm-x nxm-y)) - - (if (not (= nxm-x nxm-y)) - (error "arguments X and Y to AXPY! not of the same size")))) - -(defmethod axpy! ((alpha number) (x complex-matrix) (y real-matrix)) - (error "cannot AXPY! a complex X to a real Y, -don't know how to coerce COMPLEX to REAL")) - -;; -(defmethod axpy! ((alpha cl:real) (x real-matrix) (y real-matrix)) - (axpy! (coerce alpha 'real-matrix-element-type) x y)) - -(defmethod axpy! ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix) (y real-matrix)) - (let* ((nxm (number-of-elements y))) - (declare (type fixnum nxm)) - - (daxpy nxm alpha (store x) 1 (store y) 1) - y)) - -;; -(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix)) - (let ((nxm (number-of-elements y)) - (c-alpha (complex-coerce alpha))) - (declare (type complex-double-float c-alpha) - (type fixnum nxm)) - - (daxpy (* 2 nxm) c-alpha (store x) 1 (store y) 1) - y)) - -;; -(defmethod axpy! ((alpha number) (x real-matrix) (y complex-matrix)) - (let* ((nxm (number-of-elements y)) - (store-x (store x)) - (store-y (store y)) - (c-alpha (complex-coerce alpha))) - (declare (type complex-double-float c-alpha) - (type fixnum nxm) - (type (real-matrix-store-type *) store-x) - (type (complex-matrix-store-type *) store-y)) - - (daxpy nxm (realpart c-alpha) store-x 1 store-y 2) - (with-vector-data-addresses ((addr-y store-y) - (addr-x store-x)) - (incf-sap :double-float addr-y) - (daxpy nxm (imagpart c-alpha) addr-x 1 addr-y 2)) - y)) +(defmethod axpy ((alpha number) (x real-matrix) (y complex-matrix)) + (let ((result (copy y))) + (axpy! alpha x result))) -;; -(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix)) - (let ((nxm (number-of-elements y)) - (c-alpha (complex-coerce alpha))) - (declare (type complex-double-float c-alpha) - (type fixnum nxm)) - - (zaxpy nxm c-alpha (store x) 1 (store y) 1) - y)) \ No newline at end of file +(defmethod axpy ((alpha number) (x complex-matrix) (y complex-matrix)) + (let ((result (copy y))) + (axpy! alpha x result))) \ No newline at end of file diff --git a/src/blas.lisp b/src/blas.lisp index 2875a5f..901b2fc 100644 --- a/src/blas.lisp +++ b/src/blas.lisp @@ -110,9 +110,9 @@ " (n :integer :input) (da :double-float :input) - (dx (* :double-float)) + (dx (* :double-float :inc head-x)) (incx :integer :input) - (dy (* :double-float) :output) + (dy (* :double-float :inc head-y) :output) (incy :integer :input) ) @@ -286,11 +286,11 @@ " (n :integer :input) (za :complex-double-float) - (zx (* :complex-double-float)) + (zx (* :complex-double-float :inc head-x)) (incx :integer :input) - (zy (* :complex-double-float) :output) + (zy (* :complex-double-float :inc head-y) :output) (incy :integer :input) - ) +) (def-fortran-routine zcopy :void " @@ -632,9 +632,9 @@ Y(0),Y(2*INCY), ... , Y(2*(N-1)*INCY) " (n :integer :input) - (dx (* :double-float) :input) + (dx (* :double-float :inc head-x) :input) (incx :integer :input) - (dy (* :double-float) :input) + (dy (* :double-float :inc head-y) :input) (incy :integer :input) ) @@ -744,12 +744,12 @@ (m :integer ) (n :integer ) (alpha :double-float ) - (a (* :double-float) ) + (a (* :double-float :inc head-a) ) (lda :integer ) - (x (* :double-float) ) + (x (* :double-float :inc head-x) ) (incx :integer ) (beta :double-float ) - (y (* :double-float) :output) + (y (* :double-float :inc head-y) :output) (incy :integer ) ) diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp index 64382d7..0e02c9b 100644 --- a/src/complex-matrix.lisp +++ b/src/complex-matrix.lisp @@ -10,6 +10,9 @@ (deftype complex-matrix-store-type (size) "The type of the storage structure for a COMPLEX-MATRIX" `(simple-array double-float (,size))) + +(deftype complex-double-float () + '(cl:complex (double-float * *))) ) ;; @@ -298,7 +301,8 @@ (defun imagpart! (mat) (cond + ((typep mat 'real-matrix) nil) ((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat) :nrows (nrows mat) :ncols (ncols mat) :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat)) - :head (+ 1 (* 2 (head mat))))))) \ No newline at end of file + :head (+ 1 (* 2 (head mat))))))) diff --git a/src/copy.lisp b/src/copy.lisp index 68dff10..cf987d9 100644 --- a/src/copy.lisp +++ b/src/copy.lisp @@ -204,7 +204,7 @@ don't know how to coerce a COMPLEX to a REAL")) (generate-typed-copy!-func complex-double-copy!-typed complex-matrix-store-type complex-matrix blas:zcopy) (generate-typed-num-copy!-func complex-double-num-copy!-typed - (complex (double-float * *)) complex-matrix-store-type complex-matrix + complex-double-float complex-matrix-store-type complex-matrix blas:zcopy (num (1x1-z-array diff --git a/src/gemm.lisp b/src/gemm.lisp index c9e8736..c42b892 100644 --- a/src/gemm.lisp +++ b/src/gemm.lisp @@ -202,8 +202,8 @@ (defmethod gemm! ((alpha cl:real) (a real-matrix) (b real-matrix) (beta cl:real) (c real-matrix) &optional (job :nn)) - (real-double-gemm!-typed (coerce alpha 'real-matrix-element-type) a b - (coerce beta 'real-matrix-element-type) c + (real-double-gemm!-typed (coerce alpha 'double-float) a b + (coerce beta 'double-float) c job)) ;; @@ -215,6 +215,26 @@ (complex-double-gemm!-typed (complex-coerce alpha) a b (complex-coerce beta) c job)) + +(defmethod gemm! ((alpha cl:real) (a real-matrix) (b complex-matrix) + (beta number) (c complex-matrix) + &optional (job :nn)) + (scal! beta c) + (real-double-gemm!-typed (coerce alpha 'double-float) a (realpart! b) + 1d0 (realpart! c) job) + (real-double-gemm!-typed (coerce alpha 'double-float) a (imagpart! b) + 1d0 (imagpart! c) job)) + +(defmethod gemm! ((alpha complex) (a real-matrix) (b complex-matrix) + (beta number) (c complex-matrix) + &optional (job :nn)) + (scal! beta c) + (real-double-gemm!-typed (coerce alpha 'double-float) a (realpart! b) + 1d0 (realpart! c) job) + (real-double-gemm!-typed (coerce alpha 'double-float) a (imagpart! b) + 1d0 (imagpart! c) job)) + + ;; (defmethod gemm! ((alpha number) (a standard-matrix) (b standard-matrix) (beta number) (c complex-matrix) @@ -306,4 +326,4 @@ (gemm! (complex-coerce alpha) a b (complex-coerce beta) c - job))) + job))) \ No newline at end of file diff --git a/src/gemv.lisp b/src/gemv.lisp new file mode 100644 index 0000000..9adb6cc --- /dev/null +++ b/src/gemv.lisp @@ -0,0 +1,35 @@ +(in-package :matlisp) + +(defun real-double-gemv!-typed (alpha A x beta y job) + (declare ((type double-float alpha beta) + (type real-matrix A x y) + (type symbol job))) + (symbol-macrolet + ((common-block + (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 ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum)) + ((st-x hd-x rs-x) (slot-values x '(store head row-stride)) + :type ((real-matrix-store-type *) fixnum fixnum)) + ((st-y hd-y rs-y) (slot-values y '(store head row-stride)) + :type ((real-matrix-store-type *) fixnum fixnum)) + ((sym lda tf-job) (blas-matrix-compatible-p A) :type (nil fixnum (string 1)))) + (if (not (string= tf-job "?")) + (blas:dgemv tf-job 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) + (dotimes (i nr-a) + (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot 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)))) + + (mlet* ((fortran-job (ecase job (:n "N") (:t "T")) :type ((string 1))) + ((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 ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum)) + ((st-x hd-x rs-x) (slot-values x '(store head row-stride)) + :type ((real-matrix-store-type *) fixnum fixnum)) + ((st-y hd-y rs-y) (slot-values y '(store head row-stride)) + :type ((real-matrix-store-type *) fixnum fixnum)) + ((sym lda tf-job) (blas-matrix-compatible-p A fortran-job) :type (nil fixnum (string 1)))) + (if (not (string= tf-job "?")) + (blas:dgemv tf-job 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) + (dotimes (i nr-a) + (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot 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)))) + \ No newline at end of file diff --git a/src/scal.lisp b/src/scal.lisp index 8548a26..88da573 100644 --- a/src/scal.lisp +++ b/src/scal.lisp @@ -123,7 +123,7 @@ how to coerce COMPLEX to REAL")) ;; (generate-typed-scal!-func complex-double-dscal!-typed double-float complex-matrix-store-type complex-matrix blas:zdscal) -(generate-typed-scal!-func complex-double-zscal!-typed (complex (double-float * *)) complex-matrix-store-type complex-matrix blas:zscal) +(generate-typed-scal!-func complex-double-zscal!-typed complex-double-float complex-matrix-store-type complex-matrix blas:zscal) (defmethod scal! ((alpha cl:real) (x complex-matrix)) (complex-double-dscal!-typed (coerce alpha 'double-float) x)) diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp index 47c4ff9..3c2cfce 100644 --- a/src/standard-matrix.lisp +++ b/src/standard-matrix.lisp @@ -277,21 +277,6 @@ matrix and a number")) (format stream "~%"))) ;; -(defun get-order-stride (matrix &optional (fortran-op "N")) - (check-type matrix standard-matrix) - (let ((rs (row-stride matrix)) - (cs (col-stride matrix))) - (declare (type fixnum rs cs)) - (cond - ((= cs 1) (values :row-major rs (cond - ((string= fortran-op "N" ) "T") - ((string= fortran-op "T" ) "N")))) - ((= rs 1) (values :col-major cs (cond - ((string= fortran-op "N" ) "N") - ((string= fortran-op "T" ) "T")))) - (t nil)))) - -;; (defun make-sub-matrix (matrix i j nrows ncols) (declare (type standard-matrix matrix) (type fixnum i j nrows ncols)) @@ -326,8 +311,20 @@ matrix and a number")) (cs (col-stride matrix) :type fixnum) (ne (number-of-elements matrix) :type fixnum)) (cond - ((= nc 1) (values t rs ne)) - ((= nr 1) (values t cs ne)) - ((= rs (* nc cs)) (values t cs ne)) - ((= cs (* nr rs)) (values t rs ne)) - (t (values nil -1 -1))))) \ No newline at end of file + ((or (= nc 1) (= cs (* nr rs))) (values t rs ne)) + ((or (= nr 1) (= rs (* nc cs))) (values t cs ne)) + (t (values nil -1 -1))))) + +;; +(defun blas-matrix-compatible-p (matrix &optional (fortran-op "N")) + (declare (optimize (safety 0) (speed 3)) + (type (or real-matrix complex-matrix) mat)) + (mlet* (((rs cs) (slot-values mat '(row-stride col-stride)) + :type (fixnum fixnum))) + (cond + ((= cs 1) (values :row-major rs (cond + ((string= fortran-op "N" ) "T") + ((string= fortran-op "T" ) "N")))) + ((= rs 1) (values :col-major cs fortran-op)) + ;;Lets not confound lisp's type declaration. + (t (values nil -1 "?"))))) \ No newline at end of file commit f4b94df81eab2264c071a5a6f592e6becf76f770 Merge: 2cbbcb6 8a2722a Author: Akshay Srinivasan <aks...@gm...> Date: Tue Mar 13 17:59:35 2012 +0530 Merge branch 'master' into matlisp-cffi Conflicts: src/blas.lisp diff --cc src/blas.lisp index 964f7b0,18c22d1..2875a5f --- a/src/blas.lisp +++ b/src/blas.lisp @@@ -494,7 -495,24 +495,24 @@@ (incy :integer :input) ) + (def-fortran-routine mzdotu :void + (result (* :complex-double-float) :output) + (n :integer :input) + (zx (* :complex-double-float) :input) + (incx :integer :input) + (zy (* :complex-double-float) :input) + (incy :integer :input) + ) -(defun zdotu (n zx incx zy incy) - (let ((result (make-array 2 :element-type 'double-float))) ++(let ((result (make-array 2 :element-type 'double-float))) ++ (defun zdotu (n zx incx zy incy) + (with-vector-data-addresses ((addr-result result) + (addr-zx zx) + (addr-zy zy)) + (mzdotu addr-result n addr-zx incx addr-zy incy) + (complex (aref result 0) (aref result 1))))) + + #-(and) (def-fortran-routine zdotc :complex-double-float " Syntax @@@ -538,6 -555,523 +556,24 @@@ (incy :integer :input) ) + (def-fortran-routine mzdotc :void + (result (* :complex-double-float) :output) + (n :integer :input) + (zx (* :complex-double-float) :input) + (incx :integer :input) + (zy (* :complex-double-float) :input) + (incy :integer :input) + ) + -(defun zdotc (n zx incx zy incy) - (let ((result (make-array 2 :element-type 'double-float))) ++(let ((result (make-array 2 :element-type 'double-float))) ++ (defun zdotc (n zx incx zy incy) + (with-vector-data-addresses ((addr-result result) + (addr-zx zx) + (addr-zy zy)) + (mzdotc addr-result n addr-zx incx addr-zy incy) + (complex (aref result 0) (aref result 1))))) + + -#| -(declaim (inline fortran-daxpy)) -(def-alien-routine ("daxpy_" fortran-daxpy) void - (n int :in-out) - (da double-float :in-out) - (dx (* double-float)) - (incx int :in-out) - (dy (* double-float)) - (incy int :in-out) -) -(defun daxpy ( - n - da - dx - incx - dy - incy - ) - (declare - (type fixnum n) - (type double-float da) - (type (simple-array double-float (*)) dx) - (type fixnum incx) - (type (simple-array double-float (*)) dy) - (type fixnum incy) - ) - (with-vector-data-addresses ( - (addr-dx dx) - (addr-dy dy) - ) - (multiple-value-bind ( return-value - new-n - new-incx - new-incy - ) - (fortran-daxpy - n - da - addr-dx - incx - addr-dy - incy - ) - (declare (ignore return-value new-n new-incx new-incy)) - (values - dy - )))) - -(declaim (inline fortran-dcopy)) -(def-alien-routine ("dcopy_" fortran-dcopy) void - (n int :in-out) - (dx (* double-float)) - (incx int :in-out) - (dy (* double-float)) - (incy int :in-out) -) -(defun dcopy ( - n - dx - incx - dy - incy - ) - (declare - (type fixnum n) - (type (simple-array double-float (*)) dx) - (type fixnum incx) - (type (simple-array double-float (*)) dy) - (type fixnum incy) - ) - (with-vector-data-addresses ( - (addr-dx dx) - (addr-dy dy) - ) - (multiple-value-bind ( return-value - new-n - new-incx - new-incy - ) - (fortran-dcopy - n - addr-dx - incx - addr-dy - incy - ) - (declare (ignore return-value new-n new-incx new-incy)) - (values - dy - )))) - - -(declaim (inline fortran-drot)) -(def-alien-routine ("drot_" fortran-drot) void - (n int :in-out) - (dx (* double-float)) - (incx int :in-out) - (dy (* double-float)) - (incy int :in-out) - (c double-float :in-out) - (s double-float :in-out) -) -(defun drot ( - n - dx - incx - dy - incy - c - s - ) - (declare - (type fixnum n) - (type (simple-array double-float (*)) dx) - (type fixnum incx) - (type (simple-array double-float (*)) dy) - (type fixnum incy) - (type double-float c) - (type double-float s) - ) - (with-vector-data-addresses ( - (addr-dx dx) - (addr-dy dy) - ) - (multiple-value-bind ( return-value - new-n - new-incx - new-incy - new-c - new-s - ) - - (fortran-drot - n - addr-dx - incx - addr-dy - incy - c - s - ) - (declare (ignore return-value new-n new-incx new-incy)) - (values - dx dy new-c new-s )))) - -(declaim (inline fortran-drotg)) -(def-alien-routine ("drotg_" fortran-drotg) void - (da double-float :in-out) - (db double-float :in-out) - (c double-float :in-out) - (s double-float :in-out) -) -(defun drotg ( - da - db - c - s - ) - (declare - (type double-float da) - (type double-float db) - (type double-float c) - (type double-float s) - ) - (with-vector-data-addresses ( - ) - (multiple-value-bind ( return-value new-da new-db new-c new-s - ) - (fortran-drotg - da - db - c - s - ) - (declare (ignore return-value)) - (values - new-da new-db new-c new-s )))) - -(declaim (inline fortran-dscal)) -(def-alien-routine ("dscal_" fortran-dscal) void - (n int :in-out) - (da double-float :in-out) - (dx (* double-float)) - (incx int :in-out) -) -(defun dscal ( - n - da - dx - incx - ) - (declare - (type fixnum n) - (type double-float da) - (type (simple-array double-float (*)) dx) - (type fixnum incx) - ) - (with-vector-data-addresses ( - (addr-dx dx) - ) - (multiple-value-bind ( return-value - new-n - new-incx - ) - (fortran-dscal - n - da - addr-dx - incx - ) - (declare (ignore return-value new-n new-incx)) - (values - dx )))) - -(declaim (inline fortran-dswap)) -(def-alien-routine ("dswap_" fortran-dswap) void - (n int :in-out) - (dx (* double-float)) - (incx int :in-out) - (dy (* double-float)) - (incy int :in-out) -) -(defun dswap ( - n - dx - incx - dy - incy - ) - (declare - (type fixnum n) - (type (simple-array double-float (*)) dx) - (type fixnum incx) - (type (simple-array double-float (*)) dy) - (type fixnum incy) - ) - (with-vector-data-addresses ( - (addr-dx dx) - (addr-dy dy) - ) - (multiple-value-bind ( return-value - new-n - new-incx - new-incy - ) - (fortran-dswap - n - addr-dx - incx - addr-dy - incy - ) - (declare (ignore return-value new-n new-incx new-incy)) - (values - dx )))) - - -(declaim (inline fortran-zaxpy)) -(def-alien-routine ("zaxpy_" fortran-zaxpy) void - (n int :in-out) - (za (* double-float)) - (zx (* double-float)) - (incx int :in-out) - (zy (* double-float)) - (incy int :in-out) -) -(defun zaxpy ( - n - za - zx - incx - zy - incy - ) - (declare - (type fixnum n) - (type (simple-array double-float (*)) za) - (type (simple-array double-float (*)) zx) - (type fixnum incx) - (type (simple-array double-float (*)) zy) - (type fixnum incy) - ) - (with-vector-data-addresses ( - (addr-za za) - (addr-zx zx) - (addr-zy zy) - ) - (multiple-value-bind ( return-value - new-n - new-incx - new-incy - ) - (fortran-zaxpy - n - addr-za - addr-zx - incx - addr-zy - incy - ) - (declare (ignore return-value new-n new-incx new-incy)) - (values - zy )))) - -(declaim (inline fortran-zcopy)) -(def-alien-routine ("zcopy_" fortran-zcopy) void - (n int :in-out) - (zx (* double-float)) - (incx int :in-out) - (zy (* double-float)) - (incy int :in-out) -) -(defun zcopy ( - n - zx - incx - zy - incy - ) - (declare - (type fixnum n) - (type (simple-array double-float (*)) zx) - (type fixnum incx) - (type (simple-array double-float (*)) zy) - (type fixnum incy) - ) - (with-vector-data-addresses ( - (addr-zx zx) - (addr-zy zy) - ) - (multiple-value-bind ( return-value - new-n - new-incx - new-incy - ) - (fortran-zcopy - n - addr-zx - incx - addr-zy - incy - ) - (declare (ignore return-value new-n new-incx new-incy)) - (values - zy )))) - - -(declaim (inline fortran-zdscal)) -(def-alien-routine ("zdscal_" fortran-zdscal) void - (n int :in-out) - (da double-float :in-out) - (zx (* double-float)) - (incx int :in-out) -) -(defun zdscal ( - n - da - zx - incx - ) - (declare - (type fixnum n) - (type double-float da) - (type (simple-array double-float (*)) zx) - (type fixnum incx) - ) - (with-vector-data-addresses ( - (addr-zx zx) - ) - (multiple-value-bind ( return-value - new-n - new-incx - ) - (fortran-zdscal - n - da - addr-zx - incx - ) - (declare (ignore return-value new-n new-incx)) - (values - zx )))) - -(declaim (inline fortran-zrotg)) -(def-alien-routine ("zrotg_" fortran-zrotg) void - (ca (* double-float)) - (cb (* double-float)) - (c double-float :in-out) - (s (* double-float)) -) -(defun zrotg ( - ca - cb - c - s - ) - (declare - (type (simple-array double-float (*)) ca) - (type (simple-array double-float (*)) cb) - (type double-float c) - (type (simple-array double-float (*)) s) - ) - (with-vector-data-addresses ( - (addr-ca ca) - (addr-cb cb) - (addr-s s) - ) - (multiple-value-bind ( return-value new-c - ) - (fortran-zrotg - addr-ca - addr-cb - c - addr-s - ) - (declare (ignore return-value)) - (values - ca cb new-c s )))) - -(declaim (inline fortran-zscal)) -(def-alien-routine ("zscal_" fortran-zscal) void - (n int :in-out) - (za (* double-float)) - (zx (* double-float)) - (incx int :in-out) -) -(defun zscal ( - n - za - zx - incx - ) - (declare - (type fixnum n) - (type (simple-array double-float (*)) za) - (type (simple-array double-float (*)) zx) - (type fixnum incx) - ) - (with-vector-data-addresses ( - (addr-za za) - (addr-zx zx) - ) - (multiple-value-bind ( return-value - new-n - new-incx - ) - (fortran-zscal - n - addr-za - addr-zx - incx - ) - (declare (ignore return-value new-n new-incx)) - (values - zx )))) - -(declaim (inline fortran-zswap)) -(def-alien-routine ("zswap_" fortran-zswap) void - (n int :in-out) - (zx (* double-float)) - (incx int :in-out) - (zy (* double-float)) - (incy int :in-out) -) -(defun zswap ( - n - zx - incx - zy - incy - ) - (declare - (type fixnum n) - (type (simple-array double-float (*)) zx) - (type fixnum incx) - (type (simple-array double-float (*)) zy) - (type fixnum incy) - ) - (with-vector-data-addresses ( - (addr-zx zx) - (addr-zy zy) - ) - (multiple-value-bind ( return-value - new-n - new-incx - new-incy - ) - (fortran-zswap - n - addr-zx - incx - addr-zy - incy - ) - (declare (ignore return-value new-n new-incx new-incy)) - (values - zx )))) - -|# (def-fortran-routine idamax :integer " " ----------------------------------------------------------------------- Summary of changes: Makefile.am | 1 + configure | 130 ++++++---------------- configure.ac | 254 +++++++++++++++++++++++++++--------------- lib-src/compat/Makefile.am | 12 ++ lib-src/compat/compat.f | 34 ++++++ lib/lazy-loader.lisp.in | 61 +++++++--- src/axpy.lisp | 269 ++++++++++++++------------------------------ src/blas.lisp | 56 ++++++++-- src/complex-matrix.lisp | 6 +- src/copy.lisp | 57 +++++----- src/gemm.lisp | 26 ++++- src/gemv.lisp | 32 +++++ src/scal.lisp | 26 +++-- src/standard-matrix.lisp | 37 +++---- tests/blas.lisp | 8 +- 15 files changed, 541 insertions(+), 468 deletions(-) create mode 100644 lib-src/compat/Makefile.am create mode 100644 lib-src/compat/compat.f create mode 100644 src/gemv.lisp hooks/post-receive -- matlisp |