matrix-product in blas3.lisp calls the gsl_blas gemm routine. It has optional parameters TransA and TransB to specify whether the matrices A and/or B should be transposed prior to multiplying. These two take values :trans or :notrans.

The bug is twofold. matrix-product calls matrix-product-dimensions (defined in blas2). To do its job properly, matrix-product-dimensions should have TransA and TransB as parameters, and matrix-product should call it with those as arguments.

Here is the modified matrix-product-dimensions:

(defun matrix-product-dimensions (a b &optional (TRANSA :NOTRANS) (TRANSB :NOTRANS))
(if (typep b 'grid:matrix)
      (list (funcall (case transa
               (:notrans #'first)
               (:trans #'second)
               (t (error "Invalid tranposition keyword, ~a" transa)))
             (grid:dimensions a))
        (funcall (case transb
               (:notrans #'second)
               (:trans #'first)
               (t (error "Invalid tranposition keyword, ~a" transb)))
             (grid:dimensions b)))
      (first (funcall (case transa
               (:notrans #'first)
               (:trans #'second)
               (t (error "Invalid tranposition keyword, ~a" transa)))
             (grid:dimensions a)))))

This can probably be improved. I applied the same fix to the case when b is not a matrix, since matrix-product in blas2 which works on matrix x vector, also accepts that TransB argument.

Finally matrix-product has to be fixed. It is defined in both blas2 (for matrix vector multiplication) and in blas3 (for matrix matrix multiplication). The modifications consist of modifying the call to matrix-product-dimensions to include TransA (and TransB)

;; blas2

(defmfun matrix-product
    ((A grid:matrix) (x vector)
     &optional
     y
     (alpha 1) (beta 1) (TransA :notrans) TransB
     &aux
     (yarr
      (grid:ensure-foreign-array
       y (matrix-product-dimensions A x TransA) element-type 0)))
("gsl_blas_" :type "gemv")
((transa cblas-transpose) (alpha :element-c-type) ((mpointer A) :pointer)
   ((mpointer x) :pointer) (beta :element-c-type) ((mpointer yarr) :pointer))
:definition :generic
:element-types #+fsbv :float-complex #-fsbv :float
:inputs (A x)
:outputs (yarr)
:documentation            ; FDL
"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.
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.")

;; blas3
(defmfun matrix-product
    ((A grid:matrix) (B grid:matrix)
     &optional
     C
     (alpha 1) (beta 1) (TransA :notrans) (TransB :notrans)
     &aux
     (Carr
      (or C
    (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A B TransA TransB)
               :initial-element 0))))
("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 Carr) :pointer))
:definition :methods
:element-types #+fsbv :float-complex #-fsbv :float
:inputs (A B Carr)
:outputs (Carr))

I did not test this exhaustively, but the modifications did work my a few matrix multiplications that I needed.

Mirko