[Matlisp-commit] [matlisp-git]matlisp branch matlisp-cffi updated. 2012-02-23-55-g54341c2

[Akshay S. <ak...@us...>]

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  54341c25f149263190e4ffad1c516d93a79ad3ed (commit)
       via  146847f922138620cb4b1ad064cc2ad8f80bb304 (commit)
      from  517949e31b41b303dab91670e80b207bc45d3256 (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 54341c25f149263190e4ffad1c516d93a79ad3ed
Author: Akshay Srinivasan <aks...@gm...>
Date:   Sat Mar 17 11:32:32 2012 +0530

    -> Gemv works!

diff --git a/packages.lisp b/packages.lisp
index 566dbd6..74385fa 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -154,8 +154,10 @@
 ;;; Define the packages and symbols for Matlisp.
 
 (defpackage :utilities
-  (:use :common-lisp)
-  (:export #:lcase
+  (:use :common-lisp)  
+  (:export #:ensure-list
+	   #:zip
+	   #:zip-eq
 	   #:cut-cons-chain!
 	   #:when-let
 	   #:if-ret
diff --git a/src/gemv.lisp b/src/gemv.lisp
index 0e49fac..c58f746 100644
--- a/src/gemv.lisp
+++ b/src/gemv.lisp
@@ -39,8 +39,6 @@
 						      (* beta (matrix-ref-2d y i 0)))))))))
     y))
 
-
-
 ;;
 (defgeneric gemv! (alpha A x beta y &optional job)
   (:documentation
@@ -100,6 +98,7 @@
 
 (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))
 
@@ -110,57 +109,126 @@
 
 (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 cl:real) (y complex-matrix) &optional (job :n))
-  (real-double-gemv!-typed (coerce alpha 'double-float) A x
-			   (coerce beta 'double-float) (realpart! y) job))
-
-(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
+(defmethod gemv! ((alpha number) (A real-matrix) (x real-matrix)
 		  (beta complex) (y complex-matrix) &optional (job :n))
+  ;; y <- \beta * y  
   (scal! (complex-coerce beta) y)
-  (real-double-gemv!-typed (coerce alpha 'double-float) A x
-			   1d0 (realpart! y) job))
+  ;; 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 (realpart! y)))
+    (declare (type double-float r-be r-al)
+	     (type real-matrix r-y))
+    ;; y <- \beta * y
+    (scal! r-be y)
+    ;; (realpart! y) <- (realpart! y) + \alpha * A o x
+    (real-double-gemv!-typed r-al A x 1d0 r-y job))
+  y)
 
 (defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
 		  (beta cl:real) (y complex-matrix) &optional (job :n))
-  (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
-			   (coerce beta 'double-float) (realpart! y) job)
-  (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
-			   (coerce beta 'double-float) (imagpart! y) job))
+  (let ((r-al (coerce (realpart alpha) 'double-float))
+	(i-al (coerce (imagpart alpha) 'double-float))
+	(r-be (coerce beta 'double-float))
+	(r-y (realpart! y))
+	(i-y (imagpart! y)))
+    (declare (type double-float r-al i-al r-be)
+	     (type real-matrix r-y i-y))
+    ;; (realpart! y) <- \beta * (realpart! y) + (realpart \alpha) . A o x
+    (real-double-gemv!-typed r-al A x r-be r-y job)
+    ;; (imagpart! y) <- \beta * (imagpart! y) + (imagpart \alpha) . A o x
+    (real-double-gemv!-typed i-al A x r-be i-y job)) 
+  y)
 
-(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
-		  (beta complex) (y complex-matrix) &optional (job :n))
-  (scal! (complex-coerce beta) y)
-  (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
-			   1d0 (realpart! y) job)
-  (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
-			   1d0 (imagpart! y) job))
 ;
-
 (defmethod gemv! ((alpha number) (A real-matrix) (x complex-matrix)
-		  (beta number) (y complex-matrix) &optional (job :n))
-  (gemv! alpha A (realpart! x)
-	 beta y job)
-  (gemv! (* #c(0d0 1d0) alpha) A (imagpart! x)
-	 beta y job))
-;
-(defmethod gemv! ((alpha number) (A complex-matrix) (x real-matrix)
-		  (beta number) (y complex-matrix) &optional (job :n))
-  (gemv! alpha (realpart! A) x
-	 beta y job)
-  (gemv! (* #c(0d0 1d0) alpha) (imagpart! A) x
-	 beta y job))
+		  (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 (realpart! x))
+	(i-x (imagpart! x))
+	(r-y (realpart! y))
+	(i-y (imagpart! 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))
+    ;; (realpart! y) <- \beta * (realpart! y) + \alpha . A o (realpart! x)
+    (real-double-gemv!-typed r-al A r-x r-be r-y job)
+    ;; (imagpart! y) <- \beta * (imagpart! y) + \alpha . A o (realpart! 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 (realpart! x))
+	(i-x (imagpart! x))
+	(r-y (realpart! y))
+	(i-y (imagpart! 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)
 
- 
 ;
-(defun gemv! ((alpha number) (A complex-matrix) (x real-matrix)
-	      (beta number) (y complex-matrix))
-  )
+(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))
 
-;;
\ No newline at end of file
+(defmethod gemv! ((alpha cl:real) (A complex-matrix) (x real-matrix)
+		  (beta cl:real) (y complex-matrix) &optional (job :n))
+  (let ((r-A (realpart! A))
+	(i-A (imagpart! A))
+	(r-y (realpart! y))
+	(i-y (imagpart! 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))
+    ;; (realpart! y) <- \beta * (realpart! y) + \alpha . A o (realpart! x)
+    (real-double-gemv!-typed r-al r-A x r-be r-y job)
+    ;; (imagpart! y) <- \beta * (imagpart! y) + \alpha . A o (realpart! x)
+    (real-double-gemv!-typed r-al i-A x r-be i-y job))
+  y)
+
+(defmethod gemv! ((alpha complex) (A complex-matrix) (x real-matrix)
+		  (beta cl:real) (y complex-matrix) &optional (job :n))
+  (let ((r-A (realpart! A))
+	(i-A (imagpart! A))
+	(r-y (realpart! y))
+	(i-y (imagpart! 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))
+  y)
\ No newline at end of file
diff --git a/src/utilities.lisp b/src/utilities.lisp
index 76caf8c..6cbf2bc 100644
--- a/src/utilities.lisp
+++ b/src/utilities.lisp
@@ -145,15 +145,23 @@ else run else-body"
     (cut-cons-chain-tin lst test lst)))
 
 ;;
-(defmacro lcase (keyform &rest body)
-  (let ((key-eval (gensym)))
-    (labels ((app-equal (lst)
-	       (if (null lst)
-		   nil
-		   `(((equal ,key-eval ,(caar lst)) ,@(cdar lst))
-		     ,@(app-equal (cdr lst))))))
-      `(let ((,key-eval ,keyform))
-	 (cond
-	   ,@(app-equal body))))))
+(defun zip (&rest args)  
+  (apply #'map 'list #'list args))
+
+;;
+(defmacro mcase (keyform &rest body)
+  (labels ((app-equal (lst)
+	     (if (null lst)
+		 nil
+		 `(((and ,@(mapcar (lambda (pair) (cons 'eq pair))
+				   (zip keyform (caar lst))))
+		    ,@(cdar lst))
+		   ,@(app-equal (cdr lst))))))
+    `(cond
+       ,@(app-equal body))))
+
+(defmacro zip-eq (a b)
+  `(and ,@(mapcar (lambda (pair) (cons 'eq pair))
+		  (zip (ensure-list a) (ensure-list b)))))
 
 ;;

commit 146847f922138620cb4b1ad064cc2ad8f80bb304
Author: Akshay Srinivasan <aks...@gm...>
Date:   Sat Mar 17 08:53:49 2012 +0530

    Gemv core function works reasonably well.

diff --git a/src/axpy.lisp b/src/axpy.lisp
index d6cee37..f001fd8 100644
--- a/src/axpy.lisp
+++ b/src/axpy.lisp
@@ -85,20 +85,18 @@
 	     ((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)))
+		(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 (> nr-a nc-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 alpha st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
+     mat-b))
 
 ;;
 (defgeneric axpy! (alpha x y)
diff --git a/src/blas.lisp b/src/blas.lisp
index 901b2fc..e7f76d9 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -498,19 +498,16 @@
 (def-fortran-routine mzdotu :void
   (result (* :complex-double-float) :output)
   (n :integer :input)
-  (zx (* :complex-double-float) :input)
+  (zx (* :complex-double-float :inc head-x) :input)
   (incx :integer :input)
-  (zy (* :complex-double-float) :input)
+  (zy (* :complex-double-float :inc head-y) :input)
   (incy :integer :input)
   )
 
 (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)))))
+  (defun zdotu (n zx incx zy incy &key (head-x 0) (head-y 0))
+    (mzdotu result n zx incx zy incy :head-x head-x :head-y head-y)
+    (complex (aref result 0) (aref result 1))))
 
 #-(and)
 (def-fortran-routine zdotc :complex-double-float
@@ -559,19 +556,16 @@
 (def-fortran-routine mzdotc :void
   (result (* :complex-double-float) :output)
   (n :integer :input)
-  (zx (* :complex-double-float) :input)
+  (zx (* :complex-double-float :inc head-x) :input)
   (incx :integer :input)
-  (zy (* :complex-double-float) :input)
+  (zy (* :complex-double-float :inc head-y) :input)
   (incy :integer :input)
   )
 
 (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)))))
+  (defun zdotc (n zx incx zy incy &key (head-x 0) (head-y 0))
+    (mzdotc result n zx incx zy incy :head-x head-x :head-y head-y)
+    (complex (aref result 0) (aref result 1))))
 
 
 (def-fortran-routine idamax :integer
@@ -2143,12 +2137,12 @@
   (m :integer )
   (n :integer )
   (alpha :complex-double-float )
-  (a (* :complex-double-float) )
+  (a (* :complex-double-float :inc head-a) )
   (lda :integer )
-  (x (* :complex-double-float) )
+  (x (* :complex-double-float :inc head-x) )
   (incx :integer )
   (beta :complex-double-float )
-  (y (* :complex-double-float) :output)
+  (y (* :complex-double-float :inc head-y) :output)
   (incy :integer )
 )
 
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 0e02c9b..704f823 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -97,8 +97,8 @@
 	 (store (allocate-complex-store size)))
     (multiple-value-bind (row-stride col-stride)
 	(ecase order
-	  (:row-major (values n 1))
-	  (:col-major (values 1 m)))
+	  (:row-major (values m 1))
+	  (:col-major (values 1 n)))
       (let ((matrix
 	     (make-instance 'complex-matrix
 			    :nrows n :ncols m
@@ -306,3 +306,17 @@
 						:nrows (nrows mat) :ncols (ncols mat)
 						:row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
 						:head (+ 1 (* 2 (head mat)))))))
+
+
+(defun conjugate! (mat)
+  (cond
+    ((typep mat 'real-matrix) nil)
+    ((typep mat 'complex-matrix) (progn
+				   (transpose! (scal! -1d0 (imagpart! mat)))
+				   mat))))
+
+(defmacro with-conjugate! (matlst &rest body)
+  `(progn
+     ,@(mapcar #'(lambda (mat) `(conjugate! ,mat)) matlst)
+     ,@body
+     ,@(mapcar #'(lambda (mat) `(conjugate! ,mat)) matlst)))
\ No newline at end of file
diff --git a/src/copy.lisp b/src/copy.lisp
index 03c6e8e..8ce687d 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -82,7 +82,6 @@
 (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)
diff --git a/src/gemv.lisp b/src/gemv.lisp
index cffa7d6..0e49fac 100644
--- a/src/gemv.lisp
+++ b/src/gemv.lisp
@@ -1,32 +1,166 @@
 (in-package :matlisp)
 
-(defun real-double-gemv!-typed (alpha A x beta y job)
-  ;;Be very careful when you use this function!
+;;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
-  ;;(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-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
+  `(defun ,func (alpha A x beta y job)
+    (declare (type ,element-type alpha beta)
+	     (type ,matrix-type 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 ((,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 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-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 (string= fort-op "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))
+
+
+
+;;
+(defgeneric gemv! (alpha A x beta y &optional job)
+  (:documentation
+"
+  Syntax
+  ======
+  (GEMV! alpha A x beta y [job])
+
+  Purpose
+  =======
+  Performs the GEneral Matrix Vector operation given by
+               --      -      -
+
+            Y <- alpha * op(A) * x + beta * y
+     
+  and returns y.
+
+  alpha,beta are scalars,
+  A is a matrix, and x,y are vectors.
+
+  op(A) means either A or A'.
+
+     JOB                    Operation
+  ---------------------------------------------------
+     :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))
+  (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))
+  (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 cl:real) (y complex-matrix) &optional (job :n))
+  (real-double-gemv!-typed (coerce alpha 'double-float) A x
+			   (coerce beta 'double-float) (realpart! y) job))
+
+(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
+		  (beta complex) (y complex-matrix) &optional (job :n))
+  (scal! (complex-coerce beta) y)
+  (real-double-gemv!-typed (coerce alpha 'double-float) A x
+			   1d0 (realpart! y) job))
+
+;
+
+(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
+		  (beta cl:real) (y complex-matrix) &optional (job :n))
+  (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
+			   (coerce beta 'double-float) (realpart! y) job)
+  (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
+			   (coerce beta 'double-float) (imagpart! y) job))
+
+(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
+		  (beta complex) (y complex-matrix) &optional (job :n))
+  (scal! (complex-coerce beta) y)
+  (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
+			   1d0 (realpart! y) job)
+  (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
+			   1d0 (imagpart! y) job))
+;
+
+(defmethod gemv! ((alpha number) (A real-matrix) (x complex-matrix)
+		  (beta number) (y complex-matrix) &optional (job :n))
+  (gemv! alpha A (realpart! x)
+	 beta y job)
+  (gemv! (* #c(0d0 1d0) alpha) A (imagpart! x)
+	 beta y job))
+;
+(defmethod gemv! ((alpha number) (A complex-matrix) (x real-matrix)
+		  (beta number) (y complex-matrix) &optional (job :n))
+  (gemv! alpha (realpart! A) x
+	 beta y job)
+  (gemv! (* #c(0d0 1d0) alpha) (imagpart! A) x
+	 beta y job))
+
+ 
+;
+(defun gemv! ((alpha number) (A complex-matrix) (x real-matrix)
+	      (beta number) (y complex-matrix))
+  )
+
+;;
\ No newline at end of file
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index 0f1ec71..089ef9b 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -228,25 +228,6 @@ Matrix only has ~A elements." idx (number-of-elements matrix))))
   (list (nrows matrix) (ncols matrix)))
 
 ;;
-(defgeneric transpose! (matrix)
-  (:documentation
-   "
-   Syntax
-   ======
-   (transpose! matrix)
-
-   Purpose
-   =======
-   Exchange row and column strides so that effectively
-   the matrix is transposed in place (without much effort).
-"))
-
-(defmethod transpose! ((matrix standard-matrix))
-  (rotatef (nrows matrix) (ncols matrix))
-  (rotatef (row-stride matrix) (col-stride matrix))
-  matrix)
-
-;;
 (defgeneric fill-matrix (matrix fill-element)
   (:documentation 
    "
@@ -277,7 +258,58 @@ matrix and a number"))
     (format stream "~%")))
 
 ;;
-(defun make-sub-matrix (matrix i j nrows ncols)
+(defun transpose-i! (matrix)
+"
+   Syntax
+   ======
+   (transpose-i! matrix)
+
+   Purpose
+   =======
+   Exchange row and column strides so that effectively
+   the matrix is transposed in place (without much effort).
+"
+  (cond
+    ((typep matrix 'standard-matrix)
+     (progn
+       (rotatef (nrows matrix) (ncols matrix))
+       (rotatef (row-stride matrix) (col-stride matrix))
+       matrix))
+    ((typep matrix 'number) matrix)
+    (t (error "Don't know how to take the transpose of ~A." matrix))))
+
+;;
+(defun transpose! (matrix)
+"
+   Syntax
+   ======
+   (transpose! matrix)
+
+   Purpose
+   =======
+   Create a new matrix object which represents the transpose of the
+   the given matrix, but shares the store with matrix.
+"
+  (cond
+    ((typep matrix 'standard-matrix)
+     (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
+	     (make-instance (class-of matrix)
+			    :nrows nc :ncols nr
+			    :store (store matrix)
+			    :head hd
+			    :row-stride cs :col-stride rs)))
+    ((typep matrix 'number) matrix)
+    (t (error "Don't know how to take the transpose of ~A." matrix))))
+
+;;
+(defmacro with-transpose! (matlst &rest body)
+  `(progn
+     ,@(mapcar #'(lambda (mat) `(transpose-i! ,mat)) matlst)
+     ,@body
+     ,@(mapcar #'(lambda (mat) `(transpose-i! ,mat)) matlst)))
+
+;;
+(defun sub! (matrix i j nrows ncols)
   (declare (type standard-matrix matrix)
 	   (type fixnum i j nrows ncols))
   (let ((hd (head matrix))
@@ -286,20 +318,49 @@ matrix and a number"))
 	(rs (row-stride matrix))
 	(cs (col-stride matrix)))
     (declare (type fixnum hd nr nc rs cs))
-    (when (or (> i (- nr nrows)) (> j (- nc ncols)))
-      (error "Parent matrix is not big enough to create sub-matrix."))
+    (unless (and (< -1 i (+ i nrows) nr) (< -1 j (+ j ncols) nc))
+      (error "Bad index and/or size.
+Cannot create a sub-matrix of size (~a ~a) starting at (~a ~a)" nrows ncols i j))
     (make-instance (class-of matrix)
 		   :nrows nrows :ncols ncols
 		   :store (store matrix)
 		   :head (store-indexing i j hd rs cs)
 		   :row-stride rs :col-stride cs)))
 
+(defun (setf sub!) (mat-b mat-a i j nrows ncols)
+  (copy! mat-b (sub! mat-a i j nrows ncols)))
+
 ;;
-(defmacro with-transpose! (matlst &rest body)
-  `(progn
-     ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
-     ,@body
-     ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
+(defun row! (matrix i)
+  (declare (type standard-matrix matrix)
+	   (type fixnum i))
+  (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
+	 (unless (< -1 i nr)
+	   (error "Index ~a is outside the valid range for the given matrix." i))
+	 (make-instance (class-of matrix)
+			:nrows 1 :ncols nc
+			:store (store matrix)
+			:head (store-indexing i 0 hd rs cs)
+			:row-stride rs :col-stride cs)))
+
+(defun (setf row!) (mat-b mat-a i)
+  (copy! mat-b (row! mat-a i)))
+
+;;
+(defun col! (matrix j)
+  (declare (type standard-matrix matrix)
+	   (type fixnum j))
+  (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
+	 (unless (< -1 j nc)
+	   (error "Index ~a is outside the valid range for the given matrix." j))
+	 (make-instance (class-of matrix)
+			:nrows nr :ncols 1
+			:store (store matrix)
+			:head (store-indexing 0 j hd rs cs)
+			:row-stride rs :col-stride cs)))
+
+(defun (setf col!) (mat-b mat-a j)
+  (copy! mat-b (col! mat-a j)))
 
 ;;
 (defun blas-copyable-p (matrix)
@@ -325,6 +386,8 @@ matrix and a number"))
 	   ((= 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))
+	   ((= rs 1) (values :col-major cs (cond
+					     ((string= fortran-op "N" ) "N")
+					     ((string= fortran-op "N" ) "T"))))
 	   ;;Lets not confound lisp's type declaration.
 	   (t (values nil -1 "?")))))
\ No newline at end of file

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

Summary of changes:
 packages.lisp            |    6 +-
 src/axpy.lisp            |   26 ++---
 src/blas.lisp            |   32 +++----
 src/complex-matrix.lisp  |   18 +++-
 src/copy.lisp            |    1 -
 src/gemv.lisp            |  258 +++++++++++++++++++++++++++++++++++++++++-----
 src/standard-matrix.lisp |  119 ++++++++++++++++-----
 src/utilities.lisp       |   28 +++--
 8 files changed, 384 insertions(+), 104 deletions(-)


hooks/post-receive
-- 
matlisp