On 11/30/10 4:00 PM, Anton Voropaev wrote:
> Exactly, I intended a speed-up by means of using hardware
> floating-point numbers. (While integers are more desirable.)
>
> I have found the dgemm function in the lapack distribution.
> How do I properly invoke it. ?dgemm does not seem to work.
>
Save the lisp routine below to some file and in maxima do something like
:lisp (load (compile-file "dgemm.lisp"))
(The above is done, after loading lapack!)
In maxima, define the function:
dgemm(a, b, [options]) :=
%dgemm(a, b, options);
With change, I get the following result with current CVS on my mac,
using cmucl:
(%i10) za:genmatrix(lambda([i,j],0),1000,1000)$
Evaluation took 8.8400 seconds (8.9700 elapsed) using 267.106 MB.
(%i11) za . za $
Evaluation took 40.4700 seconds (43.2100 elapsed) using 45.899 MB.
(%i12) dgemm(za,za)$
Evaluation took 1.4700 seconds (1.4900 elapsed) using 145.174 MB.
So 40 sec to multiply with maxima, 1.5 sec with lapack.
> By the way, load("lapack") does not finish in wxMaxima.
> I use Maxima 5.19.2 and Maxima 5.22.1 in Windows XP.
It takes quite some time to compile up all of lapack. If that is not
the problem you see, then you'll have to provide more details.
Ray
dgemm.lisp:
(in-package :maxima)
(defun %%dgemm (a b &key c transpose_a transpose_b (alpha 1d0) (beta 0d0))
(multiple-value-bind (a-nrows a-ncols)
(maxima-matrix-dims a)
(multiple-value-bind (b-nrows b-ncols)
(maxima-matrix-dims b)
(let ((matrix-a (lapack-lispify-matrix a a-nrows a-ncols))
(matrix-b (lapack-lispify-matrix a a-nrows a-ncols))
(matrix-c (cond ((and c (not (zerop beta)))
(lapack-lispify-matrix a a-nrows b-ncols))
(t
;; No C matrix given, or beta is zero.
;; Force beta to be zero to tell LAPACK
;; not to add C. But we still need to
;; create a matrix.
(setf beta 0d0)
(make-array (* a-nrows b-ncols) :element-type
'double-float))))
(trans-a (if transpose_a "t" "n"))
(trans-b (if transpose_b "t" "n")))
(blas::dgemm trans-a trans-b
a-nrows b-ncols a-ncols
alpha
matrix-a a-nrows
matrix-b b-nrows
beta
matrix-c a-nrows)
;; matrix-c contains the desired result.
(lapack-maxify-matrix a-nrows b-ncols matrix-c)))))
(defun $%dgemm (a b options)
(let* ((args (lispify-maxima-keyword-options
(cdr options)
'($c $transpose_a $transpose_b $alpha $beta))))
(apply #'%%dgemm a b args)))