(setf eps 1e-10)
(setf dx (- (iseq 1 10) 5))
(setf dy (* dx dx))
(setf zx (complex dx dy))
(setf zy (complex dy dx))
;(defun test-report (name val) (format t "~&~a: ~9t~f~%" name val))
(defun test-report (name val) (check #'< val eps))
(defun bv (n x incx)
(if (< incx 0)
(reverse (select x (* (abs incx) (iseq n))))
(select x (* incx (iseq n)))))
;;;;
;;;; _ASUM
;;;;
(defun asum (n x incx)
(let ((x (bv n x incx)))
(sum (abs (realpart x)) (abs (imagpart x)))))
(let ((x (coerce dx '(vector c-double)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(let ((val (blas-dasum n x 0 incx))
(v (asum n x incx)))
(setf maxdiff (max maxdiff (abs (- val v))))))
(test-report "DASUM" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(let ((val (blas-dzasum n x 0 incx))
(v (asum n x incx)))
(setf maxdiff (max maxdiff (abs (- val v))))))
(test-report "DZASUM" maxdiff))
;;;;
;;;; _AXPY
;;;;
(defun axpy (n a x incx y incy)
(let ((x (bv n x incx))
(y (bv n y incy)))
(+ (* a x) y)))
(let ((x (coerce dx '(vector c-double)))
(y (coerce dy '(vector c-double)))
(a 3.0)
(n 5)
(val nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf val (copy-array y))
(blas-daxpy n a x 0 incx val 0 incy)
(let ((diff (max (abs (- (bv n val incy) (axpy n a x incx y incy))))))
(setf maxdiff (max maxdiff diff)))))
(test-report "DAXPY" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 2 dy)) '(vector c-dcomplex)))
(a (complex 3.0 6.0))
(n 5)
(val nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf val (copy-array y))
(blas-zaxpy n a x 0 incx val 0 incy)
(let ((diff (max (abs (- (bv n val incy) (axpy n a x incx y incy))))))
(setf maxdiff (max maxdiff diff)))))
(test-report "ZAXPY" maxdiff))
;;;;
;;;; _DOT_
;;;;
(defun dot (n x incx y incy &optional (conj nil))
(let ((x (bv n x incx))
(y (bv n y incy)))
(sum (* (if conj (conjugate x) x) y))))
(let ((x (coerce dx '(vector c-double)))
(y (coerce dy '(vector c-double)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(let ((val (blas-ddot n x 0 incx y 0 incy))
(v (dot n x incx y incy)))
(setf maxdiff (max maxdiff (abs (- val v)))))))
(test-report "DDOT" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 2 dy)) '(vector c-dcomplex)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(let ((val (blas-zdotu n x 0 incx y 0 incy))
(v (dot n x incx y incy)))
(setf maxdiff (max maxdiff (abs (- val v)))))))
(test-report "ZDOTU" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 3 dy)) '(vector c-dcomplex)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(let ((val (blas-zdotc n x 0 incx y 0 incy))
(v (dot n x incx y incy t)))
(setf maxdiff (max maxdiff (abs (- val v)))))))
(test-report "ZDOTC" maxdiff))
;;;;
;;;; _NRM2
;;;;
(defun nrm2(n x incx)
(let ((x (bv n x incx)))
(sqrt (abs (sum (* (conjugate x) x))))))
(let ((x (coerce dx '(vector c-double)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(let ((val (blas-dnrm2 n x 0 incx))
(v (nrm2 n x incx)))
(setf maxdiff (max maxdiff (abs (- val v))))))
(test-report "DNRM2" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(let ((val (blas-dznrm2 n x 0 incx))
(v (nrm2 n x incx)))
(setf maxdiff (max maxdiff (abs (- val v))))))
(test-report "DZNRM2" maxdiff))
;;;;
;;;; _ROT
;;;;
(defun rot (n x incx y incy c s)
(let ((x (bv n x incx))
(y (bv n y incy)))
(values (+ (* c x) (* s y))
(- (* c y) (* s x)))))
(let* ((x (coerce dx '(vector c-double)))
(y (coerce dy '(vector c-double)))
(c .2)
(s (sqrt (- 1 (* c c))))
(n 5)
(valx nil)
(valy nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf valx (copy-array x))
(setf valy (copy-array y))
(blas-drot n valx 0 incx valy 0 incy c s)
(multiple-value-bind
(vx vy)
(rot n x incx y incy c s)
(let ((vvx (bv n valx incx))
(vvy (bv n valy incy)))
(setf maxdiff (max maxdiff (abs (- vvx vx)) (abs (- vvy vy))))))))
(test-report "DROT" maxdiff))
(let* ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 3 dy)) '(vector c-dcomplex)))
(c .2)
(s (sqrt (- 1 (* c c))))
(n 5)
(valx nil)
(valy nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf valx (copy-array x))
(setf valy (copy-array y))
(blas-zdrot n valx 0 incx valy 0 incy c s)
(multiple-value-bind
(vx vy)
(rot n x incx y incy c s)
(let ((vvx (bv n valx incx))
(vvy (bv n valy incy)))
(setf maxdiff (max maxdiff (abs (- vvx vx)) (abs (- vvy vy))))))))
(test-report "ZDROT" maxdiff))
;;;;
;;;; _ROTG
;;;;
(defun drotg (da db)
(let ((roe (if (> (abs da) (abs db)) da db))
(scale (+ (abs da) (abs db))))
(if (= scale 0.0)
(values 0.0 0.0 1.0 0.0)
(let* ((r (* scale (sqrt (+ (^ (/ da scale) 2) (^ (/ db scale) 2)))))
(r (* (signum roe) r))
(c (/ da r))
(s (/ db r))
(z 1.0))
(if (> (abs da) (abs db)) (seft z s))
(if (and (>= (abs db) (abs da)) (/= c 0.0)) (setf z (/ 1.0 c)))
(values r z c s)))))
(defun zrotg (ca cb)
(if (= (abs ca) 0.0)
(values cb cb 0.0 1.0)
(let* ((scale (+ (abs ca) (abs cb)))
(norm (* scale (sqrt (+ (^ (abs (/ ca scale)) 2)
(^ (abs (/ cb scale)) 2)))))
(alpha (/ ca (abs ca)))
(c (/ (abs ca) norm))
(s (/ (* alpha (conjugate cb)) norm))
(ca (* alpha norm)))
(values ca cb c s))))
(let ((a 1)
(b 7))
(test-report "DROTG"
(max (abs (- (multiple-value-list (blas-drotg a b))
(multiple-value-list (drotg a b)))))))
(let ((a #c(1 2))
(b #c(7 3)))
(test-report "ZROTG"
(max (abs (- (multiple-value-list (blas-zrotg a b))
(multiple-value-list (zrotg a b)))))))
;;;;
;;;; _SCAL
;;;;
(defun scal (n a x incx)
(let ((x (bv n x incx)))
(* a x)))
(let ((x (coerce dx '(vector c-double)))
(a 3.0)
(n 5)
(val nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(setf val (copy-array x))
(blas-dscal n a val 0 incx)
(let ((v (scal n a x incx))
(vv (bv n val incx)))
(setf maxdiff (max maxdiff (max (abs (- vv v)))))))
(test-report "DSCAL" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(a 3.0)
(n 5)
(val nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(setf val (copy-array x))
(blas-zdscal n a val 0 incx)
(let ((v (scal n a x incx))
(vv (bv n val incx)))
(setf maxdiff (max maxdiff (max (abs (- vv v)))))))
(test-report "ZDSCAL" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(a (complex 3.0 6.0))
(n 5)
(val nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(setf val (copy-array x))
(blas-zscal n a val 0 incx)
(let ((v (scal n a x incx))
(vv (bv n val incx)))
(setf maxdiff (max maxdiff (max (abs (- vv v)))))))
(test-report "ZSCAL" maxdiff))
;;;;
;;;; _SWAP
;;;;
(let ((x (coerce dx '(vector c-double)))
(y (coerce dy '(vector c-double)))
(n 5)
(valx nil)
(valy nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf valx (copy-array x))
(setf valy (copy-array y))
(blas-dswap n valx 0 incx valy 0 incy)
(let ((vx (bv n x incx))
(vy (bv n y incy))
(vvx (bv n valx incx))
(vvy (bv n valy incy)))
(setf maxdiff
(max maxdiff (max (abs (- vvx vy)) (abs (- vvy vx))))))))
(test-report "DSWAP" maxdiff))
(let ((x (coerce dx '(vector c-dcomplex)))
(y (coerce dy '(vector c-dcomplex)))
(n 5)
(valx nil)
(valy nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf valx (copy-array x))
(setf valy (copy-array y))
(blas-zswap n valx 0 incx valy 0 incy)
(let ((vx (bv n x incx))
(vy (bv n y incy))
(vvx (bv n valx incx))
(vvy (bv n valy incy)))
(setf maxdiff
(max maxdiff (max (abs (- vvx vy)) (abs (- vvy vx))))))))
(test-report "ZSWAP" maxdiff))
;;;;
;;;; I_AMAX
;;;;
(defun iamax (n x incx)
(let* ((x (bv n x incx))
(ax (+ (abs (realpart x)) (abs (imagpart x)))))
(position (max ax) ax)))
(let ((x (coerce dx '(vector c-double)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(let ((val (blas-idamax n x 0 incx))
(v (iamax n x incx)))
(setf maxdiff (max maxdiff (abs (- val v))))))
(test-report "IDAMAX" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(n 5)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(let ((val (blas-izamax n x 0 incx))
(v (iamax n x incx)))
(setf maxdiff (max maxdiff (abs (- val v))))))
(test-report "IZAMAX" maxdiff))
;;;;
;;;; _COPY
;;;;
(let ((x (coerce dx '(vector c-double)))
(y (coerce dy '(vector c-double)))
(n 5)
(val nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf val (copy-array y))
(blas-dcopy n x 0 incx val 0 incy)
(let ((v (bv n x incx))
(vv (bv n val incy)))
(setf maxdiff (max maxdiff (max (abs (- v vv))))))))
(test-report "DCOPY" maxdiff))
(let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 2 dy)) '(vector c-dcomplex)))
(n 5)
(val nil)
(maxdiff 0))
(dolist (incx '(1 2 -1 -2))
(dolist (incy '(1 2 -1 -2))
(setf val (copy-array y))
(blas-zcopy n x 0 incx val 0 incy)
(let ((v (bv n x incx))
(vv (bv n val incy)))
(setf maxdiff (max maxdiff (max (abs (- v vv))))))))
(test-report "ZCOPY" maxdiff))
(setf da '(1 2 3 4 5 6))
(setf db '(1 2 3 4 5 6 13 14 15))
(setf dc '(16 17 18 19 20 21))
(setf dx '(7 8 9))
(setf dy '(10 11 12))
(defvar *verbose-tests* nil)
;(defun test-report (name val) (format t "~&~a: ~8t~f~%" name val))
;;;;
;;;; _GEMV
;;;;
(defun gemv (trans alpha a x incx beta y incy)
(when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
(when (/= (abs incy) 1) (error "increment has to be +1 or -1"))
(when (equalp trans "n") (setf a (transpose a)))
(when (equalp trans "c") (setf a (conjugate a)))
(let* ((m (array-dimension a 0))
(n (array-dimension a 1))
(z (make-array m :initial-element 0)))
(when (minusp incx) (setf x (reverse (subseq x 0 n))))
(when (minusp incy) (setf y (reverse (subseq y 0 m))))
(replace z (* beta y))
(dotimes (i m)
(dotimes (j n)
(incf (aref z i) (* alpha (aref a i j) (aref x j)))))
(if (minusp incy) (reverse z) z)))
(let ((a (coerce da '(vector c-double)))
(x (coerce dx '(vector c-double)))
(y (coerce dy '(vector c-double)))
(val nil)
(alpha 2)
(beta 3)
(maxdiff 0))
(dolist (trans '("n" "t"))
(dolist (incx '(1 -1))
(dolist (incy '(1 -1))
(setf val (copy-seq y))
(blas-dgemv trans 3 2 alpha a 0 3 x 0 incx beta val 0 incy)
(let ((v (gemv trans alpha (matrix '(2 3) a) x incx beta y incy))
(vv (if (equalp trans "n") val (subseq val 0 2))))
(when *verbose-tests*
(format t "DGEMV ~a ~2d ~2d: ~f~%"
trans
incx
incy
(max (abs (- vv v)))))
(setf maxdiff (max maxdiff (max (abs (- vv v)))))))))
(test-report "DGEMV" maxdiff))
(let ((a (coerce (complex da (* 2 da)) '(vector c-dcomplex)))
(x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 4 dy)) '(vector c-dcomplex)))
(val nil)
(alpha (complex 2 4))
(beta (complex 3 6))
(maxdiff 0))
(dolist (trans '("n" "t" "c"))
(dolist (incx '(1 -1))
(dolist (incy '(1 -1))
(setf val (copy-seq y))
(blas-zgemv trans 3 2 alpha a 0 3 x 0 incx beta val 0 incy)
(let ((v (gemv trans alpha (matrix '(2 3) a) x incx beta y incy))
(vv (if (equalp trans "n") val (subseq val 0 2))))
(when *verbose-tests*
(format t "ZGEMV ~a ~2d ~2d: ~f~%"
trans
incx
incy
(max (abs (- vv v)))))
(setf maxdiff (max maxdiff (max (abs (- vv v)))))))))
(test-report "ZGEMV" maxdiff))
;;;;
;;;; _GER_
;;;;
(defun ger (alpha x incx y incy a &optional conj)
(when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
(when (/= (abs incy) 1) (error "increment has to be +1 or -1"))
(let ((x (subseq x 0 (array-dimension a 1)))
(y (subseq y 0 (array-dimension a 0))))
(when (minusp incx) (setf x (reverse x)))
(when (minusp incy) (setf y (reverse y)))
(+ (* alpha
(outer-product (if conj (conjugate y) y) x))
a)))
(let ((a (coerce da '(vector c-double)))
(x (coerce dx '(vector c-double)))
(y (coerce dy '(vector c-double)))
(val nil)
(alpha 2)
(maxdiff 0))
(dolist (incx '(1 -1))
(dolist (incy '(1 -1))
(setf val (copy-seq a))
(blas-dger 3 2 alpha x 0 incx y 0 incy val 0 3)
(let ((v (ger alpha x incx y incy (matrix '(2 3) a))))
(when *verbose-tests*
(format t "DGER ~2d ~2d: ~f~%"
incx
incy
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))
(test-report "DGER" maxdiff))
(let ((a (coerce (complex da (* 2 da)) '(vector c-dcomplex)))
(x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 4 dy)) '(vector c-dcomplex)))
(val nil)
(alpha (complex 2 4))
(maxdiff 0))
(dolist (incx '(1 -1))
(dolist (incy '(1 -1))
(setf val (copy-seq a))
(blas-zgerc 3 2 alpha x 0 incx y 0 incy val 0 3)
(let ((v (ger alpha x incx y incy (matrix '(2 3) a) t)))
(when *verbose-tests*
(format t "ZGERC ~2d ~2d: ~f~%"
incx
incy
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))
(test-report "ZGERC" maxdiff))
(let ((a (coerce (complex da (* 2 da)) '(vector c-dcomplex)))
(x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
(y (coerce (complex dy (* 4 dy)) '(vector c-dcomplex)))
(val nil)
(alpha (complex 2 4))
(maxdiff 0))
(dolist (incx '(1 -1))
(dolist (incy '(1 -1))
(setf val (copy-seq a))
(blas-zgeru 3 2 alpha x 0 incx y 0 incy val 0 3)
(let ((v (ger alpha x incx y incy (matrix '(2 3) a))))
(when *verbose-tests*
(format t "ZGERU ~2d ~2d: ~f~%"
incx
incy
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))
(test-report "ZGERU" maxdiff))
;;;;
;;;; _TRMV
;;;;
(defun trmv (a x &optional upper transpose unit incx)
(when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
(cond
((equalp transpose "n")
(setf a (transpose a)))
(t (setf upper (if (equalp upper "u") "l" "u"))))
(when (equalp transpose "c") (setf a (conjugate a)))
(when (minusp incx) (setf x (reverse x)))
(let* ((n (array-dimension a 0))
(y (make-array n :initial-element 0)))
(when (not (equalp unit "n"))
(setf a (copy-array a))
(dotimes (i n) (setf (aref a i i) 1)))
(if (equalp upper "u")
(dotimes (i n)
(do ((j i (+ j 1)))
((>= j n))
(incf (aref y i) (* (aref a i j) (aref x j)))))
(dotimes (i n)
(dotimes (j (+ i 1))
(incf (aref y i) (* (aref a i j) (aref x j))))))
(if (minusp incx) (reverse y) y)))
(let ((b (coerce db '(vector c-double)))
(x (coerce dx '(vector c-double)))
(val nil)
(maxdiff 0))
(dolist (upper '("u" "l"))
(dolist (trans '("n" "t"))
(dolist (diag '("n" "u"))
(dolist (incx '(1 -1))
(setf val (copy-seq x))
(blas-dtrmv upper trans diag 3 b 0 3 val 0 incx)
(let ((v (trmv (matrix '(3 3) b) x upper trans diag incx)))
(when *verbose-tests*
(format t "DTRMV ~a ~a ~a ~2d: ~f~%"
upper
trans
diag
incx
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))))
(test-report "DTRMV" maxdiff))
(let ((b (coerce (complex db (* 2 db)) '(vector c-dcomplex)))
(x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
(val nil)
(maxdiff 0))
(dolist (upper '("u" "l"))
(dolist (trans '("n" "t" "c"))
(dolist (diag '("n" "u"))
(dolist (incx '(1 -1))
(setf val (copy-seq x))
(blas-ztrmv upper trans diag 3 b 0 3 val 0 incx)
(let ((v (trmv (matrix '(3 3) b) x upper trans diag incx)))
(when *verbose-tests*
(format t "ZTRMV ~a ~a ~a ~2d: ~f~%"
upper
trans
diag
incx
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))))
(test-report "ZTRMV" maxdiff))
(defun trsv (a x &optional upper transpose unit incx)
(when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
(cond
((equalp transpose "n")
(setf a (transpose a)))
(t (setf upper (if (equalp upper "u") "l" "u"))))
(when (equalp transpose "c") (setf a (conjugate a)))
(when (minusp incx) (setf x (reverse x)))
(let* ((n (array-dimension a 0))
(y (make-array n :initial-contents x)))
(if (equalp upper "u")
(do ((i (- n 1) (- i 1)))
((< i 0))
(do ((j (+ i 1) (+ j 1)))
((>= j n))
(decf (aref y i) (* (aref a i j) (aref y j))))
(when (equalp unit "n")
(setf (aref y i) (/ (aref y i) (aref a i i)))))
(dotimes (i n)
(dotimes (j i)
(decf (aref y i) (* (aref a i j) (aref y j))))
(when (equalp unit "n")
(setf (aref y i) (/ (aref y i) (aref a i i))))))
(if (minusp incx) (reverse y) y)))
(let ((b (coerce db '(vector c-double)))
(x (coerce dx '(vector c-double)))
(val nil)
(maxdiff 0))
(dolist (upper '("u" "l"))
(dolist (trans '("n" "t"))
(dolist (diag '("n" "u"))
(dolist (incx '(1 -1))
(setf val (copy-seq x))
(blas-dtrsv upper trans diag 3 b 0 3 val 0 incx)
(let ((v (trsv (matrix '(3 3) b) x upper trans diag incx)))
(when *verbose-tests*
(format t "DTRSV ~a ~a ~a ~2d: ~f~%"
upper
trans
diag
incx
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))))
(test-report "DTRSV" maxdiff))
(let ((b (coerce (complex db (* 2 db)) '(vector c-dcomplex)))
(x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
(val nil)
(maxdiff 0))
(dolist (upper '("u" "l"))
(dolist (trans '("n" "t"))
(dolist (diag '("n" "u"))
(dolist (incx '(1 -1))
(setf val (copy-seq x))
(blas-ztrsv upper trans diag 3 b 0 3 val 0 incx)
(let ((v (trsv (matrix '(3 3) b) x upper trans diag incx)))
(when *verbose-tests*
(format t "ZTRSV ~a ~a ~a ~2d: ~f~%"
upper
trans
diag
incx
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))))
(test-report "ZTRSV" maxdiff))
;;;;
;;;; _GEMM
;;;;
(defun gemm (transa transb alpha a b beta c)
(when (equalp transa "n") (setf a (transpose a)))
(when (equalp transa "c") (setf a (conjugate a)))
(when (equalp transb "n") (setf b (transpose b)))
(when (equalp transb "c") (setf b (conjugate b)))
(let* ((m (array-dimension c 1))
(n (array-dimension c 0))
(k (array-dimension a 1))
(z (* beta c)))
(dotimes (i m)
(dotimes (j n)
(dotimes (l k)
(incf (aref z j i) (* alpha (aref a i l) (aref b l j))))))
z))
(setf aa (transpose (matrix '(4 2) (iseq 1 8))))
(setf bb (transpose (matrix '(2 3) (iseq 9 14))))
(setf cc (transpose (matrix '(4 3) (iseq 15 26))))
(let ((a (coerce aa '(array c-double)))
(b (coerce bb '(array c-double)))
(c (coerce cc '(array c-double)))
(m 4)
(n 3)
(k 2)
(alpha 2)
(beta 3)
(val nil)
(maxdiff 0))
(dolist (transa '("n" "t"))
(dolist (transb '("n" "t"))
(let ((aa (if (equalp transa "n") a (transpose a)))
(bb (if (equalp transb "n") b (transpose b)))
(lda (if (equalp transa "n") m k))
(ldb (if (equalp transb "n") k n)))
(setf val (copy-array c))
(blas-dgemm transa transb m n k alpha aa 0 lda bb 0 ldb beta val 0 m)
(let ((v (gemm transa transb alpha aa bb beta c)))
(when *verbose-tests*
(format t "DGEMM ~a ~a: ~f~%"
transa
transb
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v)))))))))
(test-report "DGEMM" maxdiff))
(let ((a (coerce (complex aa (* 2 aa)) '(array c-dcomplex)))
(b (coerce (complex bb (* 3 bb)) '(array c-dcomplex)))
(c (coerce (complex cc (* 4 cc)) '(array c-dcomplex)))
(m 4)
(n 3)
(k 2)
(alpha #c(2 4))
(beta #c(3 6))
(val nil)
(maxdiff 0))
(dolist (transa '("n" "t" "c"))
(dolist (transb '("n" "t" "c"))
(let ((aa (if (equalp transa "n") a (transpose a)))
(bb (if (equalp transb "n") b (transpose b)))
(lda (if (equalp transa "n") m k))
(ldb (if (equalp transb "n") k n)))
(setf val (copy-array c))
(blas-zgemm transa transb m n k alpha aa 0 lda bb 0 ldb beta val 0 m)
(let ((v (gemm transa transb alpha aa bb beta c)))
(when *verbose-tests*
(format t "ZGEMM ~a ~a: ~f~%"
transa
transb
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v)))))))))
(test-report "ZGEMM" maxdiff))
;;;;
;;;; TRSM
;;;;
(defun trsm (side upper transpose unit a alpha x)
(cond
((equalp side "l")
(let ((cx (* alpha (row-list x))))
(apply #'bind-rows
(mapcar #'(lambda (c) (trsv a c upper transpose unit 1))
cx))))
(t
(let ((cx (* alpha (column-list x)))
(conj (equalp transpose "c")))
(setf transpose (if (equalp transpose "n") "t" "n"))
(when conj (setf cx (conjugate cx)))
(apply #'bind-columns
(mapcar #'(lambda (c)
(let ((v (trsv a c upper transpose unit 1)))
(if conj (conjugate v) v)))
cx))))))
(let ((b (coerce db '(vector c-double)))
(c (coerce dc '(vector c-double)))
(alpha 2)
(val nil)
(maxdiff 0.0))
(dolist (side '("l" "r"))
(dolist (uplo '("u" "l"))
(dolist (trans '("n" "t"))
(dolist (diag '("n" "u"))
(setf val (copy-array c))
(if (equalp side "l")
(blas-dtrsm side uplo trans diag 3 2 alpha b 0 3 val 0 3)
(blas-dtrsm side uplo trans diag 2 3 alpha b 0 3 val 0 2))
(let* ((mb (matrix '(3 3) b))
(mc (matrix (if (equalp side "l") '(2 3) '(3 2)) c))
(v (trsm side uplo trans diag mb alpha mc)))
(when *verbose-tests*
(format t "DTRSM ~a ~a ~a ~a: ~f~%"
side
uplo
trans
diag
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))))
(test-report "DTRSM" maxdiff))
(let ((b (coerce (complex db (* 2 db)) '(vector c-dcomplex)))
(c (coerce (complex dc (* 3 dc)) '(vector c-dcomplex)))
(alpha (complex 2 4))
(val nil)
(maxdiff 0))
(dolist (side '("l" "r"))
(dolist (uplo '("u" "l"))
(dolist (trans '("n" "t" "c"))
(dolist (diag '("n" "u"))
(setf val (copy-array c))
(if (equalp side "l")
(blas-ztrsm side uplo trans diag 3 2 alpha b 0 3 val 0 3)
(blas-ztrsm side uplo trans diag 2 3 alpha b 0 3 val 0 2))
(let* ((mb (matrix '(3 3) b))
(mc (matrix (if (equalp side "l") '(2 3) '(3 2)) c))
(v (trsm side uplo trans diag mb alpha mc)))
(when *verbose-tests*
(format t "ZTRSM ~a ~a ~a ~a: ~f~%"
side
uplo
trans
diag
(max (abs (- val v)))))
(setf maxdiff (max maxdiff (max (abs (- val v))))))))))
(test-report "ZTRSM" maxdiff))
(defun nax+y (a x y)
(let* ((m (array-dimension a 0))
(n (array-dimension a 1))
(da (coerce a '(array c-double)))
(dx (coerce x '(vector c-double)))
(dy (make-array m :element-type 'c-double :initial-contents y)))
(blas-dgemv "t" n m 1 da 0 n dx 0 1 1 dy 0 1)
(coerce dy (if (vectorp x) '(vector t) 'list))))
syntax highlighted by Code2HTML, v. 0.9.1