share/lapack/dgemm.lisp

   1 ;;; Simple raw interface to LAPACK dgemm, general real matrix
   2 ;;; multiplication.
   3 (in-package :maxima)
   4
   5 (defmfun $dgemm (a b &key (c nil cp) transpose_a transpose_b (alpha 1e0) (beta 0e0 betap))
   6   (flet ((maybe-transpose-dims (transp row col)
   7            (if transp
   8                (values col row)
   9                (values row col))))
  10     (multiple-value-bind (a-nrows a-ncols)
  11         (maxima-matrix-dims a)
  12       (multiple-value-bind (at-nrows at-ncols)
  13           (maybe-transpose-dims transpose_a a-nrows a-ncols)
  14         (multiple-value-bind (b-nrows b-ncols)
  15             (maxima-matrix-dims b)
  16           (multiple-value-bind (bt-nrows bt-ncols)
  17               (maybe-transpose-dims transpose_b b-nrows b-ncols)
  18             ;; Compute the dimensions of C.  If C is not given, we
  19             ;; need to create a space of the correct dimension to hold
  20             ;; the result.
  21             (multiple-value-bind (c-nrows c-ncols)
  22                 (if c
  23                     (maxima-matrix-dims c)
  24                     (values at-nrows bt-ncols))
  25               ;; Do something nice if C is given but beta is not.  We
  26               ;; take this to mean that we just want to add C,
  27               ;; implying that beta = 1.
  28               (when (and cp (not betap))
  29                 (setf beta 0e0))
  30               ;; Now for some error checking.  This is a bit redundant
  31               ;; since dgemm does some error checking, but I think we
  32               ;; prefer not to see messages from dgemm.  It's better
  33               ;; to have the messages come from maxima.
  34
  35               (unless (= at-ncols bt-nrows)
  36                 (merror "Cannot multiply a ~m by ~m matrix by a ~m by ~m matrix"
  37                         at-nrows at-ncols bt-nrows bt-ncols))
  38               (when (and c
  39                          (or (/= at-nrows c-nrows)
  40                              (/= bt-ncols c-ncols)))
  41                 (merror "Cannot add a ~m by ~m matrix to a ~m by ~m C matrix"
  42                         at-nrows bt-ncols c-nrows c-ncols))
  43
  44               ;; But it doesn't make sense to supply beta without C.
  45               ;; (Well, we could assume that C is zero, but then why
  46               ;; bother specifying beta?)
  47               (when (and betap (not cp))
  48                 (merror "beta given, but no C matrix supplied?"))
  49               (let ((alpha ($float alpha))
  50                     (beta ($float beta))
  51                     (matrix-a (lapack-lispify-matrix a a-nrows a-ncols))
  52                     (matrix-b (lapack-lispify-matrix b b-nrows b-ncols))
  53                     (matrix-c (cond ((and c (not (zerop beta)))
  54                                      (lapack-lispify-matrix c a-nrows b-ncols))
  55                                     (t
  56                                      ;; No C matrix given, or beta is zero.
  57                                      ;; Force beta to be zero to tell LAPACK
  58                                      ;; not to add C.  But we still need to
  59                                      ;; create a matrix.
  60                                      (setf beta 0e0)
  61                                      (make-array (* c-nrows c-ncols) :element-type 'flonum))))
  62                     (trans-a (if transpose_a "t" "n"))
  63                     (trans-b (if transpose_b "t" "n")))
  64                 (blas::dgemm trans-a trans-b
  65                              at-nrows bt-ncols at-ncols
  66                              alpha
  67                              matrix-a a-nrows
  68                              matrix-b b-nrows
  69                              beta
  70                              matrix-c c-nrows)
  71                 ;; matrix-c contains the desired result.
  72                 (lapack-maxify-matrix c-nrows c-ncols matrix-c)))))))))
  73
  74