transposed matrix multiplication (in blas2,3) bugfixes

[Liam Healy]

Thanks for the report, fixed in 1107fc49bc, currently on the
antik-multiple-systems branch but merged into master soon.

Liam


On Wed, Jul 17, 2013 at 2:35 PM, Mirko Vukovic <mirko.vukovic at gmail.com>wrote:

> 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/20131225/36a6abd3/attachment.html>