transposed matrix multiplication (in blas2,3) bugfixes

[Mirko Vukovic]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.common-lisp.net/pipermail/gsll-devel/attachments/20130717/23f79dfc/attachment.html>