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.vuko...@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 >