tests/blas.lsp

   1 (setf eps 1e-10)
   2 (setf dx (- (iseq 1 10) 5))
   3 (setf dy (* dx dx))
   4 (setf zx (complex dx dy))
   5 (setf zy (complex dy dx))
   6
   7 ;(defun test-report (name val) (format t "~&~a: ~9t~f~%" name val))
   8 (defun test-report (name val) (check #'< val eps))
   9
  10 (defun bv (n x incx)
  11   (if (< incx 0)
  12       (reverse (select x (* (abs incx) (iseq n))))
  13       (select x (* incx (iseq n)))))
  14
  15
  16 ;;;;
  17 ;;;; _ASUM
  18 ;;;;
  19
  20 (defun asum (n x incx)
  21   (let ((x (bv n x incx)))
  22     (sum (abs (realpart x)) (abs (imagpart x)))))
  23
  24 (let ((x (coerce dx '(vector c-double)))
  25       (n 5)
  26       (maxdiff 0))
  27   (dolist (incx '(1 2 -1 -2))
  28     (let ((val (blas-dasum n x 0 incx))
  29           (v (asum n x incx)))
  30       (setf maxdiff (max maxdiff (abs (- val v))))))
  31   (test-report "DASUM" maxdiff))
  32
  33 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
  34       (n 5)
  35       (maxdiff 0))
  36   (dolist (incx '(1 2 -1 -2))
  37     (let ((val (blas-dzasum n x 0 incx))
  38           (v (asum n x incx)))
  39       (setf maxdiff (max maxdiff (abs (- val v))))))
  40   (test-report "DZASUM" maxdiff))
  41
  42
  43 ;;;;
  44 ;;;; _AXPY
  45 ;;;;
  46
  47 (defun axpy (n a x incx y incy)
  48   (let ((x (bv n x incx))
  49         (y (bv n y incy)))
  50     (+ (* a x) y)))
  51
  52 (let ((x (coerce dx '(vector c-double)))
  53       (y (coerce dy '(vector c-double)))
  54       (a 3.0)
  55       (n 5)
  56       (val nil)
  57       (maxdiff 0))
  58   (dolist (incx '(1 2 -1 -2))
  59     (dolist (incy '(1 2 -1 -2))
  60       (setf val (copy-array y))
  61       (blas-daxpy n a x 0 incx val 0 incy)
  62       (let ((diff (max (abs (- (bv n val incy) (axpy n a x incx y incy))))))
  63         (setf maxdiff (max maxdiff diff)))))
  64   (test-report "DAXPY" maxdiff))
  65
  66 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
  67       (y (coerce (complex dy (* 2 dy)) '(vector c-dcomplex)))
  68       (a (complex 3.0 6.0))
  69       (n 5)
  70       (val nil)
  71       (maxdiff 0))
  72   (dolist (incx '(1 2 -1 -2))
  73     (dolist (incy '(1 2 -1 -2))
  74       (setf val (copy-array y))
  75       (blas-zaxpy n a x 0 incx val 0 incy)
  76       (let ((diff (max (abs (- (bv n val incy) (axpy n a x incx y incy))))))
  77         (setf maxdiff (max maxdiff diff)))))
  78   (test-report "ZAXPY" maxdiff))
  79
  80
  81 ;;;;
  82 ;;;; _DOT_
  83 ;;;;
  84
  85 (defun dot (n x incx y incy &optional (conj nil))
  86   (let ((x (bv n x incx))
  87         (y (bv n y incy)))
  88     (sum (* (if conj (conjugate x) x) y))))
  89
  90 (let ((x (coerce dx '(vector c-double)))
  91       (y (coerce dy '(vector c-double)))
  92       (n 5)
  93       (maxdiff 0))
  94   (dolist (incx '(1 2 -1 -2))
  95     (dolist (incy '(1 2 -1 -2))
  96       (let ((val (blas-ddot n x 0 incx y 0 incy))
  97             (v (dot n x incx y incy)))
  98         (setf maxdiff (max maxdiff (abs (- val v)))))))
  99   (test-report "DDOT" maxdiff))
 100
 101 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 102       (y (coerce (complex dy (* 2 dy)) '(vector c-dcomplex)))
 103       (n 5)
 104       (maxdiff 0))
 105   (dolist (incx '(1 2 -1 -2))
 106     (dolist (incy '(1 2 -1 -2))
 107       (let ((val (blas-zdotu n x 0 incx y 0 incy))
 108             (v (dot n x incx y incy)))
 109         (setf maxdiff (max maxdiff (abs (- val v)))))))
 110   (test-report "ZDOTU" maxdiff))
 111
 112 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 113       (y (coerce (complex dy (* 3 dy)) '(vector c-dcomplex)))
 114       (n 5)
 115       (maxdiff 0))
 116   (dolist (incx '(1 2 -1 -2))
 117     (dolist (incy '(1 2 -1 -2))
 118       (let ((val (blas-zdotc n x 0 incx y 0 incy))
 119             (v (dot n x incx y incy t)))
 120         (setf maxdiff (max maxdiff (abs (- val v)))))))
 121   (test-report "ZDOTC" maxdiff))
 122
 123 ;;;;
 124 ;;;; _NRM2
 125 ;;;;
 126
 127 (defun nrm2(n x incx)
 128   (let ((x (bv n x incx)))
 129     (sqrt (abs (sum (* (conjugate x) x))))))
 130
 131 (let ((x (coerce dx '(vector c-double)))
 132       (n 5)
 133       (maxdiff 0))
 134   (dolist (incx '(1 2 -1 -2))
 135     (let ((val (blas-dnrm2 n x 0 incx))
 136           (v (nrm2 n x incx)))
 137       (setf maxdiff (max maxdiff (abs (- val v))))))
 138   (test-report "DNRM2" maxdiff))
 139
 140 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 141       (n 5)
 142       (maxdiff 0))
 143   (dolist (incx '(1 2 -1 -2))
 144     (let ((val (blas-dznrm2 n x 0 incx))
 145           (v (nrm2 n x incx)))
 146       (setf maxdiff (max maxdiff (abs (- val v))))))
 147   (test-report "DZNRM2" maxdiff))
 148
 149
 150 ;;;;
 151 ;;;; _ROT
 152 ;;;;
 153
 154 (defun rot (n x incx y incy c s)
 155   (let ((x (bv n x incx))
 156         (y (bv n y incy)))
 157     (values (+ (* c x) (* s y))
 158             (- (* c y) (* s x)))))
 159
 160 (let* ((x (coerce dx '(vector c-double)))
 161        (y (coerce dy '(vector c-double)))
 162        (c .2)
 163        (s (sqrt (- 1 (* c c))))
 164        (n 5)
 165        (valx nil)
 166        (valy nil)
 167        (maxdiff 0))
 168   (dolist (incx '(1 2 -1 -2))
 169     (dolist (incy '(1 2 -1 -2))
 170       (setf valx (copy-array x))
 171       (setf valy (copy-array y))
 172       (blas-drot n valx 0 incx valy 0 incy c s)
 173       (multiple-value-bind
 174        (vx vy)
 175        (rot n x incx y incy c s)
 176        (let ((vvx (bv n valx incx))
 177              (vvy (bv n valy incy)))
 178          (setf maxdiff (max maxdiff (abs (- vvx vx)) (abs (- vvy vy))))))))
 179   (test-report "DROT" maxdiff))
 180
 181 (let* ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 182        (y (coerce (complex dy (* 3 dy)) '(vector c-dcomplex)))
 183        (c .2)
 184        (s (sqrt (- 1 (* c c))))
 185        (n 5)
 186        (valx nil)
 187        (valy nil)
 188        (maxdiff 0))
 189   (dolist (incx '(1 2 -1 -2))
 190     (dolist (incy '(1 2 -1 -2))
 191       (setf valx (copy-array x))
 192       (setf valy (copy-array y))
 193       (blas-zdrot n valx 0 incx valy 0 incy c s)
 194       (multiple-value-bind
 195        (vx vy)
 196        (rot n x incx y incy c s)
 197        (let ((vvx (bv n valx incx))
 198              (vvy (bv n valy incy)))
 199          (setf maxdiff (max maxdiff (abs (- vvx vx)) (abs (- vvy vy))))))))
 200   (test-report "ZDROT" maxdiff))
 201
 202
 203 ;;;;
 204 ;;;; _ROTG
 205 ;;;;
 206
 207 (defun drotg (da db)
 208   (let ((roe (if (> (abs da) (abs db)) da db))
 209         (scale (+ (abs da) (abs db))))
 210     (if (= scale 0.0)
 211         (values 0.0 0.0 1.0 0.0)
 212         (let* ((r (* scale (sqrt (+ (^ (/ da scale) 2) (^ (/ db scale) 2)))))
 213                (r (* (signum roe) r))
 214                (c (/ da r))
 215                (s (/ db r))
 216                (z 1.0))
 217           (if (> (abs da) (abs db)) (seft z s))
 218           (if (and (>= (abs db) (abs da)) (/= c 0.0)) (setf z (/ 1.0 c)))
 219           (values r z c s)))))
 220
 221 (defun zrotg (ca cb)
 222   (if (= (abs ca) 0.0)
 223       (values cb cb 0.0 1.0)
 224       (let* ((scale (+ (abs ca) (abs cb)))
 225              (norm (* scale (sqrt (+ (^ (abs (/ ca scale)) 2)
 226                                      (^ (abs (/ cb scale)) 2)))))
 227              (alpha (/ ca (abs ca)))
 228              (c (/ (abs ca) norm))
 229              (s (/ (* alpha (conjugate cb)) norm))
 230              (ca (* alpha norm)))
 231         (values ca cb c s))))
 232
 233 (let ((a 1)
 234       (b 7))
 235   (test-report "DROTG"
 236                (max (abs (- (multiple-value-list (blas-drotg a b))
 237                             (multiple-value-list (drotg a b)))))))
 238
 239 (let ((a #c(1 2))
 240       (b #c(7 3)))
 241   (test-report "ZROTG"
 242                (max (abs (- (multiple-value-list (blas-zrotg a b))
 243                             (multiple-value-list (zrotg a b)))))))
 244
 245 ;;;;
 246 ;;;; _SCAL
 247 ;;;;
 248
 249 (defun scal (n a x incx)
 250   (let ((x (bv n x incx)))
 251     (* a x)))
 252
 253 (let ((x (coerce dx '(vector c-double)))
 254       (a 3.0)
 255       (n 5)
 256       (val nil)
 257       (maxdiff 0))
 258   (dolist (incx '(1 2 -1 -2))
 259     (setf val (copy-array x))
 260     (blas-dscal n a val 0 incx)
 261     (let ((v (scal n a x incx))
 262           (vv (bv n val incx)))
 263       (setf maxdiff (max maxdiff (max (abs (- vv v)))))))
 264   (test-report "DSCAL" maxdiff))
 265
 266 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 267       (a 3.0)
 268       (n 5)
 269       (val nil)
 270       (maxdiff 0))
 271   (dolist (incx '(1 2 -1 -2))
 272     (setf val (copy-array x))
 273     (blas-zdscal n a val 0 incx)
 274     (let ((v (scal n a x incx))
 275           (vv (bv n val incx)))
 276       (setf maxdiff (max maxdiff (max (abs (- vv v)))))))
 277   (test-report "ZDSCAL" maxdiff))
 278
 279 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 280       (a (complex 3.0 6.0))
 281       (n 5)
 282       (val nil)
 283       (maxdiff 0))
 284   (dolist (incx '(1 2 -1 -2))
 285     (setf val (copy-array x))
 286     (blas-zscal n a val 0 incx)
 287     (let ((v (scal n a x incx))
 288           (vv (bv n val incx)))
 289       (setf maxdiff (max maxdiff (max (abs (- vv v)))))))
 290   (test-report "ZSCAL" maxdiff))
 291
 292
 293 ;;;;
 294 ;;;; _SWAP
 295 ;;;;
 296
 297 (let ((x (coerce dx '(vector c-double)))
 298       (y (coerce dy '(vector c-double)))
 299       (n 5)
 300       (valx nil)
 301       (valy nil)
 302       (maxdiff 0))
 303   (dolist (incx '(1 2 -1 -2))
 304     (dolist (incy '(1 2 -1 -2))
 305       (setf valx (copy-array x))
 306       (setf valy (copy-array y))
 307       (blas-dswap n valx 0 incx valy 0 incy)
 308       (let ((vx (bv n x incx))
 309             (vy (bv n y incy))
 310             (vvx (bv n valx incx))
 311             (vvy (bv n valy incy)))
 312         (setf maxdiff
 313               (max maxdiff (max (abs (- vvx vy)) (abs (- vvy vx))))))))
 314   (test-report "DSWAP" maxdiff))
 315
 316 (let ((x (coerce dx '(vector c-dcomplex)))
 317       (y (coerce dy '(vector c-dcomplex)))
 318       (n 5)
 319       (valx nil)
 320       (valy nil)
 321       (maxdiff 0))
 322   (dolist (incx '(1 2 -1 -2))
 323     (dolist (incy '(1 2 -1 -2))
 324       (setf valx (copy-array x))
 325       (setf valy (copy-array y))
 326       (blas-zswap n valx 0 incx valy 0 incy)
 327       (let ((vx (bv n x incx))
 328             (vy (bv n y incy))
 329             (vvx (bv n valx incx))
 330             (vvy (bv n valy incy)))
 331         (setf maxdiff
 332               (max maxdiff (max (abs (- vvx vy)) (abs (- vvy vx))))))))
 333   (test-report "ZSWAP" maxdiff))
 334
 335 ;;;;
 336 ;;;; I_AMAX
 337 ;;;;
 338
 339 (defun iamax (n x incx)
 340   (let* ((x (bv n x incx))
 341          (ax (+ (abs (realpart x)) (abs (imagpart x)))))
 342     (position (max ax) ax)))
 343
 344 (let ((x (coerce dx '(vector c-double)))
 345       (n 5)
 346       (maxdiff 0))
 347   (dolist (incx '(1 2 -1 -2))
 348     (let ((val (blas-idamax n x 0 incx))
 349           (v (iamax n x incx)))
 350       (setf maxdiff (max maxdiff (abs (- val v))))))
 351   (test-report "IDAMAX" maxdiff))
 352
 353 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 354       (n 5)
 355       (maxdiff 0))
 356   (dolist (incx '(1 2 -1 -2))
 357     (let ((val (blas-izamax n x 0 incx))
 358           (v (iamax n x incx)))
 359       (setf maxdiff (max maxdiff (abs (- val v))))))
 360   (test-report "IZAMAX" maxdiff))
 361
 362
 363 ;;;;
 364 ;;;; _COPY
 365 ;;;;
 366
 367 (let ((x (coerce dx '(vector c-double)))
 368       (y (coerce dy '(vector c-double)))
 369       (n 5)
 370       (val nil)
 371       (maxdiff 0))
 372   (dolist (incx '(1 2 -1 -2))
 373     (dolist (incy '(1 2 -1 -2))
 374       (setf val (copy-array y))
 375       (blas-dcopy n x 0 incx val 0 incy)
 376       (let ((v (bv n x incx))
 377             (vv (bv n val incy)))
 378         (setf maxdiff (max maxdiff (max (abs (- v vv))))))))
 379   (test-report "DCOPY" maxdiff))
 380
 381 (let ((x (coerce (complex dx (* 2 dx)) '(vector c-dcomplex)))
 382       (y (coerce (complex dy (* 2 dy)) '(vector c-dcomplex)))
 383       (n 5)
 384       (val nil)
 385       (maxdiff 0))
 386   (dolist (incx '(1 2 -1 -2))
 387     (dolist (incy '(1 2 -1 -2))
 388       (setf val (copy-array y))
 389       (blas-zcopy n x 0 incx val 0 incy)
 390       (let ((v (bv n x incx))
 391             (vv (bv n val incy)))
 392         (setf maxdiff (max maxdiff (max (abs (- v vv))))))))
 393   (test-report "ZCOPY" maxdiff))
 394
 395 (setf da '(1 2 3 4 5 6))
 396 (setf db '(1 2 3 4 5 6 13 14 15))
 397 (setf dc '(16 17 18 19 20 21))
 398 (setf dx '(7 8 9))
 399 (setf dy '(10 11 12))
 400
 401 (defvar *verbose-tests* nil)
 402
 403 ;(defun test-report (name val) (format t "~&~a: ~8t~f~%" name val))
 404
 405
 406 ;;;;
 407 ;;;; _GEMV
 408 ;;;;
 409
 410 (defun gemv (trans alpha a x incx beta y incy)
 411   (when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
 412   (when (/= (abs incy) 1) (error "increment has to be +1 or -1"))
 413   (when (equalp trans "n") (setf a (transpose a)))
 414   (when (equalp trans "c") (setf a (conjugate a)))
 415   (let* ((m (array-dimension a 0))
 416          (n (array-dimension a 1))
 417          (z (make-array m :initial-element 0)))
 418     (when (minusp incx) (setf x (reverse (subseq x 0 n))))
 419     (when (minusp incy) (setf y (reverse (subseq y 0 m))))
 420     (replace z (* beta y))
 421     (dotimes (i m)
 422       (dotimes (j n)
 423         (incf (aref z i) (* alpha (aref a i j) (aref x j)))))
 424     (if (minusp incy) (reverse z) z)))
 425
 426 (let ((a (coerce da '(vector c-double)))
 427       (x (coerce dx '(vector c-double)))
 428       (y (coerce dy '(vector c-double)))
 429       (val nil)
 430       (alpha 2)
 431       (beta 3)
 432       (maxdiff 0))
 433   (dolist (trans '("n" "t"))
 434     (dolist (incx '(1 -1))
 435       (dolist (incy '(1 -1))
 436         (setf val (copy-seq y))
 437         (blas-dgemv trans 3 2 alpha a 0 3 x 0 incx beta val 0 incy)
 438         (let ((v (gemv trans alpha (matrix '(2 3) a) x incx beta y incy))
 439               (vv (if (equalp trans "n") val (subseq val 0 2))))
 440           (when *verbose-tests*
 441                 (format t "DGEMV ~a ~2d ~2d: ~f~%"
 442                         trans
 443                         incx
 444                         incy
 445                         (max (abs (- vv v)))))
 446           (setf maxdiff (max maxdiff (max (abs (- vv v)))))))))
 447   (test-report "DGEMV" maxdiff))
 448
 449 (let ((a (coerce (complex da (* 2 da)) '(vector c-dcomplex)))
 450       (x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
 451       (y (coerce (complex dy (* 4 dy)) '(vector c-dcomplex)))
 452       (val nil)
 453       (alpha (complex 2 4))
 454       (beta (complex 3 6))
 455       (maxdiff 0))
 456   (dolist (trans '("n" "t" "c"))
 457     (dolist (incx '(1 -1))
 458       (dolist (incy '(1 -1))
 459         (setf val (copy-seq y))
 460         (blas-zgemv trans 3 2 alpha a 0 3 x 0 incx beta val 0 incy)
 461         (let ((v (gemv trans alpha (matrix '(2 3) a) x incx beta y incy))
 462               (vv (if (equalp trans "n") val (subseq val 0 2))))
 463           (when *verbose-tests*
 464                 (format t "ZGEMV ~a ~2d ~2d: ~f~%"
 465                         trans
 466                         incx
 467                         incy
 468                         (max (abs (- vv v)))))
 469           (setf maxdiff (max maxdiff (max (abs (- vv v)))))))))
 470   (test-report "ZGEMV" maxdiff))
 471
 472 ;;;;
 473 ;;;; _GER_
 474 ;;;;
 475
 476 (defun ger (alpha x incx y incy a &optional conj)
 477   (when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
 478   (when (/= (abs incy) 1) (error "increment has to be +1 or -1"))
 479   (let ((x (subseq x 0 (array-dimension a 1)))
 480         (y (subseq y 0 (array-dimension a 0))))
 481     (when (minusp incx) (setf x (reverse x)))
 482     (when (minusp incy) (setf y (reverse y)))
 483     (+ (* alpha
 484           (outer-product (if conj (conjugate y) y) x))
 485        a)))
 486
 487 (let ((a (coerce da '(vector c-double)))
 488       (x (coerce dx '(vector c-double)))
 489       (y (coerce dy '(vector c-double)))
 490       (val nil)
 491       (alpha 2)
 492       (maxdiff 0))
 493   (dolist (incx '(1 -1))
 494     (dolist (incy '(1 -1))
 495       (setf val (copy-seq a))
 496       (blas-dger 3 2 alpha x 0 incx y 0 incy val 0 3)
 497       (let ((v (ger alpha x incx y incy (matrix '(2 3) a))))
 498         (when *verbose-tests*
 499               (format t "DGER ~2d ~2d: ~f~%"
 500                       incx
 501                       incy
 502                       (max (abs (- val v)))))
 503         (setf maxdiff (max maxdiff (max (abs (- val v))))))))
 504   (test-report "DGER" maxdiff))
 505
 506 (let ((a (coerce (complex da (* 2 da)) '(vector c-dcomplex)))
 507       (x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
 508       (y (coerce (complex dy (* 4 dy)) '(vector c-dcomplex)))
 509       (val nil)
 510       (alpha (complex 2 4))
 511       (maxdiff 0))
 512   (dolist (incx '(1 -1))
 513     (dolist (incy '(1 -1))
 514       (setf val (copy-seq a))
 515       (blas-zgerc 3 2 alpha x 0 incx y 0 incy val 0 3)
 516       (let ((v (ger alpha x incx y incy (matrix '(2 3) a) t)))
 517         (when *verbose-tests*
 518               (format t "ZGERC ~2d ~2d: ~f~%"
 519                       incx
 520                       incy
 521                       (max (abs (- val v)))))
 522         (setf maxdiff (max maxdiff (max (abs (- val v))))))))
 523   (test-report "ZGERC" maxdiff))
 524
 525 (let ((a (coerce (complex da (* 2 da)) '(vector c-dcomplex)))
 526       (x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
 527       (y (coerce (complex dy (* 4 dy)) '(vector c-dcomplex)))
 528       (val nil)
 529       (alpha (complex 2 4))
 530       (maxdiff 0))
 531   (dolist (incx '(1 -1))
 532     (dolist (incy '(1 -1))
 533       (setf val (copy-seq a))
 534       (blas-zgeru 3 2 alpha x 0 incx y 0 incy val 0 3)
 535       (let ((v (ger alpha x incx y incy (matrix '(2 3) a))))
 536         (when *verbose-tests*
 537               (format t "ZGERU ~2d ~2d: ~f~%"
 538                       incx
 539                       incy
 540                       (max (abs (- val v)))))
 541         (setf maxdiff (max maxdiff (max (abs (- val v))))))))
 542   (test-report "ZGERU" maxdiff))
 543
 544
 545 ;;;;
 546 ;;;; _TRMV
 547 ;;;;
 548
 549 (defun trmv (a x &optional upper transpose unit incx)
 550   (when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
 551   (cond
 552    ((equalp transpose "n")
 553     (setf a (transpose a)))
 554    (t (setf upper (if (equalp upper "u") "l" "u"))))
 555   (when (equalp transpose "c") (setf a (conjugate a)))
 556   (when (minusp incx) (setf x (reverse x)))
 557   (let* ((n (array-dimension a 0))
 558          (y (make-array n :initial-element 0)))
 559     (when (not (equalp unit "n"))
 560           (setf a (copy-array a))
 561           (dotimes (i n) (setf (aref a i i) 1)))
 562     (if (equalp upper "u")
 563         (dotimes (i n)
 564           (do ((j i (+ j 1)))
 565               ((>= j n))
 566               (incf (aref y i) (* (aref a i j) (aref x j)))))
 567         (dotimes (i n)
 568           (dotimes (j (+ i 1))
 569             (incf (aref y i) (* (aref a i j) (aref x j))))))
 570     (if (minusp incx) (reverse y) y)))
 571
 572 (let ((b (coerce db '(vector c-double)))
 573       (x (coerce dx '(vector c-double)))
 574       (val nil)
 575       (maxdiff 0))
 576   (dolist (upper '("u" "l"))
 577     (dolist (trans '("n" "t"))
 578       (dolist (diag '("n" "u"))
 579         (dolist (incx '(1 -1))
 580           (setf val (copy-seq x))
 581           (blas-dtrmv upper trans diag 3 b 0 3 val 0 incx)
 582           (let ((v (trmv (matrix '(3 3) b) x upper trans diag incx)))
 583             (when *verbose-tests*
 584                   (format t "DTRMV ~a ~a ~a ~2d: ~f~%"
 585                           upper
 586                           trans
 587                           diag
 588                           incx
 589                           (max (abs (- val v)))))
 590             (setf maxdiff (max maxdiff (max (abs (- val v))))))))))
 591   (test-report "DTRMV" maxdiff))
 592
 593 (let ((b (coerce (complex db (* 2 db)) '(vector c-dcomplex)))
 594       (x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
 595       (val nil)
 596       (maxdiff 0))
 597   (dolist (upper '("u" "l"))
 598     (dolist (trans '("n" "t" "c"))
 599       (dolist (diag '("n" "u"))
 600         (dolist (incx '(1 -1))
 601           (setf val (copy-seq x))
 602           (blas-ztrmv upper trans diag 3 b 0 3 val 0 incx)
 603           (let ((v (trmv (matrix '(3 3) b) x upper trans diag incx)))
 604             (when *verbose-tests*
 605                   (format t "ZTRMV ~a ~a ~a ~2d: ~f~%"
 606                           upper
 607                           trans
 608                           diag
 609                           incx
 610                           (max (abs (- val v)))))
 611             (setf maxdiff (max maxdiff (max (abs (- val v))))))))))
 612   (test-report "ZTRMV" maxdiff))
 613
 614 (defun trsv (a x &optional upper transpose unit incx)
 615   (when (/= (abs incx) 1) (error "increment has to be +1 or -1"))
 616   (cond
 617    ((equalp transpose "n")
 618     (setf a (transpose a)))
 619    (t (setf upper (if (equalp upper "u") "l" "u"))))
 620   (when (equalp transpose "c") (setf a (conjugate a)))
 621   (when (minusp incx) (setf x (reverse x)))
 622   (let* ((n (array-dimension a 0))
 623          (y (make-array n :initial-contents x)))
 624     (if (equalp upper "u")
 625         (do ((i (- n 1) (- i 1)))
 626             ((< i 0))
 627             (do ((j (+ i 1) (+ j 1)))
 628                 ((>= j n))
 629                 (decf (aref y i) (* (aref a i j) (aref y j))))
 630             (when (equalp unit "n")
 631                   (setf (aref y i) (/ (aref y i) (aref a i i)))))
 632         (dotimes (i n)
 633           (dotimes (j i)
 634             (decf (aref y i) (* (aref a i j) (aref y j))))
 635           (when (equalp unit "n")
 636                 (setf (aref y i) (/ (aref y i) (aref a i i))))))
 637     (if (minusp incx) (reverse y) y)))
 638
 639 (let ((b (coerce db '(vector c-double)))
 640       (x (coerce dx '(vector c-double)))
 641       (val nil)
 642       (maxdiff 0))
 643   (dolist (upper '("u" "l"))
 644     (dolist (trans '("n" "t"))
 645       (dolist (diag '("n" "u"))
 646         (dolist (incx '(1 -1))
 647           (setf val (copy-seq x))
 648           (blas-dtrsv upper trans diag 3 b 0 3 val 0 incx)
 649           (let ((v (trsv (matrix '(3 3) b) x upper trans diag incx)))
 650             (when *verbose-tests*
 651                   (format t "DTRSV ~a ~a ~a ~2d: ~f~%"
 652                           upper
 653                           trans
 654                           diag
 655                           incx
 656                           (max (abs (- val v)))))
 657             (setf maxdiff (max maxdiff (max (abs (- val v))))))))))
 658   (test-report "DTRSV" maxdiff))
 659
 660 (let ((b (coerce (complex db (* 2 db)) '(vector c-dcomplex)))
 661       (x (coerce (complex dx (* 3 dx)) '(vector c-dcomplex)))
 662       (val nil)
 663       (maxdiff 0))
 664   (dolist (upper '("u" "l"))
 665     (dolist (trans '("n" "t"))
 666       (dolist (diag '("n" "u"))
 667         (dolist (incx '(1 -1))
 668           (setf val (copy-seq x))
 669           (blas-ztrsv upper trans diag 3 b 0 3 val 0 incx)
 670           (let ((v (trsv (matrix '(3 3) b) x upper trans diag incx)))
 671             (when *verbose-tests*
 672                   (format t "ZTRSV ~a ~a ~a ~2d: ~f~%"
 673                           upper
 674                           trans
 675                           diag
 676                           incx
 677                           (max (abs (- val v)))))
 678             (setf maxdiff (max maxdiff (max (abs (- val v))))))))))
 679   (test-report "ZTRSV" maxdiff))
 680
 681
 682
 683 ;;;;
 684 ;;;; _GEMM
 685 ;;;;
 686
 687 (defun gemm (transa transb alpha a b beta c)
 688   (when (equalp transa "n") (setf a (transpose a)))
 689   (when (equalp transa "c") (setf a (conjugate a)))
 690   (when (equalp transb "n") (setf b (transpose b)))
 691   (when (equalp transb "c") (setf b (conjugate b)))
 692   (let* ((m (array-dimension c 1))
 693          (n (array-dimension c 0))
 694          (k (array-dimension a 1))
 695          (z (* beta c)))
 696     (dotimes (i m)
 697       (dotimes (j n)
 698         (dotimes (l k)
 699           (incf (aref z j i) (* alpha (aref a i l) (aref b l j))))))
 700     z))
 701
 702
 703 (setf aa (transpose (matrix '(4 2) (iseq 1 8))))
 704 (setf bb (transpose (matrix '(2 3) (iseq 9 14))))
 705 (setf cc (transpose (matrix '(4 3) (iseq 15 26))))
 706
 707 (let ((a (coerce aa '(array c-double)))
 708       (b (coerce bb '(array c-double)))
 709       (c (coerce cc '(array c-double)))
 710       (m 4)
 711       (n 3)
 712       (k 2)
 713       (alpha 2)
 714       (beta 3)
 715       (val nil)
 716       (maxdiff 0))
 717   (dolist (transa '("n" "t"))
 718     (dolist (transb '("n" "t"))
 719       (let ((aa (if (equalp transa "n") a (transpose a)))
 720             (bb (if (equalp transb "n") b (transpose b)))
 721             (lda (if (equalp transa "n") m k))
 722             (ldb (if (equalp transb "n") k n)))
 723         (setf val (copy-array c))
 724         (blas-dgemm transa transb m n k alpha aa 0 lda bb 0 ldb beta val 0 m)
 725         (let ((v (gemm transa transb alpha aa bb beta c)))
 726           (when *verbose-tests*
 727                 (format t "DGEMM ~a ~a: ~f~%"
 728                         transa
 729                         transb
 730                         (max (abs (- val v)))))
 731           (setf maxdiff (max maxdiff (max (abs (- val v)))))))))
 732   (test-report "DGEMM" maxdiff))
 733
 734 (let ((a (coerce (complex aa (* 2 aa)) '(array c-dcomplex)))
 735       (b (coerce (complex bb (* 3 bb)) '(array c-dcomplex)))
 736       (c (coerce (complex cc (* 4 cc)) '(array c-dcomplex)))
 737       (m 4)
 738       (n 3)
 739       (k 2)
 740       (alpha #c(2 4))
 741       (beta #c(3 6))
 742       (val nil)
 743       (maxdiff 0))
 744   (dolist (transa '("n" "t" "c"))
 745     (dolist (transb '("n" "t" "c"))
 746       (let ((aa (if (equalp transa "n") a (transpose a)))
 747             (bb (if (equalp transb "n") b (transpose b)))
 748             (lda (if (equalp transa "n") m k))
 749             (ldb (if (equalp transb "n") k n)))
 750         (setf val (copy-array c))
 751         (blas-zgemm transa transb m n k alpha aa 0 lda bb 0 ldb beta val 0 m)
 752         (let ((v (gemm transa transb alpha aa bb beta c)))
 753           (when *verbose-tests*
 754                 (format t "ZGEMM ~a ~a: ~f~%"
 755                         transa
 756                         transb
 757                         (max (abs (- val v)))))
 758           (setf maxdiff (max maxdiff (max (abs (- val v)))))))))
 759   (test-report "ZGEMM" maxdiff))
 760
 761
 762 ;;;;
 763 ;;;; TRSM
 764 ;;;;
 765
 766 (defun trsm (side upper transpose unit a alpha x)
 767   (cond
 768    ((equalp side "l")
 769     (let ((cx (* alpha (row-list x))))
 770       (apply #'bind-rows
 771              (mapcar #'(lambda (c) (trsv a c upper transpose unit 1))
 772                      cx))))
 773    (t
 774     (let ((cx (* alpha (column-list x)))
 775           (conj (equalp transpose "c")))
 776       (setf transpose (if (equalp transpose "n") "t" "n"))
 777       (when conj (setf cx (conjugate cx)))
 778       (apply #'bind-columns
 779              (mapcar #'(lambda (c)
 780                          (let ((v (trsv a c upper transpose unit 1)))
 781                            (if conj (conjugate v) v)))
 782                      cx))))))
 783
 784 (let ((b (coerce db '(vector c-double)))
 785       (c (coerce dc '(vector c-double)))
 786       (alpha 2)
 787       (val nil)
 788       (maxdiff 0.0))
 789   (dolist (side '("l" "r"))
 790     (dolist (uplo '("u" "l"))
 791       (dolist (trans '("n" "t"))
 792         (dolist (diag '("n" "u"))
 793           (setf val (copy-array c))
 794           (if (equalp side "l")
 795               (blas-dtrsm side uplo trans diag 3 2 alpha b 0 3 val 0 3)
 796               (blas-dtrsm side uplo trans diag 2 3 alpha b 0 3 val 0 2))
 797           (let* ((mb (matrix '(3 3) b))
 798                  (mc (matrix (if (equalp side "l") '(2 3) '(3 2)) c))
 799                  (v (trsm side uplo trans diag mb alpha mc)))
 800             (when *verbose-tests*
 801                   (format t "DTRSM ~a ~a ~a ~a: ~f~%"
 802                           side
 803                           uplo
 804                           trans
 805                           diag
 806                           (max (abs (- val v)))))
 807             (setf maxdiff (max maxdiff (max (abs (- val v))))))))))
 808   (test-report "DTRSM" maxdiff))
 809
 810 (let ((b (coerce (complex db (* 2 db)) '(vector c-dcomplex)))
 811       (c (coerce (complex dc (* 3 dc)) '(vector c-dcomplex)))
 812       (alpha (complex 2 4))
 813       (val nil)
 814       (maxdiff 0))
 815   (dolist (side '("l" "r"))
 816     (dolist (uplo '("u" "l"))
 817       (dolist (trans '("n" "t" "c"))
 818         (dolist (diag '("n" "u"))
 819           (setf val (copy-array c))
 820           (if (equalp side "l")
 821               (blas-ztrsm side uplo trans diag 3 2 alpha b 0 3 val 0 3)
 822               (blas-ztrsm side uplo trans diag 2 3 alpha b 0 3 val 0 2))
 823           (let* ((mb (matrix '(3 3) b))
 824                  (mc (matrix (if (equalp side "l") '(2 3) '(3 2)) c))
 825                  (v (trsm side uplo trans diag mb alpha mc)))
 826             (when *verbose-tests*
 827                   (format t "ZTRSM ~a ~a ~a ~a: ~f~%"
 828                           side
 829                           uplo
 830                           trans
 831                           diag
 832                           (max (abs (- val v)))))
 833             (setf maxdiff (max maxdiff (max (abs (- val v))))))))))
 834   (test-report "ZTRSM" maxdiff))
 835
 836
 837 (defun nax+y (a x y)
 838   (let* ((m (array-dimension a 0))
 839          (n (array-dimension a 1))
 840          (da (coerce a '(array c-double)))
 841          (dx (coerce x '(vector c-double)))
 842          (dy (make-array m :element-type 'c-double :initial-contents y)))
 843     (blas-dgemv "t" n m 1 da 0 n dx 0 1 1 dy 0 1)
 844     (coerce dy (if (vectorp x) '(vector t) 'list))))