diff --git a/gsll.asd b/gsll.asd
index 3c2b90320303cd7fe638629154352f64d9e21962..1fc4d0320a38cfc196820a6dd2e68a01aed7705b 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2008-08-06 22:41:57EDT gsll.asd>
+;; Time-stamp: <2008-08-10 17:01:08EDT gsll.asd>
 ;; $Id$
 
 (asdf:defsystem "gsll"
@@ -81,7 +81,7 @@
 	    :components
 	    ((:file "blas1")
 	     (:file "blas2")
-	     ;(:file "blas3" :depends-on (blas2))
+	     (:file "blas3" :depends-on (blas2))
 	     ;(:file "lu")
 	     ;(:file "qr")
 	     ;(:file "qrpt")
diff --git a/linear-algebra/blas2.lisp b/linear-algebra/blas2.lisp
index cf604f48a1943693306d3368635d85a63b9cfa4b..642ac4ed85887c393116bf2eeb27bd8872ddf58c 100644
--- a/linear-algebra/blas2.lisp
+++ b/linear-algebra/blas2.lisp
@@ -1,6 +1,6 @@
 ;; BLAS level 2, Matrix-vector operations
 ;; Liam Healy, Wed Apr 26 2006 - 21:08
-;; Time-stamp: <2008-08-06 22:49:49EDT blas2.lisp>
+;; Time-stamp: <2008-08-09 19:43:20EDT blas2.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -25,8 +25,9 @@
 ;;;; Functions
 ;;;;****************************************************************************
 
-(defmfun matrix-vector-product
-    (TransA alpha (A matrix) (x vector) beta (y vector))
+(defmfun matrix-product
+    ((A matrix) (x vector) (y vector)
+     &optional (alpha 1) (beta 1) (TransA :notrans) TransB)
   ("gsl_blas_" :type "gemv")
   ((transa cblas-transpose) (alpha :element-c-type) ((mpointer A) :pointer)
    ((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer))
@@ -34,32 +35,55 @@
   :element-types :float-complex
   :inputs (A x y)
   :outputs (y)
+  :return (y)
   :documentation			; FDL
-  "The matrix-vector product and sum
+  "If the second and third arguments are vectors, compute
+   the matrix-vector product and sum
     y = alpha op(A) x + beta y, where op(A) = A, A^T, A^H
-    for TransA = :notrans, :trans, :conjtrans.")
+  for TransA = :notrans, :trans, :conjtrans.
+  If the second and third arguments are matrices, compute
+  the matrix-matrix product and sum C = alpha
+  op(A) op(B) + beta C where op(A) = A, A^T, A^H for TransA =
+  :notrans, :trans, :conjtrans and similarly for the
+  parameter TransB.")
 
-(defmfun matrix-vector-product-triangular
-    (uplo TransA diag (A matrix) (x vector))
+(defmfun matrix-product-triangular
+    ((A matrix) (x vector)
+     &optional (alpha 1) (uplo :upper) (TransA :notrans)
+     (diag :nonunit) (side :left))
   ("gsl_blas_" :type "trmv")
   ((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag)
    ((mpointer A) :pointer) ((mpointer x) :pointer))
   :definition :generic
   :element-types :float-complex
-  :inputs (TransA A x)
+  :inputs (A x)
   :outputs (x)
+  :return (x)
   :documentation			; FDL
-  "The matrix-vector product x =\alpha op(A) x
+  "If the second argument is a vector, compute
+   the matrix-vector product x =alpha op(A) x
    for the triangular matrix A, where op(A) = A, A^T, A^H for
    TransA = :NoTrans, :Trans, :ConjTrans. When Uplo
    is :Upper then the upper triangle of A is used, and when
    Uplo is :Lower then the lower triangle of A is used. If
    Diag is :NonUnit then the diagonal of the matrix is used,
    but if Diag is :Unit then the diagonal elements of the
-   matrix A are taken as unity and are not referenced.")
+   matrix A are taken as unity and are not referenced.
+   If the second argument is a matrix, compute
+   the matrix-matrix product B = alpha op(A) B
+   if Side is :Left and B = \alpha B op(A) if Side is
+   :Right. The matrix A is triangular and op(A) = A, A^T, A^H
+   for TransA = :NoTrans, :Trans, :ConjTrans When Uplo
+   is :Upper then the upper triangle of A is used, and when
+   Uplo is :Lower then the lower triangle of A is used. If
+   Diag is :NonUnit then the diagonal of A is used, but if
+   Diag is :Unit then the diagonal elements of the matrix A
+   are taken as unity and are not referenced.")
 
-(defmfun inverse-matrix-vector-product
-    (uplo TransA diag (A matrix) (x vector))
+(defmfun inverse-matrix-product
+    ((A matrix) (x vector)
+     &optional (alpha 1) (uplo :upper) (TransA :notrans)
+     (diag :nonunit) (side :left))
   ("gsl_blas_" :type "trsv")
   ((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag)
    ((mpointer A) :pointer) ((mpointer x) :pointer))
@@ -67,17 +91,30 @@
   :element-types :float-complex
   :inputs (A x)
   :outputs (x)
+  :return (x)
   :documentation			; FDL
-  "Compute inv(op(A)) x for x, where op(A) = A, A^T, A^H for
-    TransA = :NoTrans, :Trans, :ConjTrans. When Uplo is
-    :Upper then the upper triangle of A is used, and when Uplo is
-    :Lower then the lower triangle of A is used. If Diag is
-    :NonUnit then the diagonal of the matrix is used, but if Diag
-    is :Unit then the diagonal elements of the matrix A are taken
-    as unity and are not referenced.")
-
-(defmfun matrix-vector-product-symmetric
-    (uplo alpha (A matrix) (x vector) beta (y vector))
+   "If the second argument is a vector, compute
+   inv(op(A)) x for x, where op(A) = A, A^T, A^H for
+   TransA = :NoTrans, :Trans, :ConjTrans. When Uplo is
+   :Upper then the upper triangle of A is used, and when Uplo is
+   :Lower then the lower triangle of A is used. If Diag is
+   :NonUnit then the diagonal of the matrix is used, but if Diag
+   is :Unit then the diagonal elements of the matrix A are taken
+   as unity and are not referenced.
+   If the second argument is a matrix, compute
+   the inverse-matrix matrix product B = alpha op(inv(A))B if
+   Side is :Left and B = alpha B op(inv(A)) if Side is
+   :Right. The matrix A is triangular and op(A) = A, A^T, A^H
+   for TransA = :NoTrans, :Trans, :ConjTrans When
+   Uplo is :Upper then the upper triangle of A is used, and
+   when Uplo is :Lower then the lower triangle of A is
+   used. If Diag is :NonUnit then the diagonal of A is used,
+   but if Diag is :Unit then the diagonal elements of the
+   matrix A are taken as unity and are not referenced.")
+
+(defmfun matrix-product-symmetric
+    ((A matrix) (x vector) (y vector)
+     &optional (alpha 1) (beta 1) (uplo :upper) (side :left))
   ("gsl_blas_" :type "symv")
   ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer A) :pointer)
    ((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer))
@@ -85,16 +122,26 @@
   :element-types :float
   :inputs (A x y)
   :outputs (y)
+  :return (y)
   :documentation			; FDL
-  "The matrix-vector product and sum y = alpha A
+  "If the second and third arguments are vectors, compute
+  the matrix-vector product and sum y = alpha A
   x + beta y for the symmetric matrix A. Since the matrix A is
   symmetric only its upper half or lower half need to be
   stored. When Uplo is :Upper then the upper triangle and
   diagonal of A are used, and when Uplo is :Lower then the
-  lower triangle and diagonal of A are used.")
+  lower triangle and diagonal of A are used.
+  If the second and third arguments are matrices, compute
+  the matrix-matrix product and sum C = alpha A
+  B + beta C for Side is :Left and C = alpha B A + beta C
+  for Side is :Right, where the matrix A is symmetric. When
+  Uplo is :Upper then the upper triangle and diagonal of A
+  are used, and when Uplo is :Lower then the lower triangle
+  and diagonal of A are used.")
 
-(defmfun matrix-vector-product-hermitian 
-  (uplo alpha (A matrix) (x vector) beta (y vector))
+(defmfun matrix-product-hermitian
+    ((A matrix) (x vector) (y vector)
+     &optional (alpha 1) (beta 1) (uplo :upper) (side :left))
   ("gsl_blas_" :type "hemv")
   ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer A) :pointer)
    ((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer))
@@ -102,14 +149,23 @@
   :element-types :complex
   :inputs (A x y)
   :outputs (y)
+  :return (y)
   :documentation			; FDL
-  "The matrix-vector product and sum y = alpha A x + beta y for the
+  "If the second and third arguments are vectors, compute the
+  matrix-vector product and sum y = alpha A x + beta y for the
   hermitian matrix A. Since the matrix A is hermitian only its upper
-  half or lower half need to be stored. When Uplo is CblasUpper then
+  half or lower half need to be stored. When Uplo is :upper then
   the upper triangle and diagonal of A are used, and when Uplo is
-  CblasLower then the lower triangle and diagonal of A are used. The
+  :lower then the lower triangle and diagonal of A are used. The
   imaginary elements of the diagonal are automatically assumed to be
-  zero and are not referenced.")
+  zero and are not referenced.  If the second and third arguments are
+  matrices, compute the matrix-matrix product and sum C = alpha A B +
+  beta C if Side is :left and C = \alpha B A + \beta C if Side
+  is :right, where the matrix A is hermitian. When Uplo is
+  :upper then the upper triangle and diagonal of A are used, and
+  when Uplo is :lower then the lower triangle and diagonal of A
+  are used. The imaginary elements of the diagonal are automatically
+  set to zero.")
 
 (defmfun rank-1-update (alpha (x vector) (y vector) (A matrix))
   ("gsl_blas_" :type "ger" :suffix)
@@ -119,6 +175,7 @@
   :element-types :float-complex
   :inputs (x y A)
   :outputs (A)
+  :return (A)
   :documentation			; FDL
    "The rank-1 update A = alpha x y^T + A of the matrix A.")
 
@@ -130,10 +187,13 @@
   :element-types :complex
   :inputs (x y A)
   :outputs (A)
+  :return (A)
   :documentation			; FDL
-   "The conjugate rank-1 update A = alpha x y^H + A of the matrix A.")
+  "The conjugate rank-1 update A = alpha x y^H + A of the matrix A.")
 
-(defmfun symmetric-rank-1-update (uplo alpha (x vector) (A matrix))
+(defmfun symmetric-rank-1-update
+    ((x vector) (A matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "syr")
   ((uplo cblas-uplo) (alpha :element-c-type)
    ((mpointer x) :pointer) ((mpointer A) :pointer))
@@ -141,15 +201,25 @@
   :element-types :float
   :inputs (x A)
   :outputs (A)
+  :return (A)
   :documentation			; FDL
-  "The symmetric rank-1 update A = \alpha x x^T +
-  A of the symmetric matrix A. Since the matrix A is symmetric
-  only its upper half or lower half need to be stored. When Uplo
-  is :Upper then the upper triangle and diagonal of A are
-  used, and when Uplo is :Lower then the lower triangle and
-  diagonal of A are used.")
-
-(defmfun hermitian-rank-1-update (uplo alpha (x vector) (A matrix))
+  "If the first argument is a vector,
+  the symmetric rank-1 update A = \alpha x x^T + A of the symmetric
+  matrix A. Since the matrix A is symmetric only its upper half or
+  lower half need to be stored. When Uplo is :Upper then the upper
+  triangle and diagonal of A are used, and when Uplo is :Lower then
+  the lower triangle and diagonal of A are used.  If the first
+  argument is a matrix, a rank-k update of the symmetric matrix C, C =
+  \alpha A A^T + \beta C when Trans is CblasNoTrans and C = \alpha A^T
+  A + \beta C when Trans is CblasTrans. Since the matrix C is
+  symmetric only its upper half or lower half need to be stored. When
+  Uplo is CblasUpper then the upper triangle and diagonal of C are
+  used, and when Uplo is CblasLower then the lower triangle and
+  diagonal of C are used.")
+
+(defmfun hermitian-rank-1-update
+    ((x vector) (A matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "her")
   ((uplo cblas-uplo) (alpha :element-c-type)
    ((mpointer x) :pointer) ((mpointer A) :pointer))
@@ -157,15 +227,27 @@
   :element-types :complex
   :inputs (x A)
   :outputs (A)
+  :return (A)
   :documentation			; FDL
-  "Compute the hermitian rank-1 update A = alpha x x^H + A of the
+  "If the first argument is a vector,
+  compute the hermitian rank-1 update A = alpha x x^H + A of the
   hermitian matrix A. Since the matrix A is hermitian only its upper
-  half or lower half need to be stored. When Uplo is :upper then
-  the upper triangle and diagonal of A are used, and when Uplo is
+  half or lower half need to be stored. When Uplo is :upper then the
+  upper triangle and diagonal of A are used, and when Uplo is
   :lower then the lower triangle and diagonal of A are used. The
-  imaginary elements of the diagonal are automatically set to zero.")
+  imaginary elements of the diagonal are automatically set to zero.
+  If the first argument is a matrix, compute a rank-k update of the
+  hermitian matrix C, C = \alpha A A^H + \beta C when Trans is
+  :notrans and C = \alpha A^H A + \beta C when Trans is
+  :trans. Since the matrix C is hermitian only its upper half or
+  lower half need to be stored. When Uplo is :upper then the upper
+  triangle and diagonal of C are used, and when Uplo is :lower
+  then the lower triangle and diagonal of C are used. The imaginary
+  elements of the diagonal are automatically set to zero.")
 
-(defmfun symmetric-rank-2-update (uplo alpha (x vector) (y vector) (A matrix))
+(defmfun symmetric-rank-2-update
+    ((x vector) (y vector) (A matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "syr2")
   ((uplo cblas-uplo) (alpha :element-c-type)
    ((mpointer x) :pointer)  ((mpointer y) :pointer)
@@ -174,15 +256,25 @@
   :element-types :float
   :inputs (x y A)
   :outputs (A)
+  :return (A)
   :documentation			; FDL
-  "The symmetric rank-2 update A = alpha x y^T +
-  alpha y x^T + A of the symmetric matrix A. Since the matrix A
-  is symmetric only its upper half or lower half need to be
-  stored. When Uplo is :Upper then the upper triangle and
-  diagonal of A are used, and when Uplo is :Lower then the
-  lower triangle and diagonal of A are used.")
+  "If the first two arguments are vectors, compute
+  the symmetric rank-2 update A = alpha x y^T + alpha y x^T + A of the
+  symmetric matrix A. Since the matrix A is symmetric only its upper
+  half or lower half need to be stored. When Uplo is :upper then the
+  upper triangle and diagonal of A are used, and when Uplo is :lower
+  then the lower triangle and diagonal of A are used.  If the first
+  two arguments are matrices, compute a rank-2k update of the
+  symmetric matrix C, C = \alpha A B^T + \alpha B A^T + \beta C when
+  Trans is :notrans and C = \alpha A^T B + \alpha B^T A + \beta C
+  when Trans is :trans. Since the matrix C is symmetric only its
+  upper half or lower half need to be stored. When Uplo is :upper
+  then the upper triangle and diagonal of C are used, and when Uplo is
+  :lower then the lower triangle and diagonal of C are used.")
 
-(defmfun hermitian-rank-2-update (uplo alpha (x vector) (y vector) (A matrix))
+(defmfun hermitian-rank-2-update
+    ((x vector) (y vector) (A matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "her2")
   ((uplo cblas-uplo) (alpha :element-c-type)
    ((mpointer x) :pointer) ((mpointer A) :pointer))
@@ -190,13 +282,24 @@
   :element-types :complex
   :inputs (x A)
   :outputs (A)
+  :return (A)
   :documentation			; FDL
-  "The hermitian rank-2 update A = alpha x y^H + alpha^* y x^H A of
+  "If the first two arguments are vectors, compute the
+  hermitian rank-2 update A = alpha x y^H + alpha^* y x^H A of
   the hermitian matrix A. Since the matrix A is hermitian only its
   upper half or lower half need to be stored. When uplo is :upper
   then the upper triangle and diagonal of A are used, and when uplo is
   :lower then the lower triangle and diagonal of A are used. The
-  imaginary elements of the diagonal are automatically set to zero.")
+  imaginary elements of the diagonal are automatically set to zero.
+  If the first two arguments are matrices, compute a rank-2k update of
+  the hermitian matrix C, C = \alpha A B^H + \alpha^* B A^H + \beta C
+  when Trans is :notrans and C = \alpha A^H B + \alpha^* B^H A +
+  \beta C when Trans is :conjtrans. Since the matrix C is
+  hermitian only its upper half or lower half need to be stored. When
+  Uplo is :upper then the upper triangle and diagonal of C are
+  used, and when Uplo is :lower then the lower triangle and
+  diagonal of C are used. The imaginary elements of the diagonal are
+  automatically set to zero.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
diff --git a/linear-algebra/blas3.lisp b/linear-algebra/blas3.lisp
index 32aeb5497d3447bf95ddde238e93c18361427195..4d733a8df4fe30c59a0635db552ed42dc8076d5e 100644
--- a/linear-algebra/blas3.lisp
+++ b/linear-algebra/blas3.lisp
@@ -1,6 +1,6 @@
 ;; BLAS level 3, Matrix-matrix operations
 ;; Liam Healy, Wed Apr 26 2006 - 21:08
-;; Time-stamp: <2008-03-09 14:13:42EDT blas3.lisp>
+;; Time-stamp: <2008-08-09 19:46:23EDT blas3.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -14,185 +14,127 @@
   (:left 141) :right)
 
 ;;;;****************************************************************************
-;;;; Generic
+;;;; Functions
 ;;;;****************************************************************************
 
-(export '(gemm symm trmm trsm syrk syr2k))
-
-(defgeneric gemm (TransA TransB alpha A B beta C)
-  (:method :after (TransA TransB alpha A B beta C)
-	   (cl-invalidate C))
-  (:documentation			; FDL
-   "The matrix-matrix product and sum C = alpha
-  op(A) op(B) + beta C where op(A) = A, A^T, A^H for TransA =
-  :NoTrans, :Trans, :ConjTrans and similarly for the
-  parameter TransB."))
-
-(defgeneric symm (side uplo alpha A B beta C)
-  (:method :after (side uplo alpha A B beta C)
-	   (cl-invalidate C))
-  (:documentation			; FDL
-   "The matrix-matrix product and sum C = alpha A
-  B + beta C for Side is :Left and C = alpha B A + beta C
-  for Side is :Right, where the matrix A is symmetric. When
-  Uplo is :Upper then the upper triangle and diagonal of A
-  are used, and when Uplo is :Lower then the lower triangle
-  and diagonal of A are used."))
-
-(defgeneric trmm (side uplo TransA diag alpha A B)
-  (:method :after (side uplo TransA diag alpha A B)
-	   (cl-invalidate B))
-  (:documentation			; FDL
-   "The matrix-matrix product B = \alpha op(A) B
-  for Side is :Left and B = \alpha B op(A) for Side is
-  :Right. The matrix A is triangular and op(A) = A, A^T, A^H
-  for TransA = :NoTrans, :Trans, :ConjTrans When Uplo
-  is :Upper then the upper triangle of A is used, and when
-  Uplo is :Lower then the lower triangle of A is used. If
-  Diag is :NonUnit then the diagonal of A is used, but if
-  Diag is :Unit then the diagonal elements of the matrix A
-  are taken as unity and are not referenced."))
-
-(defgeneric trsm (side uplo TransA diag alpha A B)
-  (:method :after (side uplo TransA diag alpha A B)
-	   (cl-invalidate B))
-  (:documentation
-   "The inverse-matrix matrix product B = \alpha op(inv(A))B for
-   Side is :Left and B = \alpha B op(inv(A)) for Side is
-   :Right. The matrix A is triangular and op(A) = A, A^T, A^H
-   for TransA = :NoTrans, :Trans, :ConjTrans When
-   Uplo is :Upper then the upper triangle of A is used, and
-   when Uplo is :Lower then the lower triangle of A is
-   used. If Diag is :NonUnit then the diagonal of A is used,
-   but if Diag is :Unit then the diagonal elements of the
-   matrix A are taken as unity and are not referenced."))
-
-(defgeneric syrk (uplo trans alpha A beta C)
-  (:method :after (uplo trans alpha A beta C)
-	   (cl-invalidate C))
-  (:documentation "A rank-k update of the symmetric matrix C, C =
-  alpha A A^T + beta C when Trans is :NoTrans and C =
-  alpha A^T A + beta C when Trans is :Trans. Since the
-  matrix C is symmetric only its upper half or lower half need to
-  be stored. When Uplo is :Upper then the upper triangle and
-  diagonal of C are used, and when Uplo is :Lower then the
-  lower triangle and diagonal of C are used."))
-
-(defgeneric syr2k (uplo trans alpha A B beta C)
-  (:method :after (uplo trans alpha A B beta C)
-	   (cl-invalidate C))
-  (:documentation "A rank-2k update of the symmetric matrix C, C
-  = alpha A B^T + alpha B A^T + beta C when Trans is
-  :NoTrans and C = alpha A^T B + alpha B^T A + \beta C when
-  Trans is :Trans. Since the matrix C is symmetric only its
-  upper half or lower half need to be stored. When Uplo is
-  :Upper then the upper triangle and diagonal of C are used,
-  and when Uplo is :Lower then the lower triangle and
-  diagonal of C are used."))
-
-;;;;****************************************************************************
-;;;; Single
-;;;;****************************************************************************
-
-(defmfun gemm
-    (TransA TransB alpha (A matrix-single-float) (B matrix-single-float)
-	    beta (C matrix-single-float))
-  "gsl_blas_sgemm"
-  ((transa cblas-transpose) (transb cblas-transpose)
-   (alpha :float) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)
-   (beta :float) ((pointer C) gsl-matrix-c))
-  :type :method)
-
-(defmfun symm
-    (side uplo alpha (A matrix-single-float) (B matrix-single-float)
-	  beta (C matrix-single-float))
-  "gsl_blas_ssymm"
-  ((side cblas-side) (uplo cblas-uplo) (alpha :float)
-   ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)
-   (beta :float) ((pointer C) gsl-matrix-c))
-  :type :method)
-
-(defmfun trmm
-    (side uplo transa diag alpha (A matrix-single-float) (B matrix-single-float))
-  "gsl_blas_strmm"
-  ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
-   (alpha :float) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c))
-  :type :method)
-
-(defmfun trsm
-    (side uplo transa diag alpha (A matrix-single-float) (B matrix-single-float))
-  "gsl_blas_strsm"
-  ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
-   (alpha :float) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c))
-  :type :method)
-
-(defmfun syrk
-    (uplo trans alpha (A matrix-single-float) beta (C matrix-single-float))
-  "gsl_blas_ssyrk"
-  ((uplo cblas-uplo) (trans cblas-transpose) (alpha :float)
-   ((pointer A) gsl-matrix-c) (beta :float) ((pointer C) gsl-matrix-c))
-  :type :method)
-
-(defmfun syr2k
-    (uplo trans alpha (A matrix-single-float) (B matrix-single-float)
-	  beta (C matrix-single-float))
-  "gsl_blas_ssyr2k"
-  ((uplo cblas-uplo) (trans cblas-transpose) (alpha :float)
-   ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)
-   (beta :float) ((pointer C) gsl-matrix-c))
-  :type :method)
-
-;;;;****************************************************************************
-;;;; Double
-;;;;****************************************************************************
-
-(defmfun gemm
-    (TransA TransB alpha (A matrix-double-float) (B matrix-double-float)
-	    beta (C matrix-double-float))
-  "gsl_blas_dgemm"
-  ((transa cblas-transpose) (transb cblas-transpose)
-   (alpha :double) (A gsl-matrix-c) (B gsl-matrix-c)
-   (beta :double) (C gsl-matrix-c))
-  :type :method)
-
-(defmfun symm
-    (side uplo alpha (A matrix-double-float) (B matrix-double-float)
-	  beta (C matrix-double-float))
-  "gsl_blas_dsymm"
-  ((side cblas-side) (uplo cblas-uplo) (alpha :double)
-   ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)
-   (beta :double) ((pointer C) gsl-matrix-c))
-  :type :method)
-
-(defmfun trmm
-    (side uplo transa diag alpha (A matrix-double-float) (B matrix-double-float))
-  "gsl_blas_dtrmm"
-  ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
-   (alpha :double) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c))
-  :type :method)
-
-(defmfun trsm
-    (side uplo transa diag alpha (A matrix-double-float) (B matrix-double-float))
-  "gsl_blas_dtrsm"
-  ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
-   (alpha :double) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c))
-  :type :method)
-
-(defmfun syrk
-    (uplo trans alpha (A matrix-double-float) beta (C matrix-double-float))
-  "gsl_blas_dsyrk"
-  ((uplo cblas-uplo) (trans cblas-transpose) (alpha :double)
-   ((pointer A) gsl-matrix-c) (alpha :double) ((pointer C) gsl-matrix-c))
-  :type :method)
-
-(defmfun syr2k
-    (uplo trans (alpha float) (A matrix-double-float) (B matrix-double-float)
-	  beta (C matrix-double-float))
-  "gsl_blas_dsyr2k"
-  ((uplo cblas-uplo) (trans cblas-transpose) (alpha :double)
-   ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)
-   (beta :double) ((pointer C) gsl-matrix-c))
-  :type :method)
+(defmfun matrix-product
+    ((A matrix) (B matrix) (C matrix)
+     &optional (alpha 1) (beta 1) (TransA :notrans) (TransB :notrans))
+  ("gsl_blas_" :type "gemm")
+  ((TransA cblas-transpose) (TransB cblas-transpose)
+   (alpha :element-c-type) ((mpointer A) :pointer)
+   ((mpointer B) :pointer) (beta :element-c-type) ((mpointer C) :pointer))
+  :definition :methods
+  :element-types :float-complex
+  :inputs (A B C)
+  :outputs (C)
+  :return (C))
+
+(defmfun matrix-product-symmetric
+    ((A matrix) (B matrix) (C matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (side :left))
+  ("gsl_blas_" :type "symm")
+  ((side cblas-side) (uplo cblas-uplo) (alpha :element-c-type)
+   (A :pointer) (B :pointer) (beta :element-c-type) (C :pointer)) 
+  :definition :generic
+  :element-types :float-complex
+  :inputs (A B C)
+  :outputs (C)
+  :return (C))
+
+(defmfun matrix-product-hermitian
+    ((A matrix) (B matrix) (C matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (side :left))
+  ("gsl_blas_" :type "hemm")
+  ((side cblas-side) (uplo cblas-uplo) (alpha :element-c-type)
+   ((mpointer A) :pointer) ((mpointer B) :pointer)
+   (beta :element-c-type) ((mpointer C) :pointer))
+  :definition :methods
+  :element-types :complex
+  :inputs (A B C)
+  :outputs (C)
+  :return (C))
+
+(defmfun matrix-product-triangular
+    ((A matrix) (B matrix)
+     &optional (alpha 1) (uplo :upper) (TransA :notrans)
+     (diag :nonunit) (side :left))
+  ("gsl_blas_" :type "trmm")
+  ((side cblas-side) (uplo cblas-uplo) (TransA cblas-transpose)
+   (diag cblas-diag)
+   (alpha :element-c-type) ((mpointer A) :pointer) ((mpointer B) :pointer))
+  :definition :methods
+  :element-types :float-complex
+  :inputs (A B)
+  :outputs (B)
+  :return (B))
+
+(defmfun inverse-matrix-product
+    ((A matrix) (B matrix)
+     &optional (alpha 1) (uplo :upper) (TransA :notrans)
+     (diag :nonunit) (side :left))
+  ("gsl_blas_" :type "trsm")
+  ((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag)
+   ((mpointer A) :pointer) ((mpointer B) :pointer))
+  :definition :methods
+  :element-types :float-complex
+  :inputs (A B)
+  :outputs (B)
+  :return (B))
+
+(defmfun symmetric-rank-1-update
+    ((A matrix) (C matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
+  ("gsl_blas_" :type "syrk")
+  ((uplo cblas-uplo) (trans cblas-transpose)
+   (alpha :element-c-type) ((mpointer A) :pointer)
+   (beta :element-c-type) ((mpointer C) :pointer))
+  :definition :methods
+  :element-types :float-complex
+  :inputs (A C)
+  :outputs (C)
+  :return (C))
+
+(defmfun hermitian-rank-1-update
+    ((A matrix) (C matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
+  ("gsl_blas_" :type "herk")
+  ((uplo cblas-uplo) (trans cblas-transpose)
+   (alpha :element-c-type) ((mpointer A) :pointer)
+   (beta :element-c-type) ((mpointer C) :pointer))
+  :definition :methods
+  :element-types :complex
+  :inputs (A C)
+  :outputs (C)
+  :return (C))
+
+(defmfun symmetric-rank-2-update
+    ((A matrix) (B matrix) (C matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
+  ("gsl_blas_" :type "syr2k")
+  ((uplo cblas-uplo) (trans cblas-transpose)
+   (alpha :element-c-type) ((mpointer A) :pointer)
+   ((mpointer B) :pointer) (beta :element-c-type)
+   ((mpointer C) :pointer))
+  :definition :methods
+  :element-types :float
+  :inputs (A B C)
+  :outputs (C)
+  :return (C))
+
+(defmfun hermitian-rank-2-update
+    ((A matrix) (B matrix) (C matrix)
+     &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
+  ("gsl_blas_" :type "her2k")
+  ((uplo cblas-uplo) (trans cblas-transpose)
+   (alpha :element-c-type) ((mpointer A) :pointer)
+   ((mpointer B) :pointer) (beta :element-c-type)
+   ((mpointer C) :pointer))
+  :definition :methods
+  :element-types :complex
+  :inputs (A B C)
+  :outputs (C)
+  :return (C))
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test