Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-164-g938bc52
|
From: Akshay S. <ak...@us...> - 2012-08-04 15:15:33
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, tensor has been updated
via 938bc521fe3431d9a4cbcc0c7bab9c4bb616aaf4 (commit)
via 9b50685d6952b3be9ff29473595b2694ea234b08 (commit)
via 0b4fcdfe7d12f45c1d46f3b42589a5f2ff54e8dc (commit)
from 7424dca3c956b58d494e938ed7acf90abc79d086 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 938bc521fe3431d9a4cbcc0c7bab9c4bb616aaf4
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Aug 4 20:38:14 2012 +0530
o Changed the order of application in the ediv functions.
o Made most of the generated functions inline.
diff --git a/lib/matlisp/dediv.f b/lib/matlisp/dediv.f
index 8151b5f..6b51547 100644
--- a/lib/matlisp/dediv.f
+++ b/lib/matlisp/dediv.f
@@ -26,7 +26,7 @@
* code for both increments equal to 1
*
20 do 30 i = 1,n
- dy(i) = dy(i) / dx(i)
+ dy(i) = dx(i) / dy(i)
30 continue
diff --git a/lib/matlisp/zediv.f b/lib/matlisp/zediv.f
index c1308d9..b0e8b21 100644
--- a/lib/matlisp/zediv.f
+++ b/lib/matlisp/zediv.f
@@ -26,7 +26,7 @@
* code for both increments equal to 1
*
20 do 30 i = 1,n
- dy(i) = dy(i) / dx(i)
+ dy(i) = dx(i) / dy(i)
30 continue
diff --git a/src/foreign-core/libmatlisp.lisp b/src/foreign-core/libmatlisp.lisp
index 8e15a2f..4aceebf 100644
--- a/src/foreign-core/libmatlisp.lisp
+++ b/src/foreign-core/libmatlisp.lisp
@@ -5,7 +5,7 @@
(DESCAL n dx incx dy incy)
Multiplies the vector X and Y element-wise.
- Y <- Y .* E
+ Y <- Y .* X
"
(n :integer :input)
(dx (* :double-float :inc head-x) :input)
@@ -18,7 +18,7 @@
(ZESCAL n dx incx dy incy)
Multiplies the vector X and Y element-wise.
- Y <- Y .* E
+ Y <- Y .* X
"
(n :integer :input)
(dx (* :complex-double-float :inc head-x) :input)
@@ -31,7 +31,7 @@
(DEDIV n dx incx dy incy)
Divides the vector Y by X element-wise.
- Y <- Y .* E
+ Y <- X ./ Y
"
(n :integer :input)
(dx (* :double-float :inc head-x) :input)
@@ -44,7 +44,7 @@
(ZEDIV n dx incx dy incy)
Divide the vector Y by X element-wise.
- Y <- Y .* E
+ Y <- X ./ Y
"
(n :integer :input)
(dx (* :complex-double-float :inc head-x) :input)
diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp
index 3e13987..ec77357 100644
--- a/src/level-1/axpy.lisp
+++ b/src/level-1/axpy.lisp
@@ -34,7 +34,7 @@
;;Use only after checking the arguments for compatibility.
(let* ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (alpha from to)
+ `(definline ,func (alpha from to)
(declare (type ,tensor-class from to)
(type ,(getf opt :element-type) alpha))
(let ((strd-p (blas-copyable-p from to))
@@ -73,7 +73,7 @@
;;Use only after checking the arguments for compatibility.
(let* ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (num-from to)
+ `(definline ,func (num-from to)
(declare (type ,tensor-class to)
(type ,(getf opt :element-type) num-from))
(let ((min-strd (consecutive-store-p to))
@@ -82,7 +82,8 @@
((and min-strd call-fortran?)
(let ((num-array (,(getf opt :store-allocator) 1)))
(declare (type ,(linear-array-type (getf opt :store-type)) num-array))
- ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) 1) 'num-array 0)
+ (let-typed ((id (,(getf opt :coercer) 1) :type ,(getf opt :element-type)))
+ ,(funcall (getf opt :value-writer) `id 'num-array 0))
(,blas-func (number-of-elements to) num-from
num-array 0
(store to) min-strd
@@ -148,6 +149,8 @@
(real-typed-axpy! (coerce-real alpha) x y))
(defmethod axpy! ((alpha number) (x real-tensor) (y complex-tensor))
+ ;;Weird, shouldn't SBCL know this already ?
+ (declare (type complex-tensor y))
(let ((tmp (tensor-realpart~ y)))
(declare (type real-tensor tmp))
(etypecase alpha
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index a981a66..a9ee5ee 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -34,7 +34,7 @@
;;Use only after checking the arguments for compatibility.
(let* ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (from to)
+ `(definline ,func (from to)
(declare (type ,tensor-class from to))
(let ((strd-p (blas-copyable-p from to))
(call-fortran? (> (number-of-elements to) ,fortran-lb)))
@@ -66,7 +66,7 @@
;;Use only after checking the arguments for compatibility.
(let* ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (num-from to)
+ `(definline ,func (num-from to)
(declare (type ,tensor-class to)
(type ,(getf opt :element-type) num-from))
(let ((min-stride (consecutive-store-p to))
diff --git a/src/level-1/dot.lisp b/src/level-1/dot.lisp
index 8b777ee..dfb7d42 100644
--- a/src/level-1/dot.lisp
+++ b/src/level-1/dot.lisp
@@ -27,7 +27,7 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package #:matlisp)
-(defun real-typed-dot (x y conjugate-p)
+(definline real-typed-dot (x y conjugate-p)
(declare (type real-vector x y)
(ignore conjugate-p))
(let ((call-fortran? (> (number-of-elements x)
@@ -52,7 +52,7 @@
summing (* (aref sto-x of-x) (aref sto-y of-y)) into dot of-type real-type
finally (return dot))))))))
-(defun complex-typed-dot (x y conjugate-p)
+(definline complex-typed-dot (x y conjugate-p)
(declare (type complex-vector x y))
(let ((call-fortran? (> (number-of-elements x)
*complex-l1-fcall-lb*)))
@@ -144,7 +144,8 @@
(real-typed-dot x y nil))
(defmethod dot ((x real-vector) (y complex-vector) &optional (conjugate-p t))
- (declare (ignore conjugate-p))
+ (declare (ignore conjugate-p)
+ (type complex-vector y))
(let ((vw.y (tensor-realpart~ y)))
(declare (type real-vector vw.y))
(let ((rpart (prog1 (real-typed-dot x vw.y nil)
diff --git a/src/level-1/realimag.lisp b/src/level-1/realimag.lisp
index ec51f8a..adbe4e2 100644
--- a/src/level-1/realimag.lisp
+++ b/src/level-1/realimag.lisp
@@ -28,7 +28,7 @@
(in-package #:matlisp)
-(defun tensor-realpart~ (tensor)
+(definline tensor-realpart~ (tensor)
"
Syntax
======
@@ -50,7 +50,7 @@
:head (the index-type (* 2 (head tensor)))))
(number (realpart tensor))))
-(defun tensor-imagpart~ (tensor)
+(definline tensor-imagpart~ (tensor)
"
Syntax
======
diff --git a/src/level-1/scal.lisp b/src/level-1/scal.lisp
index 2687144..49e1cbc 100644
--- a/src/level-1/scal.lisp
+++ b/src/level-1/scal.lisp
@@ -31,7 +31,7 @@
(defmacro generate-typed-scal! (func (tensor-class fortran-func fortran-lb))
(let* ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (from to)
+ `(definline ,func (from to)
(declare (type ,tensor-class from to))
(let ((strd-p (blas-copyable-p from to))
(call-fortran? (> (number-of-elements to) ,fortran-lb)))
@@ -56,10 +56,32 @@
,(funcall (getf opt :value-writer) 'mul 't-sto 't-of))))))))
to)))
+(defmacro generate-typed-num-scal! (func (tensor-class blas-func fortran-lb))
+ (let ((opt (get-tensor-class-optimization tensor-class)))
+ (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
+ `(definline ,func (alpha to)
+ (declare (type ,tensor-class to)
+ (type ,(getf opt :element-type) alpha))
+ (let ((min-stride (consecutive-store-p to))
+ (call-fortran? (> (number-of-elements to) ,fortran-lb)))
+ (cond
+ ((and min-stride call-fortran?)
+ (,blas-func (number-of-elements to) alpha (store to) min-stride (head to)))
+ (t
+ (let ((t-sto (store to)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) t-sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions to))
+ with (linear-sums
+ (t-of (strides to) (head to)))
+ do (let ((scal-val (* ,(funcall (getf opt :reader) 't-sto 't-of) alpha)))
+ ,(funcall (getf opt :value-writer) 'scal-val 't-sto 't-of))))))))
+ to)))
+
(defmacro generate-typed-div! (func (tensor-class fortran-func fortran-lb))
(let* ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (from to)
+ `(definline ,func (from to)
(declare (type ,tensor-class from to))
(let ((strd-p (blas-copyable-p from to))
(call-fortran? (> (number-of-elements to) ,fortran-lb)))
@@ -80,21 +102,25 @@
(t-of (strides to) (head to)))
do (let*-typed ((val-f ,(funcall (getf opt :reader) 'f-sto 'f-of) :type ,(getf opt :element-type))
(val-t ,(funcall (getf opt :reader) 't-sto 't-of) :type ,(getf opt :element-type))
- (mul (/ val-t val-f) :type ,(getf opt :element-type)))
+ (mul (/ val-f val-t) :type ,(getf opt :element-type)))
,(funcall (getf opt :value-writer) 'mul 't-sto 't-of))))))))
to)))
-(defmacro generate-typed-num-scal! (func (tensor-class blas-func fortran-lb))
+(defmacro generate-typed-num-div! (func (tensor-class fortran-func fortran-lb))
(let ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (alpha to)
+ `(definline ,func (alpha to)
(declare (type ,tensor-class to)
(type ,(getf opt :element-type) alpha))
(let ((min-stride (consecutive-store-p to))
(call-fortran? (> (number-of-elements to) ,fortran-lb)))
(cond
((and min-stride call-fortran?)
- (,blas-func (number-of-elements to) alpha (store to) min-stride (head to)))
+ (let ((num-array (,(getf opt :store-allocator) 1)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) num-array))
+ (let-typed ((id (,(getf opt :coercer) 1) :type ,(getf opt :element-type)))
+ ,(funcall (getf opt :value-writer) `id 'num-array 0))
+ (,fortran-func (number-of-elements to) num-array 0 (store to) min-stride (head to))))
(t
(let ((t-sto (store to)))
(declare (type ,(linear-array-type (getf opt :store-type)) t-sto))
@@ -102,7 +128,7 @@
(mod-dotimes (idx (dimensions to))
with (linear-sums
(t-of (strides to) (head to)))
- do (let ((scal-val (* ,(funcall (getf opt :reader) 't-sto 't-of) alpha)))
+ do (let-typed ((scal-val (/ alpha ,(funcall (getf opt :reader) 't-sto 't-of)) :type ,(getf opt :element-type)))
,(funcall (getf opt :value-writer) 'scal-val 't-sto 't-of))))))))
to)))
@@ -116,6 +142,9 @@
(generate-typed-div! real-typed-div!
(real-tensor dediv *real-l1-fcall-lb*))
+(generate-typed-num-div! real-typed-num-div!
+ (real-tensor dediv *real-l1-fcall-lb*))
+
;;Complex
(definline zordscal (nele alpha x incx &optional hd-x)
(if (zerop (imagpart alpha))
@@ -128,7 +157,10 @@
(generate-typed-scal! complex-typed-scal!
(complex-tensor zescal *complex-l1-fcall-lb*))
-(generate-typed-scal! complex-typed-div!
+(generate-typed-div! complex-typed-div!
+ (complex-tensor zediv *complex-l1-fcall-lb*))
+
+(generate-typed-num-div! complex-typed-num-div!
(complex-tensor zediv *complex-l1-fcall-lb*))
;;---------------------------------------------------------------;;
@@ -175,32 +207,28 @@
Purpose
=======
- X <- X ./ alpha
-
- Yes the calling order is twisted.
+ X <- alpha ./ X
")
(:method :before ((x standard-tensor) (y standard-tensor))
(assert (lvec-eq (dimensions x) (dimensions y) #'=) nil
'tensor-dimension-mismatch)))
(defmethod div! ((alpha number) (x real-tensor))
- (real-typed-num-scal! (coerce-real (/ 1 alpha)) x))
+ (real-typed-num-div! (coerce-real alpha) x))
(defmethod div! ((x real-tensor) (y real-tensor))
(real-typed-div! x y))
(defmethod div! ((alpha number) (x complex-tensor))
- (complex-typed-num-scal! (coerce-complex (/ 1 alpha)) x))
+ (complex-typed-num-div! (coerce-complex alpha) x))
(defmethod div! ((x complex-tensor) (y complex-tensor))
(complex-typed-div! x y))
(defmethod div! ((x real-tensor) (y complex-tensor))
- (let ((tmp (tensor-realpart~ y)))
- (real-typed-div! x tmp)
- ;;Move view to the imaginary part
- (incf (head tmp))
- (real-typed-div! x tmp)))
+ ;;The alternative is worse!
+ (let ((tmp (copy! x (apply #'make-complex-tensor (lvec->list (dimensions x))))))
+ (complex-typed-div! tmp y)))
;;
(defgeneric scal (alpha x)
@@ -223,6 +251,9 @@
(defmethod scal ((alpha number) (x number))
(* alpha x))
+(defmethod scal ((x standard-tensor) (alpha number))
+ (scal alpha x))
+
(defmethod scal ((alpha number) (x real-tensor))
(let ((result (if (complexp alpha)
(copy! x (apply #'make-complex-tensor (lvec->list (dimensions x))))
@@ -240,7 +271,11 @@
(let ((result (copy x)))
(scal! alpha result)))
-(defmethod scal ((x standard-tensor) (y complex-tensor))
+(defmethod scal ((x real-tensor) (y complex-tensor))
+ (let ((result (copy y)))
+ (scal! x result)))
+
+(defmethod scal ((x complex-tensor) (y complex-tensor))
(let ((result (copy y)))
(scal! x result)))
@@ -253,7 +288,7 @@
Purpose
=======
- X <- X ./ alpha
+ alpha ./ X
Yes the calling order is twisted.
"))
@@ -261,15 +296,23 @@
(defmethod div ((alpha number) (x number))
(/ x alpha))
+(defmethod div ((x standard-tensor) (y number))
+ (let ((result (copy x)))
+ (scal! (/ 1 y) result)))
+
+(defmethod div ((x (eql nil)) (y standard-tensor))
+ (let ((result (copy y)))
+ (div! 1 result)))
+
+(defmethod div ((x real-tensor) (y real-tensor))
+ (div! x (copy y)))
+
(defmethod div ((alpha number) (x real-tensor))
(let ((result (if (complexp alpha)
(copy! x (apply #'make-complex-tensor (lvec->list (dimensions x))))
(copy x))))
(div! alpha result)))
-(defmethod div ((x real-tensor) (y real-tensor))
- (div! x (copy y)))
-
(defmethod div ((x complex-tensor) (y real-tensor))
(let ((result (copy! y (apply #'make-complex-tensor (lvec->list (dimensions x))))))
(div! x result)))
@@ -278,14 +321,11 @@
(let ((result (copy x)))
(div! alpha result)))
-(defmethod div ((x standard-tensor) (y complex-tensor))
+(defmethod div ((x real-tensor) (y complex-tensor))
(let ((result (copy y)))
(div! x result)))
-(defmethod div ((x real-tensor) (y (eql nil)))
- (let ((result (copy! 1 (apply #'make-real-tensor (lvec->list (dimensions x))))))
+(defmethod div ((x complex-tensor) (y complex-tensor))
+ (let ((result (copy y)))
(div! x result)))
-(defmethod div ((x complex-tensor) (y (eql nil)))
- (let ((result (copy! 1 (apply #'make-complex-tensor (lvec->list (dimensions x))))))
- (div! x result)))
diff --git a/src/level-1/swap.lisp b/src/level-1/swap.lisp
index 081f8c3..8c4a790 100644
--- a/src/level-1/swap.lisp
+++ b/src/level-1/swap.lisp
@@ -34,7 +34,7 @@
;;Use only after checking the arguments for compatibility.
(let* ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
- `(defun ,func (x y)
+ `(definline ,func (x y)
(declare (type ,tensor-class x y))
(let ((strd-p (blas-copyable-p x y))
(call-fortran? (> (number-of-elements x) ,fortran-lb)))
diff --git a/src/level-1/tensor-maker.lisp b/src/level-1/tensor-maker.lisp
index 541d1d2..6f7cb3b 100644
--- a/src/level-1/tensor-maker.lisp
+++ b/src/level-1/tensor-maker.lisp
@@ -5,48 +5,50 @@
(cocl (get-tensor-counterclass tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
(assert cocl nil 'tensor-cannot-find-counter-class :tensor-class tensor-class)
- `(defun ,func-name (&rest args)
- (labels ((make-dims (dims)
- (declare (type cons dims))
- (let*-typed ((vdim (make-index-store dims) :type index-store-vector)
- (ss (very-quickly (lvec-foldl #'(lambda (x y) (the index-type (* x y))) vdim)))
- (store (,(getf opt :store-allocator) ss))
- (rnk (length vdim)))
- (make-instance (case rnk (2 ',(getf cocl :matrix)) (1 ',(getf cocl :vector)) (t ',tensor-class))
- :store store :dimensions vdim)))
- (make-from-array (arr)
- (declare (type (array * *) arr))
- (let* ((ret (make-dims (array-dimensions arr)))
- (st-r (store ret))
- (lst (make-list (rank ret))))
- (declare (type ,tensor-class ret)
- (type ,(linear-array-type (getf opt :store-type)) st-r))
- (mod-dotimes (idx (dimensions ret))
- with (linear-sums
- (of-r (strides ret) (head ret)))
- do ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) (apply #'aref arr (lvec->list! idx lst))) 'st-r 'of-r))
- ret))
- (make-from-list (lst)
- (let* ((ret (make-dims (list-dimensions lst)))
- (st-r (store ret)))
- (declare (type ,tensor-class ret)
- (type ,(linear-array-type (getf opt :store-type)) st-r))
- (list-loop (idx ele lst)
- with (linear-sums
- (of-r (strides ret) (head ret)))
- do ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) ele) 'st-r 'of-r))
- ret)))
- (let ((largs (length args)))
- (if (= largs 1)
- (etypecase (first args)
- (array
- (make-from-array (first args)))
- (cons
- (make-from-list (first args)))
- (integer
- (make-dims (list (first args)))))
- (make-dims args)))))))
-
+ `(progn
+ (declaim (ftype (function (&rest t) ,tensor-class) ,func-name))
+ (defun ,func-name (&rest args)
+ (labels ((make-dims (dims)
+ (declare (type cons dims))
+ (let*-typed ((vdim (make-index-store dims) :type index-store-vector)
+ (ss (very-quickly (lvec-foldl #'(lambda (x y) (the index-type (* x y))) vdim)))
+ (store (,(getf opt :store-allocator) ss))
+ (rnk (length vdim)))
+ (make-instance (case rnk (2 ',(getf cocl :matrix)) (1 ',(getf cocl :vector)) (t ',tensor-class))
+ :store store :dimensions vdim)))
+ (make-from-array (arr)
+ (declare (type (array * *) arr))
+ (let* ((ret (make-dims (array-dimensions arr)))
+ (st-r (store ret))
+ (lst (make-list (rank ret))))
+ (declare (type ,tensor-class ret)
+ (type ,(linear-array-type (getf opt :store-type)) st-r))
+ (mod-dotimes (idx (dimensions ret))
+ with (linear-sums
+ (of-r (strides ret) (head ret)))
+ do ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) (apply #'aref arr (lvec->list! idx lst))) 'st-r 'of-r))
+ ret))
+ (make-from-list (lst)
+ (let* ((ret (make-dims (list-dimensions lst)))
+ (st-r (store ret)))
+ (declare (type ,tensor-class ret)
+ (type ,(linear-array-type (getf opt :store-type)) st-r))
+ (list-loop (idx ele lst)
+ with (linear-sums
+ (of-r (strides ret) (head ret)))
+ do ,(funcall (getf opt :value-writer) `(,(getf opt :coercer) ele) 'st-r 'of-r))
+ ret)))
+ (let ((largs (length args)))
+ (if (= largs 1)
+ (etypecase (first args)
+ (array
+ (make-from-array (first args)))
+ (cons
+ (make-from-list (first args)))
+ (integer
+ (make-dims (list (first args)))))
+ (make-dims args))))))))
+
(make-tensor-maker make-real-tensor (real-tensor))
(make-tensor-maker make-complex-tensor (complex-tensor))
diff --git a/src/level-1/trans.lisp b/src/level-1/trans.lisp
index d5c0087..b0de83a 100644
--- a/src/level-1/trans.lisp
+++ b/src/level-1/trans.lisp
@@ -27,7 +27,7 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package #:matlisp)
-(defun transpose! (A &optional permutation)
+(definline transpose! (A &optional permutation)
"
Syntax
======
@@ -58,10 +58,10 @@
(rotatef (aref strd-A (1- rnk)) (aref strd-A 0))))
A)
-(defun (setf transpose!) (value A &optional permutation)
+(definline (setf transpose!) (value A &optional permutation)
(copy! value (transpose! A permutation)))
-(defun transpose~ (A &optional permutation)
+(definline transpose~ (A &optional permutation)
"
Syntax
======
@@ -86,7 +86,7 @@
:parent-tensor A)))
(transpose! displaced permutation)))
-(defun (setf transpose~) (value A &optional permutation)
+(definline (setf transpose~) (value A &optional permutation)
(declare (type standard-tensor A))
(copy! value (transpose~ A permutation)))
@@ -109,7 +109,7 @@
(declare (type standard-tensor A))
(copy (transpose~ A permutation)))
-(defun (setf transpose) (value A &optional permutation)
+(definline (setf transpose) (value A &optional permutation)
(declare (type standard-tensor A))
(copy! value (transpose~ A permutation)))
@@ -125,7 +125,7 @@
;;
-(defun mconjugate! (A)
+(definline mconjugate! (A)
"
Syntax
======
@@ -154,7 +154,7 @@
Like mconjugate!, but non-destructive."
(etypecase A
(standard-tensor (mconjugate! (copy A)))
- (number (conjugate A))))
+ (number (cl:conjugate A))))
;;
(defun htranspose! (A &optional permutation)
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index 78004a5..8891d33 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -9,7 +9,7 @@
;;Use only after checking the arguments for compatibility.
(let* ((opt (get-tensor-class-optimization matrix-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class)
- `(defun ,func (alpha A x beta y job)
+ `(definline ,func (alpha A x beta y job)
(declare (type ,(getf opt :element-type) alpha beta)
(type ,matrix-class A)
(type ,vector-class x y)
@@ -70,7 +70,7 @@
(generate-typed-gemv! complex-base-typed-gemv!
(complex-matrix complex-vector zgemv *complex-l2-fcall-lb*))
-(defun complex-typed-gemv! (alpha A x beta y job)
+(definline complex-typed-gemv! (alpha A x beta y job)
(declare (type complex-matrix A)
(type complex-vector x y)
(type complex-type alpha beta)
@@ -78,10 +78,14 @@
(if (member job '(:n :t))
(complex-base-typed-gemv! alpha A x beta y job)
;;The CBLAS way.
- (let ((cx (mconjugate x)))
+ (let-typed ((cx (let-typed ((ret (apply #'make-real-tensor (lvec->list (dimensions x))) :type complex-vector))
+ (complex-typed-axpy! #c(-1d0 0d0) x ret))
+ :type complex-vector))
+ (complex-typed-num-scal! #c(-1d0 0d0) (tensor-realpart~ y))
(complex-base-typed-gemv! (cl:conjugate alpha) A cx
- (cl:conjugate beta) (mconjugate! y) (ecase job (:h :t) (:c :n)))
- (mconjugate! y))))
+ (cl:conjugate beta) y (ecase job (:h :t) (:c :n)))
+ (complex-typed-num-scal! #c(-1d0 0d0) (tensor-realpart~ y))
+ y)))
;;---------------------------------------------------------------;;
diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp
index b8124c3..386c193 100644
--- a/src/level-3/gemm.lisp
+++ b/src/level-3/gemm.lisp
@@ -31,7 +31,7 @@
(defmacro generate-typed-gemm! (func (matrix-class blas-gemm-func blas-gemv-func fortran-lb-parameter))
(let* ((opt (get-tensor-class-optimization matrix-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class)
- `(defun ,func (alpha A B beta C job)
+ `(definline ,func (alpha A B beta C job)
(declare (type ,(getf opt :element-type) alpha beta)
(type ,matrix-class A B C)
(type symbol job))
@@ -167,7 +167,7 @@
(generate-typed-gemm! complex-base-typed-gemm!
(complex-matrix zgemm zgemv *complex-l3-fcall-lb*))
-(defun complex-typed-gemm! (alpha A B beta C job)
+(definline complex-typed-gemm! (alpha A B beta C job)
(declare (type complex-matrix A B C)
(type complex-type alpha beta)
(type symbol job))
@@ -175,8 +175,16 @@
(if (and (member job-A '(:n :t))
(member job-B '(:n :t)))
(complex-base-typed-gemm! alpha A B beta C job)
- (let ((A (ecase job-A ((:h :c) (mconjugate A)) ((:n :t) A)))
- (B (ecase job-B ((:h :c) (mconjugate B)) ((:n :t) B)))
+ (let ((A (ecase job-A
+ ((:n :t) A)
+ ((:h :c)
+ (let ((ret (apply #'make-complex-tensor (lvec->list (dimensions A)))))
+ (complex-typed-axpy! #c(-1d0 0d0) A ret)))))
+ (B (ecase job-B
+ ((:n :t) B)
+ ((:h :c)
+ (let ((ret (apply #'make-complex-tensor (lvec->list (dimensions B)))))
+ (complex-typed-axpy! #c(-1d0 0d0) B ret)))))
(tjob (combine-jobs (ecase job-A ((:n :t) job-A) (:h :t) (:c :n))
(ecase job-B ((:n :t) job-B) (:h :t) (:c :n)))))
(complex-base-typed-gemm! alpha A B
commit 9b50685d6952b3be9ff29473595b2694ea234b08
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Aug 4 18:18:30 2012 +0530
o Added element-wise scal and div methods.
diff --git a/src/level-1/scal.lisp b/src/level-1/scal.lisp
index 15dfde7..2687144 100644
--- a/src/level-1/scal.lisp
+++ b/src/level-1/scal.lisp
@@ -28,7 +28,63 @@
(in-package #:matlisp)
-(defmacro generate-typed-scal! (func (tensor-class blas-func fortran-lb))
+(defmacro generate-typed-scal! (func (tensor-class fortran-func fortran-lb))
+ (let* ((opt (get-tensor-class-optimization tensor-class)))
+ (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
+ `(defun ,func (from to)
+ (declare (type ,tensor-class from to))
+ (let ((strd-p (blas-copyable-p from to))
+ (call-fortran? (> (number-of-elements to) ,fortran-lb)))
+ (cond
+ ((and strd-p call-fortran?)
+ (,fortran-func (number-of-elements from)
+ (store from) (first strd-p)
+ (store to) (second strd-p)
+ (head from) (head to)))
+ (t
+ (let ((f-sto (store from))
+ (t-sto (store to)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions from))
+ with (linear-sums
+ (f-of (strides from) (head from))
+ (t-of (strides to) (head to)))
+ do (let*-typed ((val-f ,(funcall (getf opt :reader) 'f-sto 'f-of) :type ,(getf opt :element-type))
+ (val-t ,(funcall (getf opt :reader) 't-sto 't-of) :type ,(getf opt :element-type))
+ (mul (* val-f val-t) :type ,(getf opt :element-type)))
+ ,(funcall (getf opt :value-writer) 'mul 't-sto 't-of))))))))
+ to)))
+
+(defmacro generate-typed-div! (func (tensor-class fortran-func fortran-lb))
+ (let* ((opt (get-tensor-class-optimization tensor-class)))
+ (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
+ `(defun ,func (from to)
+ (declare (type ,tensor-class from to))
+ (let ((strd-p (blas-copyable-p from to))
+ (call-fortran? (> (number-of-elements to) ,fortran-lb)))
+ (cond
+ ((and strd-p call-fortran?)
+ (,fortran-func (number-of-elements from)
+ (store from) (first strd-p)
+ (store to) (second strd-p)
+ (head from) (head to)))
+ (t
+ (let ((f-sto (store from))
+ (t-sto (store to)))
+ (declare (type ,(linear-array-type (getf opt :store-type)) f-sto t-sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions from))
+ with (linear-sums
+ (f-of (strides from) (head from))
+ (t-of (strides to) (head to)))
+ do (let*-typed ((val-f ,(funcall (getf opt :reader) 'f-sto 'f-of) :type ,(getf opt :element-type))
+ (val-t ,(funcall (getf opt :reader) 't-sto 't-of) :type ,(getf opt :element-type))
+ (mul (/ val-t val-f) :type ,(getf opt :element-type)))
+ ,(funcall (getf opt :value-writer) 'mul 't-sto 't-of))))))))
+ to)))
+
+(defmacro generate-typed-num-scal! (func (tensor-class blas-func fortran-lb))
(let ((opt (get-tensor-class-optimization tensor-class)))
(assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class)
`(defun ,func (alpha to)
@@ -51,17 +107,29 @@
to)))
;;Real
-(generate-typed-scal! real-typed-scal!
+(generate-typed-num-scal! real-typed-num-scal!
(real-tensor dscal *real-l1-fcall-lb*))
+(generate-typed-scal! real-typed-scal!
+ (real-tensor descal *real-l1-fcall-lb*))
+
+(generate-typed-div! real-typed-div!
+ (real-tensor dediv *real-l1-fcall-lb*))
+
;;Complex
(definline zordscal (nele alpha x incx &optional hd-x)
(if (zerop (imagpart alpha))
(zdscal nele (realpart alpha) x incx hd-x)
(zscal nele alpha x incx hd-x)))
-(generate-typed-scal! complex-typed-scal!
+(generate-typed-num-scal! complex-typed-num-scal!
(complex-tensor zordscal *complex-l1-fcall-lb*))
+
+(generate-typed-scal! complex-typed-scal!
+ (complex-tensor zescal *complex-l1-fcall-lb*))
+
+(generate-typed-scal! complex-typed-div!
+ (complex-tensor zediv *complex-l1-fcall-lb*))
;;---------------------------------------------------------------;;
(defgeneric scal! (alpha x)
@@ -74,13 +142,65 @@
Purpose
=======
X <- alpha .* X
-"))
+")
+ (:method :before ((x standard-tensor) (y standard-tensor))
+ (assert (lvec-eq (dimensions x) (dimensions y) #'=) nil
+ 'tensor-dimension-mismatch)))
(defmethod scal! ((alpha number) (x real-tensor))
- (real-typed-scal! (coerce-real alpha) x))
+ (real-typed-num-scal! (coerce-real alpha) x))
+
+(defmethod scal! ((x real-tensor) (y real-tensor))
+ (real-typed-scal! x y))
(defmethod scal! ((alpha number) (x complex-tensor))
- (complex-typed-scal! (coerce-complex alpha) x))
+ (complex-typed-num-scal! (coerce-complex alpha) x))
+
+(defmethod scal! ((x complex-tensor) (y complex-tensor))
+ (complex-typed-scal! x y))
+
+(defmethod scal! ((x real-tensor) (y complex-tensor))
+ (let ((tmp (tensor-realpart~ y)))
+ (real-typed-scal! x tmp)
+ ;;Move view to the imaginary part
+ (incf (head tmp))
+ (real-typed-scal! x tmp)))
+
+;;
+(defgeneric div! (alpha x)
+ (:documentation "
+ Syntax
+ ======
+ (div! alpha x)
+
+ Purpose
+ =======
+ X <- X ./ alpha
+
+ Yes the calling order is twisted.
+")
+ (:method :before ((x standard-tensor) (y standard-tensor))
+ (assert (lvec-eq (dimensions x) (dimensions y) #'=) nil
+ 'tensor-dimension-mismatch)))
+
+(defmethod div! ((alpha number) (x real-tensor))
+ (real-typed-num-scal! (coerce-real (/ 1 alpha)) x))
+
+(defmethod div! ((x real-tensor) (y real-tensor))
+ (real-typed-div! x y))
+
+(defmethod div! ((alpha number) (x complex-tensor))
+ (complex-typed-num-scal! (coerce-complex (/ 1 alpha)) x))
+
+(defmethod div! ((x complex-tensor) (y complex-tensor))
+ (complex-typed-div! x y))
+
+(defmethod div! ((x real-tensor) (y complex-tensor))
+ (let ((tmp (tensor-realpart~ y)))
+ (real-typed-div! x tmp)
+ ;;Move view to the imaginary part
+ (incf (head tmp))
+ (real-typed-div! x tmp)))
;;
(defgeneric scal (alpha x)
@@ -104,15 +224,68 @@
(* alpha x))
(defmethod scal ((alpha number) (x real-tensor))
- (let ((result (copy x)))
+ (let ((result (if (complexp alpha)
+ (copy! x (apply #'make-complex-tensor (lvec->list (dimensions x))))
+ (copy x))))
(scal! alpha result)))
-(defmethod scal ((alpha complex) (x real-tensor))
- (let* ((result (apply #'make-complex-tensor (lvec->list (dimensions x)))))
- (declare (type complex-tensor result))
- (copy! x result)
- (scal! alpha result)))
+(defmethod scal ((x real-tensor) (y real-tensor))
+ (scal! x (copy y)))
+
+(defmethod scal ((x complex-tensor) (y real-tensor))
+ (let ((result (copy! y (apply #'make-complex-tensor (lvec->list (dimensions x))))))
+ (scal! x result)))
(defmethod scal ((alpha number) (x complex-tensor))
(let ((result (copy x)))
(scal! alpha result)))
+
+(defmethod scal ((x standard-tensor) (y complex-tensor))
+ (let ((result (copy y)))
+ (scal! x result)))
+
+;;
+(defgeneric div (x y)
+ (:documentation "
+ Syntax
+ ======
+ (div! alpha x)
+
+ Purpose
+ =======
+ X <- X ./ alpha
+
+ Yes the calling order is twisted.
+"))
+
+(defmethod div ((alpha number) (x number))
+ (/ x alpha))
+
+(defmethod div ((alpha number) (x real-tensor))
+ (let ((result (if (complexp alpha)
+ (copy! x (apply #'make-complex-tensor (lvec->list (dimensions x))))
+ (copy x))))
+ (div! alpha result)))
+
+(defmethod div ((x real-tensor) (y real-tensor))
+ (div! x (copy y)))
+
+(defmethod div ((x complex-tensor) (y real-tensor))
+ (let ((result (copy! y (apply #'make-complex-tensor (lvec->list (dimensions x))))))
+ (div! x result)))
+
+(defmethod div ((alpha number) (x complex-tensor))
+ (let ((result (copy x)))
+ (div! alpha result)))
+
+(defmethod div ((x standard-tensor) (y complex-tensor))
+ (let ((result (copy y)))
+ (div! x result)))
+
+(defmethod div ((x real-tensor) (y (eql nil)))
+ (let ((result (copy! 1 (apply #'make-real-tensor (lvec->list (dimensions x))))))
+ (div! x result)))
+
+(defmethod div ((x complex-tensor) (y (eql nil)))
+ (let ((result (copy! 1 (apply #'make-complex-tensor (lvec->list (dimensions x))))))
+ (div! x result)))
commit 0b4fcdfe7d12f45c1d46f3b42589a5f2ff54e8dc
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Aug 4 16:22:25 2012 +0530
o Added escal ediv Fortran declarations.
diff --git a/matlisp.asd b/matlisp.asd
index 43f09d1..da36f7a 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -98,7 +98,8 @@
:depends-on ("foreign-interface")
:components ((:file "blas")
(:file "lapack")
- (:file "dfftpack")))
+ (:file "dfftpack")
+ (:file "libmatlisp")))
(:module "matlisp-base"
:depends-on ("foreign-core")
:pathname "base"
diff --git a/packages.lisp b/packages.lisp
index 1693289..8d71112 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -145,10 +145,17 @@
(:export #:zffti #:zfftf #:zfftb #:zffti #:zfftf #:zfftb)
(:documentation "FFT routines"))
+(defpackage "MATLISP-LIBMATLISP"
+ (:use #:common-lisp #:matlisp-ffi)
+ (:export
+ #:descal #:dediv
+ #:zescal #:zediv)
+ (:documentation "BLAS routines"))
+
(defpackage "MATLISP"
(:use #:common-lisp
#:matlisp-conditions #:matlisp-utilities #:matlisp-ffi
- #:matlisp-blas #:matlisp-lapack #:matlisp-dfftpack)
+ #:matlisp-blas #:matlisp-lapack #:matlisp-dfftpack #:matlisp-libmatlisp)
(:export #:index-type #:index-array #:allocate-index-store #:make-index-store
;;Standard-tensor
#:standard-tensor
diff --git a/src/foreign-core/libmatlisp.lisp b/src/foreign-core/libmatlisp.lisp
new file mode 100644
index 0000000..8e15a2f
--- /dev/null
+++ b/src/foreign-core/libmatlisp.lisp
@@ -0,0 +1,53 @@
+(in-package #:matlisp-libmatlisp)
+
+(def-fortran-routine descal :void
+ "
+ (DESCAL n dx incx dy incy)
+
+ Multiplies the vector X and Y element-wise.
+ Y <- Y .* E
+ "
+ (n :integer :input)
+ (dx (* :double-float :inc head-x) :input)
+ (incx :integer :input)
+ (dy (* :double-float :inc head-y) :output)
+ (incy :integer :output))
+
+(def-fortran-routine zescal :void
+ "
+ (ZESCAL n dx incx dy incy)
+
+ Multiplies the vector X and Y element-wise.
+ Y <- Y .* E
+ "
+ (n :integer :input)
+ (dx (* :complex-double-float :inc head-x) :input)
+ (incx :integer :input)
+ (dy (* :complex-double-float :inc head-y) :output)
+ (incy :integer :output))
+
+(def-fortran-routine dediv :void
+ "
+ (DEDIV n dx incx dy incy)
+
+ Divides the vector Y by X element-wise.
+ Y <- Y .* E
+ "
+ (n :integer :input)
+ (dx (* :double-float :inc head-x) :input)
+ (incx :integer :input)
+ (dy (* :double-float :inc head-y) :output)
+ (incy :integer :output))
+
+(def-fortran-routine zediv :void
+ "
+ (ZEDIV n dx incx dy incy)
+
+ Divide the vector Y by X element-wise.
+ Y <- Y .* E
+ "
+ (n :integer :input)
+ (dx (* :complex-double-float :inc head-x) :input)
+ (incx :integer :input)
+ (dy (* :complex-double-float :inc head-y) :output)
+ (incy :integer :output))
diff --git a/tests/loopy-tests.lisp b/tests/loopy-tests.lisp
index 2ed4ed1..67db2be 100644
--- a/tests/loopy-tests.lisp
+++ b/tests/loopy-tests.lisp
@@ -122,7 +122,7 @@
(hd-a (head t-a))
(hd-b (head t-b))
(hd-c (head t-c)))
- (declare (type (real-array *) st-a st-b st-c)
+ (declare (type real-store-vector st-a st-b st-c)
(type index-type rstrd-a cstrd-a rstrd-b cstrd-b rstrd-c cstrd-c nr-c
nc-c nc-a hd-a hd-b hd-c))
(mod-dotimes (idx (dimensions t-a))
-----------------------------------------------------------------------
Summary of changes:
lib/matlisp/dediv.f | 2 +-
lib/matlisp/zediv.f | 2 +-
matlisp.asd | 3 +-
packages.lisp | 9 ++-
src/foreign-core/libmatlisp.lisp | 53 +++++++++
src/level-1/axpy.lisp | 9 +-
src/level-1/copy.lisp | 4 +-
src/level-1/dot.lisp | 7 +-
src/level-1/realimag.lisp | 4 +-
src/level-1/scal.lisp | 239 +++++++++++++++++++++++++++++++++++--
src/level-1/swap.lisp | 2 +-
src/level-1/tensor-maker.lisp | 86 +++++++-------
src/level-1/trans.lisp | 14 +-
src/level-2/gemv.lisp | 14 ++-
src/level-3/gemm.lisp | 16 ++-
tests/loopy-tests.lisp | 2 +-
16 files changed, 379 insertions(+), 87 deletions(-)
create mode 100644 src/foreign-core/libmatlisp.lisp
hooks/post-receive
--
matlisp
|
