Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch classy updated. 2012-02-23-268-gc6c440e
|
From: Akshay S. <ak...@us...> - 2013-08-24 03:33:49
|
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, classy has been updated
via c6c440e0043ee6633cb729a0bc590e9ca97d5eff (commit)
via 6fb6102eb0a39b2bf48ba2dfbe98f2d7a1966935 (commit)
via 60f2b8490d5a31c90886c51d081bdb802f5431d6 (commit)
via 050548a939e93208a4990b6de248e3e39b3caa45 (commit)
via 7aa696273c37819ebc5b9ee1040d0f194dd8145a (commit)
via 570ae7eb80324580ee27cfa7ba1d20b11f779e41 (commit)
via 79a87b8605dba8ae97c8f354e42f1f081b127771 (commit)
via fc524fd099c95abfc3af0280e8a200e461cc9493 (commit)
via a3ea5898fd82c16b66bf0ce3d5615e370deb40a8 (commit)
via 69ca54a98c4e4a03e004268297644094b5541cae (commit)
via 6c6f96e88fab82f42a2cd563c53e90c48eb8da24 (commit)
via 1f3bc86e39b2b9d1a23946434486a99faf9d7eaf (commit)
via fd0546544cd3c21641688e03ff221b031ac01ae4 (commit)
via 62f126ef664982e0c8ffd132ade8bb5308833f56 (commit)
via 9afd000d1a6b497e3bf4fdc0318884b412773de7 (commit)
via a7bad8cec909a69bb312917406d3dfc1626f8c12 (commit)
via 831a7c79908907c3702b623d0eae3a0a1f746a58 (commit)
via 05cd941088d8c303f3b5f81d34a6fb336f1033a9 (commit)
via f051c33ed570af222cff1fbf93802dc8844034ad (commit)
via f314da645424005358b1921156218063d8ea64c2 (commit)
via f47066fa8877ca56d4d58f36a5a08515593b8a2f (commit)
via 1894d2f6b21c756ec5bb2ddd443e6d38a7087f2b (commit)
from 808353d428ddc07d365bf1de8abcc86f0179ee08 (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 c6c440e0043ee6633cb729a0bc590e9ca97d5eff
Author: Akshay Srinivasan <aks...@gm...>
Date: Fri Aug 23 20:33:50 2013 -0700
Added seq, map functions.
diff --git a/matlisp.asd b/matlisp.asd
index 44384f4..b4b9ac0 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -168,7 +168,9 @@
(:module "matlisp-special"
:pathname "special"
:depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
- :components ((:file "random")))
+ :components ((:file "random")
+ (:file "map")
+ (:file "seq")))
#+nil
(:module "matlisp-sugar"
:pathname "sugar"
diff --git a/src/base/generic-copy.lisp b/src/base/generic-copy.lisp
index 543fa5c..aa30346 100644
--- a/src/base/generic-copy.lisp
+++ b/src/base/generic-copy.lisp
@@ -78,6 +78,10 @@
y))
(copy! x y)))
+(defmethod copy! ((x cons) (y standard-tensor))
+ ;;You seriously weren't expecting efficiency were you :) ?
+ (let ((arr (make-array (list-dimensions x) :initial-contents x)))
+ (copy! arr y)))
;;
(defgeneric copy (object)
(:documentation
diff --git a/src/classes/numeric.lisp b/src/classes/numeric.lisp
index 6f3e87e..e266cfc 100644
--- a/src/classes/numeric.lisp
+++ b/src/classes/numeric.lisp
@@ -53,8 +53,9 @@
(deft/method t/l3-lb (sym complex-numeric-tensor) ()
'*complex-l3-fcall-lb*)
-;;SBCL uses specialized arrays for floating complex arrays.
-#-sbcl
+;;Comment this block if you want to use (simple-array (complex double-float) (*))
+;;as the underlying store. This will make Lisp-implementations of gemm .. faster
+;;but you'll lose the ability to use tensor-realpart~/imagpart~.
(progn
(deft/method t/store-element-type (sym complex-numeric-tensor) ()
(let ((cplx-type (macroexpand-1 `(t/field-type ,sym))))
diff --git a/src/level-1/map.lisp b/src/level-1/map.lisp
deleted file mode 100644
index 3359ea9..0000000
--- a/src/level-1/map.lisp
+++ /dev/null
@@ -1,72 +0,0 @@
-(in-package #:matlisp)
-
-(defgeneric mapsor! (func x y)
- (:documentation "
- y <- func(x)
-")
- (:method :before ((func function) (x standard-tensor) (y standard-tensor))
- (assert (very-quickly (lvec-eq (dimensions x) (dimensions y))) nil 'tensor-dimension-mismatch)))
-
-(defmethod mapsor! ((func function) (x standard-tensor) (y standard-tensor))
- (let ((clx (class-name (class-of x)))
- (cly (class-name (class-of y))))
- (assert (and
- (member clx *tensor-type-leaves*)
- (member cly *tensor-type-leaves*))
- nil 'tensor-abstract-class :tensor-class (list clx cly))
- (compile-and-eval
- `(defmethod mapsor! ((func function) (x ,clx) (y ,cly))
- (let ((sto-x (store x))
- (sto-y (store y)))
- (declare (type ,(store-type clx) sto-x)
- (type ,(store-type cly) sto-y))
- (very-quickly
- (mod-dotimes (idx (dimensions x))
- :with (linear-sums
- (of-x (strides x))
- (of-y (strides y)))
- :do (t/store-set ,cly (funcall func (t/store-ref ,clx sto-x of-x)) sto-y of-y))))
- y)))
- (mapsor! func x y))
-
-(defun mapcol (func x &optional (axis 1))
- (declare (type standard-tensor x))
- (assert (tensor-matrixp x) nil 'tensor-dimension-mismatch)
- (let* ((v-x (slice~ x axis))
- (st-x (aref (the index-store-vector (strides x)) axis))
- (ret (zeros (aref (the index-store-vector (dimensions x)) axis) (class-of x))))
- (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
- :do (progn
- (setf (ref ret i) (funcall func v-x))
- (incf (slot-value v-x 'head) st-x)))
- ret))
-
-(defun tensor-min (x)
- (let ((min nil))
- (mod-dotimes (idx (dimensions x))
- :do (when (or (null min) (< (ref x idx) min))
- (setf min (ref x idx))))
- min))
-
-(defun idx-min (v)
- (let ((min-idx nil)
- (min nil))
- (loop :for i :from 0 :below (aref (the index-store-vector (dimensions v)) 0)
- :do (when (or (null min) (< (ref v i) min))
- (setf min-idx i
- min (ref v i))))
- (values min-idx min)))
-
-(defun mapxis! (func x y &optional (axis 0))
- (declare (type standard-tensor x y))
- (multiple-value-bind (view-x view-y) (let ((slst (make-list (rank x) :initial-element '(* * *))))
- (rplaca (nthcdr axis slst) (list 0 '* 1))
- (values (sub-tensor~ x slst nil) (sub-tensor~ y slst nil)))
- (let ((st-x (aref (the index-store-vector (dimensions x)) axis))
- (st-y (aref (the index-store-vector (dimensions x)) axis)))
- (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
- :do (progn
- ;;May die if you're doing fancy stuff.
- (funcall func view-x view-y)
-
-
diff --git a/src/special/map.lisp b/src/special/map.lisp
new file mode 100644
index 0000000..9cb083a
--- /dev/null
+++ b/src/special/map.lisp
@@ -0,0 +1,70 @@
+(in-package #:matlisp)
+
+(defgeneric mapsor! (func x y)
+ (:documentation "
+ Syntax
+ ======
+ (MAPSOR! func x y)
+
+ Purpose
+ =======
+ Applies the function element-wise on x, and sets the corresponding
+ elements in y to the value returned by the function.
+
+ Example
+ =======
+ > (mapsor! #'sin (randn '(2 2)) (zeros '(2 2)))
+")
+ (:method :before ((func function) (x standard-tensor) (y standard-tensor))
+ (assert (very-quickly (lvec-eq (dimensions x) (dimensions y))) nil 'tensor-dimension-mismatch)))
+
+(defmethod mapsor! ((func function) (x standard-tensor) (y standard-tensor))
+ (let ((clx (class-name (class-of x)))
+ (cly (class-name (class-of y))))
+ (assert (and
+ (member clx *tensor-type-leaves*)
+ (member cly *tensor-type-leaves*))
+ nil 'tensor-abstract-class :tensor-class (list clx cly))
+ (compile-and-eval
+ `(defmethod mapsor! ((func function) (x ,clx) (y ,cly))
+ (let ((sto-x (store x))
+ (sto-y (store y)))
+ (declare (type ,(store-type clx) sto-x)
+ (type ,(store-type cly) sto-y))
+ (very-quickly
+ (mod-dotimes (idx (dimensions x))
+ :with (linear-sums
+ (of-x (strides x))
+ (of-y (strides y)))
+ :do (t/store-set ,cly (funcall func (t/store-ref ,clx sto-x of-x)) sto-y of-y))))
+ y)))
+ (mapsor! func x y))
+;;
+(defun mapslice (func x &optional (axis 1))
+ (declare (type standard-tensor x))
+ (let* ((v-x (slice~ x axis))
+ (st-x (aref (strides x) axis)))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :collect (prog1 (funcall func (copy v-x))
+ (incf (slot-value v-x 'head) st-x)))))
+
+(defmacro tensor-foldl (type func ten init &optional (init-type (field-type type)))
+ (using-gensyms (decl (ten init))
+ (with-gensyms (sto idx of funcsym)
+ `(let* (,@decl
+ ,@(unless (symbolp func)
+ `((,funcsym ,func)))
+ (,sto (store ,ten)))
+ (declare (type ,type ,ten)
+ ,@(unless (symbolp func) `((type function ,funcsym)))
+ (type ,(store-type type) ,sto)
+ ,@(when init-type
+ `((type ,init-type ,init))))
+ (very-quickly
+ (mod-dotimes (,idx (dimensions ,ten))
+ :with (linear-sums
+ (,of (strides ,ten)))
+ :do (setf ,init (,@(if (symbolp func)
+ `(,func)
+ `(funcall ,funcsym)) ,init (t/store-ref ,type ,sto ,of)))))
+ ,init))))
diff --git a/src/sugar/seq.lisp b/src/special/seq.lisp
similarity index 70%
rename from src/sugar/seq.lisp
rename to src/special/seq.lisp
index 9c96efa..41ef589 100644
--- a/src/sugar/seq.lisp
+++ b/src/special/seq.lisp
@@ -27,22 +27,26 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package #:matlisp)
-(defun arange (start end &optional (h 1d0))
+(defun range (start end &optional (h 1d0))
(let ((quo (ceiling (if (> start end) (- start end) (- end start)) h)))
(if (= quo 0) nil
- (let*-typed ((ret (real-typed-zeros (idxv quo)) :type real-tensor)
- (sto-r (store ret) :type real-store-vector)
- (h (coerce-real-unforgiving (if (> start end) (- h) h)) :type real-type))
- (loop :for i :from 0 :below quo
- :for ori := (coerce-real-unforgiving start) :then (+ ori h)
- :do (setf (aref sto-r i) ori))
- ret))))
+ (let* ((ret (zeros quo 'real-tensor))
+ (sto (store ret))
+ (h (coerce (if (> start end) (- h) h) 'double-float)))
+ (declare (type (simple-array double-float (*)) sto)
+ (type double-float h))
+ (very-quickly
+ (loop :for i :from 0 :below quo
+ :for ori := (coerce start 'double-float) :then (+ ori h)
+ :do (t/store-set real-tensor ori sto i)))
+ ret))))
-(defun alinspace (start end &optional (num-points (1+ (abs (- start end)))))
- (let ((h (/ (- end start) (1- num-points))))
- (arange start (+ h end) (abs h))))
+(defun linspace (start end &optional (num-points (1+ (abs (- start end)))))
+ (let* ((num-points (floor num-points))
+ (h (/ (- end start) (1- num-points))))
+ (range start (+ h end) (abs h))))
-(defun range (start end &optional (h 1))
+(defun list-range (start end &optional (h 1))
(declare (type real start end h))
(let ((quo (ceiling (if (> start end) (- start end) (- end start)) h)))
(if (= quo 0) nil
@@ -51,6 +55,8 @@
:for ori := start :then (+ ori h)
:collect ori)))))
-(defun linspace (start end &optional (num-points (1+ (abs (- start end)))))
- (let ((h (/ (- end start) (1- num-points))))
- (range start (+ h end) (abs h))))
+(defun list-linspace (start end &optional (num-points (1+ (abs (- start end)))))
+ (let* ((num-points (floor num-points))
+ (h (/ (- end start) (1- num-points))))
+ (list-range start (+ h end) (abs h))))
+
diff --git a/src/utilities/lvec.lisp b/src/utilities/lvec.lisp
index 3ab1918..5e81a47 100644
--- a/src/utilities/lvec.lisp
+++ b/src/utilities/lvec.lisp
@@ -5,7 +5,7 @@
(declare (type vector))
(loop
:for i :of-type fixnum :from 0 :below (length vec)
- :for ret = (aref vec 0) :then (funcall func (aref vec i) ret)
+ :for ret = (aref vec 0) :then (funcall func ret (aref vec i))
:finally (return ret)))
(definline lvec-foldr (func vec)
commit 6fb6102eb0a39b2bf48ba2dfbe98f2d7a1966935
Author: Akshay Srinivasan <aks...@gm...>
Date: Thu Aug 22 21:58:42 2013 -0700
Cleaned up random.
diff --git a/AUTHORS b/AUTHORS
index 5f68edf..054d16e 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -9,3 +9,6 @@ Femlisp (www.femlisp.org); it has been used here
The infix reader is modified and included here with the
permission of its original author Mark Kantrowitz.
+
+The random number generators are borrowed and adapted from
+cl-random written by Tamas Papp.
\ No newline at end of file
diff --git a/matlisp.asd b/matlisp.asd
index cfb71d3..44384f4 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -165,6 +165,10 @@
:components ((:file "lu")
(:file "chol")
(:file "eig")))
+ (:module "matlisp-special"
+ :pathname "special"
+ :depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
+ :components ((:file "random")))
#+nil
(:module "matlisp-sugar"
:pathname "sugar"
diff --git a/src/packages/odepack/dlsode.lisp b/src/packages/odepack/dlsode.lisp
index dab15bb..466603c 100644
--- a/src/packages/odepack/dlsode.lisp
+++ b/src/packages/odepack/dlsode.lisp
@@ -68,13 +68,12 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun pend-field (neq time y ydot)
- (setf (fv-ref ydot 0) (fv-ref y 1)
- (fv-ref ydot 1) (- (sin (fv-ref y 0)))))
+ (setf (fv-ref ydot 0) (fv-ref y 1)
+ (fv-ref ydot 1) (- (sin (fv-ref y 0)))))
(defun pend-report (ts y)
(format t "~A ~A ~A ~%" ts (aref y 0) (aref y 1)))
-
(defun pcart-field (neq time y ydot)
(declare (ignore neq time))
(very-quickly
diff --git a/src/level-1/random.lisp b/src/special/random.lisp
similarity index 50%
rename from src/level-1/random.lisp
rename to src/special/random.lisp
index 4bdf231..997d039 100644
--- a/src/level-1/random.lisp
+++ b/src/special/random.lisp
@@ -1,17 +1,18 @@
(in-package #:matlisp)
-(declaim (inline draw-standard-exponential))
-(defun draw-standard-exponential ()
+(declaim (ftype (function () double-float) draw-standard-normal draw-standard-normal-marsaglia draw-standard-exponential))
+
+(definline draw-standard-exponential ()
"Return a random variable from the Exponential(1) distribution, which has density exp(-x)."
+ ;; Adapted from cl-random, originally written by Tamas Papp
;; need 1-random, because there is a small but nonzero chance of getting a 0.
(- (log (- 1d0 (random 1d0)))))
-(declaim (ftype (function () double-float) draw-standard-normal-leva draw-standard-normal-marsaglia))
-(definline draw-standard-normal-leva ()
+(definline draw-standard-normal ()
"Draw a random number from N(0,1)."
;; Method from Leva (1992). This is considered much better/faster than the Box-Muller method.
;; Adapted from cl-random, originally written by Tamas Papp
- ;; This tends to be just as fast as Marsaglia with storage.
+ ;; This seems to be just as fast as Marsaglia with storage.
(very-quickly
(loop
:do (let* ((u (random 1d0))
@@ -43,25 +44,33 @@
(setf prev (* y mult))
(return (* x mult))))))))))
-;;
-(defun randn (dims)
- (let* ((ret (zeros dims 'real-tensor))
- (sto (store ret)))
- (declare (type (simple-array double-float (*)) sto))
- (very-quickly
- (mod-dotimes (idx (dimensions ret))
- :with (linear-sums
- (of-ret (strides ret)))
- :do (setf (aref sto of-ret) (the double-float (draw-standard-normal-leva)))))
- ret))
+;;
+(defmacro fill-tensor (type (func tensor))
+ (using-gensyms (decl (tensor))
+ (with-gensyms (sto ofst)
+ `(let* (,@decl
+ (,sto (store ,tensor)))
+ (declare (type ,type ,tensor)
+ (type ,(store-type type) ,sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions ,tensor))
+ :with (linear-sums
+ (,ofst (strides ,tensor)))
+ :do (t/store-set ,type ,(etypecase func (symbol `(,func)) (cons func)) ,sto ,ofst)))
+ ,tensor))))
+
+(macrolet ((generate-rand (func clause)
+ (let ((clause (etypecase clause
+ (symbol `(,clause))
+ (cons clause))))
+ `(defun ,func (&optional dims)
+ (if dims
+ (fill-tensor real-tensor (,clause (zeros dims 'real-tensor)))
+ ,clause))))
+ (generate-rands ((&rest args))
+ `(progn
+ ,@(mapcar #'(lambda (x) `(generate-rand ,(car x) ,(cadr x))) args))))
+ (generate-rands ((randn (draw-standard-normal))
+ (rand (random 1d0))
+ (rande (draw-standard-exponential)))))
-(defun rand (dims)
- (let* ((ret (zeros dims 'real-tensor))
- (sto (store ret)))
- (declare (type (simple-array double-float (*)) sto))
- (very-quickly
- (mod-dotimes (idx (dimensions ret))
- :with (linear-sums
- (of-ret (strides ret)))
- :do (setf (aref sto of-ret) (random 1d0))))
- ret))
commit 60f2b8490d5a31c90886c51d081bdb802f5431d6
Author: Akshay Srinivasan <aks...@gm...>
Date: Wed Aug 21 19:50:19 2013 -0700
Made attribute hash-table allocation lazy.
diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp
index 9bfaf56..d34a360 100644
--- a/lib-src/gnuplot/gnuplot.lisp
+++ b/lib-src/gnuplot/gnuplot.lisp
@@ -8,7 +8,7 @@
ccl:run-program
gnuplot-binary nil :input :stream :wait nil :output t))
-(defun plot2d (data &key (color (list "#FF0000")))
+(defun gnuplot-send (str &rest args)
(unless *current-gnuplot-process*
(setf *current-gnuplot-process* (open-gnuplot-stream)))
(let ((stream (#+:sbcl
@@ -16,23 +16,16 @@
#+:ccl
ccl:external-process-input-stream
*current-gnuplot-process*)))
- (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create)
- (loop :for i :from 0 :below (loop :for x :in data :minimizing (number-of-elements x))
- :do (loop :for x :in data :do (format s "~a " (tensor-ref x i)) :finally (format s "~%"))))
- (format stream "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color)
- (finish-output stream)))
-
-(defun gnuplot-send (str)
- (unless *current-gnuplot-process*
- (setf *current-gnuplot-process* (open-gnuplot-stream)))
- (let ((stream (#+:sbcl
- sb-ext:process-input
- #+:ccl
- ccl:external-process-input-stream
- *current-gnuplot-process*)))
- (format stream "~a~%" str)
+ (apply #'format (append (list stream str) args))
(finish-output stream)))
+(defun plot2d (data &key (lines t) (color (list "#FF0000")))
+ (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create)
+ (loop :for i :from 0 :below (loop :for x :in data :minimizing (size x))
+ :do (loop :for x :in data :do (format s "~a " (coerce (ref x i) 'single-float)) :finally (format s "~%"))))
+ (if lines
+ (gnuplot-send "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color)
+ (gnuplot-send "plot '/tmp/matlisp-gnuplot.out'~%")))
;; (defclass gnuplot-plot-info ()
;; ((title
diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp
index c34ed6d..426536f 100644
--- a/src/base/standard-tensor.lisp
+++ b/src/base/standard-tensor.lisp
@@ -58,10 +58,22 @@
(store :initarg :store :reader store
:documentation "The actual storage for the tensor.")
;;
- (attributes :initarg :attributes :reader attributes :initform (make-hash-table)
+ (attributes :initarg :attributes :initform nil
:documentation "Place for computable attributes of an object instance."))
(:documentation "Basic tensor class."))
+;;Create hash-table only when necessary
+(definline attributes (x)
+ (declare (type standard-tensor x))
+ (or (slot-value x 'attributes)
+ (let ((htbl (make-hash-table)))
+ (setf (slot-value x 'attributes) htbl)
+ htbl)))
+
+(declaim (ftype (function (standard-tensor) index-store-vector) strides dimensions)
+ (ftype (function (standard-tensor) index-type) head)
+ (ftype (function (standard-tensor) hash-table) attributes))
+
(defmacro memoizing ((tensor name) &rest body)
(declare (type symbol name))
(with-gensyms (tens)
@@ -95,6 +107,7 @@
(declare (type standard-tensor tensor))
(lvec-foldr #'* (the index-store-vector (dimensions tensor))))
+;;
(defgeneric store-size (tensor)
(:documentation "
Syntax
@@ -424,3 +437,14 @@
:strides (very-quickly (vectorify (the index-store-vector nstds) nrank 'index-type))
:store (store tensor)
:parent-tensor tensor)))))))
+
+(definline slice~ (x axis &optional (idx 0))
+ (let ((slst (make-list (rank x) :initial-element '(* * *))))
+ (rplaca (nthcdr axis slst) (list idx '* (1+ idx)))
+ (sub-tensor~ x slst nil)))
+
+(definline row-slice~ (x idx)
+ (slice~ x 0 idx))
+
+(definline col-slice~ (x idx)
+ (slice~ x 1 idx))
diff --git a/src/level-1/map.lisp b/src/level-1/map.lisp
new file mode 100644
index 0000000..3359ea9
--- /dev/null
+++ b/src/level-1/map.lisp
@@ -0,0 +1,72 @@
+(in-package #:matlisp)
+
+(defgeneric mapsor! (func x y)
+ (:documentation "
+ y <- func(x)
+")
+ (:method :before ((func function) (x standard-tensor) (y standard-tensor))
+ (assert (very-quickly (lvec-eq (dimensions x) (dimensions y))) nil 'tensor-dimension-mismatch)))
+
+(defmethod mapsor! ((func function) (x standard-tensor) (y standard-tensor))
+ (let ((clx (class-name (class-of x)))
+ (cly (class-name (class-of y))))
+ (assert (and
+ (member clx *tensor-type-leaves*)
+ (member cly *tensor-type-leaves*))
+ nil 'tensor-abstract-class :tensor-class (list clx cly))
+ (compile-and-eval
+ `(defmethod mapsor! ((func function) (x ,clx) (y ,cly))
+ (let ((sto-x (store x))
+ (sto-y (store y)))
+ (declare (type ,(store-type clx) sto-x)
+ (type ,(store-type cly) sto-y))
+ (very-quickly
+ (mod-dotimes (idx (dimensions x))
+ :with (linear-sums
+ (of-x (strides x))
+ (of-y (strides y)))
+ :do (t/store-set ,cly (funcall func (t/store-ref ,clx sto-x of-x)) sto-y of-y))))
+ y)))
+ (mapsor! func x y))
+
+(defun mapcol (func x &optional (axis 1))
+ (declare (type standard-tensor x))
+ (assert (tensor-matrixp x) nil 'tensor-dimension-mismatch)
+ (let* ((v-x (slice~ x axis))
+ (st-x (aref (the index-store-vector (strides x)) axis))
+ (ret (zeros (aref (the index-store-vector (dimensions x)) axis) (class-of x))))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :do (progn
+ (setf (ref ret i) (funcall func v-x))
+ (incf (slot-value v-x 'head) st-x)))
+ ret))
+
+(defun tensor-min (x)
+ (let ((min nil))
+ (mod-dotimes (idx (dimensions x))
+ :do (when (or (null min) (< (ref x idx) min))
+ (setf min (ref x idx))))
+ min))
+
+(defun idx-min (v)
+ (let ((min-idx nil)
+ (min nil))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions v)) 0)
+ :do (when (or (null min) (< (ref v i) min))
+ (setf min-idx i
+ min (ref v i))))
+ (values min-idx min)))
+
+(defun mapxis! (func x y &optional (axis 0))
+ (declare (type standard-tensor x y))
+ (multiple-value-bind (view-x view-y) (let ((slst (make-list (rank x) :initial-element '(* * *))))
+ (rplaca (nthcdr axis slst) (list 0 '* 1))
+ (values (sub-tensor~ x slst nil) (sub-tensor~ y slst nil)))
+ (let ((st-x (aref (the index-store-vector (dimensions x)) axis))
+ (st-y (aref (the index-store-vector (dimensions x)) axis)))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :do (progn
+ ;;May die if you're doing fancy stuff.
+ (funcall func view-x view-y)
+
+
diff --git a/src/level-1/random.lisp b/src/level-1/random.lisp
new file mode 100644
index 0000000..4bdf231
--- /dev/null
+++ b/src/level-1/random.lisp
@@ -0,0 +1,67 @@
+(in-package #:matlisp)
+
+(declaim (inline draw-standard-exponential))
+(defun draw-standard-exponential ()
+ "Return a random variable from the Exponential(1) distribution, which has density exp(-x)."
+ ;; need 1-random, because there is a small but nonzero chance of getting a 0.
+ (- (log (- 1d0 (random 1d0)))))
+
+(declaim (ftype (function () double-float) draw-standard-normal-leva draw-standard-normal-marsaglia))
+(definline draw-standard-normal-leva ()
+ "Draw a random number from N(0,1)."
+ ;; Method from Leva (1992). This is considered much better/faster than the Box-Muller method.
+ ;; Adapted from cl-random, originally written by Tamas Papp
+ ;; This tends to be just as fast as Marsaglia with storage.
+ (very-quickly
+ (loop
+ :do (let* ((u (random 1d0))
+ (v (* 1.7156d0 (- (random 1d0) 0.5d0)))
+ (x (- u 0.449871d0))
+ (y (+ (abs v) 0.386595d0))
+ (q (+ (expt x 2) (* y (- (* 0.19600d0 y) (* 0.25472d0 x))))))
+ (declare (type double-float u v x y q))
+ (unless (and (> q 0.27597d0)
+ (or (> q 0.27846d0)
+ (plusp (+ (expt v 2) (* 4 (expt u 2) (log u))))))
+ (return (/ v u)))))))
+
+;;Not thread safe, obviously
+(let ((prev nil))
+ (defun draw-standard-normal-marsaglia ()
+ (if prev
+ (prog1 prev
+ (setf prev nil))
+ (very-quickly
+ (loop
+ :do (let* ((x (1- (random 2d0)))
+ (y (1- (random 2d0)))
+ (s (+ (* x x) (* y y))))
+ (declare (type double-float x y s))
+ (when (<= s 1d0)
+ (let ((mult (sqrt (/ (* -2 (log s)) s))))
+ (declare (type double-float mult))
+ (setf prev (* y mult))
+ (return (* x mult))))))))))
+
+;;
+(defun randn (dims)
+ (let* ((ret (zeros dims 'real-tensor))
+ (sto (store ret)))
+ (declare (type (simple-array double-float (*)) sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions ret))
+ :with (linear-sums
+ (of-ret (strides ret)))
+ :do (setf (aref sto of-ret) (the double-float (draw-standard-normal-leva)))))
+ ret))
+
+(defun rand (dims)
+ (let* ((ret (zeros dims 'real-tensor))
+ (sto (store ret)))
+ (declare (type (simple-array double-float (*)) sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions ret))
+ :with (linear-sums
+ (of-ret (strides ret)))
+ :do (setf (aref sto of-ret) (random 1d0))))
+ ret))
commit 050548a939e93208a4990b6de248e3e39b3caa45
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Aug 20 01:21:26 2013 -0700
Added sum method. Made changes to dot.
diff --git a/matlisp.asd b/matlisp.asd
index f3c236b..cfb71d3 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -148,7 +148,9 @@
(:file "realimag"
:depends-on ("copy"))
(:file "trans"
- :depends-on ("scal" "copy"))))
+ :depends-on ("scal" "copy"))
+ (:file "sum"
+ :depends-on ("dot" "copy"))))
(:module "matlisp-level-2"
:pathname "level-2"
:depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1")
diff --git a/src/level-1/dot.lisp b/src/level-1/dot.lisp
index 354c271..3a81f7e 100644
--- a/src/level-1/dot.lisp
+++ b/src/level-1/dot.lisp
@@ -34,41 +34,48 @@
(if conjp 'zdotc 'zdotu))
;;
-(deft/generic (t/blas-dot #'subtypep) sym (x y &optional conjp))
-(deft/method t/blas-dot (sym blas-numeric-tensor) (x y &optional (conjp t))
- (using-gensyms (decl (x y))
- `(let (,@decl)
- (declare (type ,sym ,x ,y))
- (,(macroexpand-1 `(t/blas-dot-func ,sym ,conjp))
- (aref (the index-store-vector (dimensions ,x)) 0)
- (the ,(store-type sym) (store ,x)) (aref (the index-store-vector (strides ,x)) 0)
- (the ,(store-type sym) (store ,y)) (aref (the index-store-vector (strides ,y)) 0)
- (head ,x) (head ,y)))))
+(deft/generic (t/blas-dot #'subtypep) sym (x y &optional conjp num-y?))
+(deft/method t/blas-dot (sym blas-numeric-tensor) (x y &optional (conjp t) (num-y? nil))
+ (using-gensyms (decl (x y))
+ (with-gensyms (sto)
+ `(let (,@decl
+ ,@(when num-y? `((,sto (t/store-allocator ,sym 1)))))
+ (declare (type ,sym ,x ,@(unless num-y? `(,y)))
+ ,@(when num-y? `((type ,(field-type sym) ,y)
+ (type ,(store-type sym) ,sto))))
+ ,@(when num-y? `((t/store-set ,sym ,y ,sto 0)))
+ (,(macroexpand-1 `(t/blas-dot-func ,sym ,conjp))
+ (aref (the index-store-vector (dimensions ,x)) 0)
+ (the ,(store-type sym) (store ,x)) (aref (the index-store-vector (strides ,x)) 0)
+ (the ,(store-type sym) ,(if num-y? sto `(store ,y))) ,(if num-y? 0 `(aref (the index-store-vector (strides ,y)) 0))
+ (head ,x) ,(if num-y? 0 `(head ,y)))))))
-(deft/generic (t/dot #'subtypep) sym (x y &optional conjp))
-(deft/method t/dot (sym standard-tensor) (x y &optional (conjp t))
+(deft/generic (t/dot #'subtypep) sym (x y &optional conjp num-y?))
+(deft/method t/dot (sym standard-tensor) (x y &optional (conjp t) (num-y? nil))
(using-gensyms (decl (x y))
(with-gensyms (sto-x sto-y of-x of-y stp-x stp-y dot)
`(let (,@decl)
- (declare (type ,sym ,x ,y))
+ (declare (type ,sym ,x ,@(unless num-y? `(,y)))
+ ,@(when num-y? `((type ,(field-type sym) ,y))))
(let ((,sto-x (store ,x))
(,stp-x (aref (the index-store-vector (strides ,x)) 0))
(,of-x (head ,x))
- (,sto-y (store ,y))
- (,stp-y (aref (the index-store-vector (strides ,y)) 0))
- (,of-y (head ,y))
+ ,@(unless num-y?
+ `((,sto-y (store ,y))
+ (,stp-y (aref (the index-store-vector (strides ,y)) 0))
+ (,of-y (head ,y))))
(,dot (t/fid+ ,(field-type sym))))
- (declare (type ,(store-type sym) ,sto-x ,sto-y)
- (type index-type ,stp-x ,stp-y ,of-x ,of-y)
+ (declare (type ,(store-type sym) ,sto-x ,@(unless `(,sto-y)))
+ (type index-type ,stp-x ,of-x ,@(unless num-y? `(,stp-y ,of-y)))
(type ,(field-type sym) ,dot))
(very-quickly
(loop :repeat (aref (the index-store-vector (dimensions ,x)) 0)
:do (setf ,dot (t/f+ ,(field-type sym) ,dot
(t/f* ,(field-type sym)
,(recursive-append (when conjp `(t/fc ,(field-type sym))) `(t/store-ref ,sym ,sto-x ,of-x))
- (t/store-ref ,sym ,sto-y ,of-y)))
+ ,(if num-y? y `(t/store-ref ,sym ,sto-y ,of-y))))
,of-x (+ ,of-x ,stp-x)
- ,of-y (+ ,of-y ,stp-y))))
+ ,@(unless num-y? `(,of-y (+ ,of-y ,stp-y))))))
,dot)))))
;;---------------------------------------------------------------;;
(defgeneric dot (x y &optional conjugate-p)
@@ -133,3 +140,24 @@
(dot x y conjugate-p))
(t
(error "Don't know how to compute the dot product of ~a , ~a." clx cly)))))
+
+(defmethod dot ((x standard-tensor) (y t) &optional (conjugate-p t))
+ (let ((clx (class-name (class-of x))))
+ (assert (member clx *tensor-type-leaves*)
+ nil 'tensor-abstract-class :tensor-class (list clx))
+ (compile-and-eval
+ `(defmethod dot ((x ,clx) (y t) &optional (conjugate-p t))
+ (let ((y (t/coerce ,(field-type clx) y)))
+ (declare (type ,(field-type clx) y))
+ ,(recursive-append
+ (when (subtypep clx 'blas-numeric-tensor)
+ `(if (call-fortran? x (t/l1-lb ,clx))
+ (if conjugate-p
+ (t/blas-dot ,clx x y t t)
+ (t/blas-dot ,clx x y nil t))))
+ `(if conjugate-p
+ ;;Please do your checks before coming here.
+ (t/dot ,clx x y t t)
+ (t/dot ,clx x y nil t))))))
+ (dot x y conjugate-p)))
+
diff --git a/src/level-1/sum.lisp b/src/level-1/sum.lisp
new file mode 100644
index 0000000..dc6dada
--- /dev/null
+++ b/src/level-1/sum.lisp
@@ -0,0 +1,65 @@
+(in-package #:matlisp)
+
+(deft/generic (t/sum #'subtypep) sym (x ret &optional axis))
+(deft/method t/sum (sym standard-tensor) (x ret &optional (axis 0))
+ (if (null ret)
+ `(dot ,x (t/fid* ,(field-type sym)) nil)
+ (using-gensyms (decl (x axis ret))
+ (with-gensyms (view argstd)
+ `(let* (,@decl)
+ (declare (type ,sym ,x ,ret)
+ (type index-type ,axis))
+ (let ((,view (let ((slst (make-list (rank ,x) :initial-element '(* * *))))
+ (rplaca (nthcdr ,axis slst) (list 0 '* 1))
+ (sub-tensor~ ,x slst nil)))
+ (,argstd (aref (the index-store-vector (strides ,x)) ,axis)))
+ (declare (type ,sym ,view)
+ (type index-type ,argstd))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions ,x)) ,axis)
+ :do (progn
+ (axpy! (t/fid* ,(field-type sym)) ,view ,ret)
+ (incf (slot-value ,view 'head) ,argstd)))
+ ,ret))))))
+
+(defgeneric sum! (x y &optional axis)
+ (:documentation "
+ (SUM! x y [axis 0])
+
+ --
+ y <- \ x(:, : ..., i, :, :..)
+ /_
+ i
+ where the index to be summed over is chosen using @arg{axis}.
+")
+ (:method :before ((x standard-tensor) (y standard-tensor) &optional (axis 0))
+ (assert (and
+ (= (1- (rank x)) (rank y))
+ (let ((dims-x (dimensions x))
+ (dims-y (dimensions y)))
+ (declare (type index-store-vector dims-x dims-y))
+ (loop
+ :for i :from 0 :below (rank x)
+ :and j := 0 :then (if (= i axis) j (1+ j))
+ :do (unless (or (= i axis) (= (aref dims-x i) (aref dims-y j)))
+ (return nil))
+ :finally (return t))))
+ nil 'tensor-dimension-mismatch))
+ (:method :before ((x standard-tensor) (y (eql nil)) &optional (axis 0))
+ (declare (ignore axis))
+ (assert (tensor-vectorp x) nil 'tensor-dimension-mismatch)))
+
+(defmethod sum! ((x standard-tensor) (y t) &optional (axis 0))
+ (let ((clx (class-name (class-of x))))
+ (assert (member clx *tensor-type-leaves*)
+ nil 'tensor-abstract-class :tensor-class (list clx))
+ (compile-and-eval
+ `(progn
+ (defmethod sum! ((x ,clx) (y ,clx) &optional (axis 0))
+ (t/sum ,clx x y axis))
+ (defmethod sum! ((x ,clx) (y (eql nil)) &optional (axis 0))
+ (declare (ignore axis))
+ (t/sum ,clx x nil))))
+ (sum! x y axis)))
+
+
+
commit 7aa696273c37819ebc5b9ee1040d0f194dd8145a
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Aug 19 20:19:10 2013 -0700
Removed the old loopy file.
diff --git a/src/old/loopy-old.lisp b/src/old/loopy-old.lisp
deleted file mode 100644
index 008dcd2..0000000
--- a/src/old/loopy-old.lisp
+++ /dev/null
@@ -1,118 +0,0 @@
-;;Very ugly inflexible code; get rid of this in some time or make use of mod-dotimes.
-(defmacro mod-loop ((idx dims) &body body)
- (check-type idx symbol)
- (let ((tensor-table (make-hash-table)))
- (labels ((get-tensors (decl)
- (if (null decl) t
- (let ((cdecl (car decl)))
- (when (and (eq (first cdecl) 'type)
- (get-tensor-class-optimization (second cdecl)))
- (dolist (sym (cddr cdecl))
- (let ((hsh (list
- :class (second cdecl)
- :stride-sym (gensym (string+ (symbol-name sym) "-stride"))
- :store-sym (gensym (string+ (symbol-name sym) "-store"))
- :offset-sym (gensym (string+ (symbol-name sym) "-offset"))
- :ref-count 0)))
- (setf (gethash sym tensor-table) hsh))))
- (get-tensors (cdr decl)))))
- (ttrans-p (code)
- (and (consp code) (eq (first code) 'tensor-ref)
- (gethash (second code) tensor-table)
- (eq (third code) idx)))
- (incref (ten)
- (incf (getf (gethash ten tensor-table) :ref-count)))
- (transform-setf-tensor-ref (snippet ret)
- (if (null snippet) ret
- (transform-setf-tensor-ref
- (cddr snippet)
- (append ret
- (destructuring-bind (to from &rest rest) snippet
- (declare (ignore rest))
- (let ((to-t? (ttrans-p to))
- (fr-t? (ttrans-p from)))
- (cond
- ((and to-t? fr-t?)
- (let ((to-opt (gethash (second to) tensor-table))
- (fr-opt (gethash (second from) tensor-table)))
- (if (eq (second (multiple-value-list (get-tensor-class-optimization (getf to-opt :class))))
- (second (multiple-value-list (get-tensor-class-optimization (getf fr-opt :class)))))
- (progn
- (incref (second to)) (incref (second from))
- (cdr (funcall (getf (get-tensor-class-optimization (getf to-opt :class)) :reader-writer)
- (getf fr-opt :store-sym) (getf fr-opt :offset-sym) (getf to-opt :store-sym) (getf to-opt :offset-sym))))
- (list to (find-tensor-refs from nil)))))
- (to-t?
- (incref (second to))
- (let ((to-opt (gethash (second to) tensor-table)))
- ;;Add type checking here!
- (cdr (funcall (getf (get-tensor-class-optimization (getf to-opt :class)) :value-writer)
- (find-tensor-refs from nil) (getf to-opt :store-sym) (getf to-opt :offset-sym)))))
- (t
- (list to (find-tensor-refs from nil))))))))))
- (transform-tensor-ref (snippet)
- (if (eq (first snippet) 'setf)
- (cons 'setf (transform-setf-tensor-ref (cdr snippet) nil))
- (destructuring-bind (tref ten index) snippet
- (assert (eq tref 'tensor-ref))
- (let ((topt (gethash ten tensor-table)))
- (if (not (and (eq index idx) topt)) snippet
- (progn
- (incref ten)
- (funcall (getf (get-tensor-class-optimization (getf topt :class)) :reader) (getf topt :store-sym) (getf topt :offset-sym))))))))
- (find-tensor-refs (code ret)
- (if (null code) (reverse ret)
- (cond
- ((consp code)
- (if (member (first code) '(tensor-ref setf))
- (transform-tensor-ref code)
- (find-tensor-refs (cdr code) (cons (find-tensor-refs (car code) nil) ret))))
- (t code)))))
- ;;
- (when (eq (caar body) 'declare)
- (get-tensors (cdar body)))
- (let ((tr-body (find-tensor-refs body nil)))
- (with-gensyms (dims-sym rank-sym count-sym)
- `(let* ((,dims-sym ,dims)
- (,rank-sym (length ,dims-sym))
- (,idx (allocate-index-store ,rank-sym))
- ,@(loop for key being the hash-keys of tensor-table
- when (> (getf (gethash key tensor-table) :ref-count) 0)
- collect (let ((hsh (gethash key tensor-table)))
- `(,(getf hsh :stride-sym) (strides ,key))))
- ,@(loop for key being the hash-keys of tensor-table
- when (> (getf (gethash key tensor-table) :ref-count) 0)
- collect (let ((hsh (gethash key tensor-table)))
- `(,(getf hsh :store-sym) (store ,key)))))
- (declare (type (index-array *) ,idx ,@(loop for key being the hash-keys of tensor-table
- when (> (getf (gethash key tensor-table) :ref-count) 0)
- collect (getf (gethash key tensor-table) :stride-sym)))
- ,@(loop for key being the hash-keys of tensor-table
- when (> (getf (gethash key tensor-table) :ref-count) 0)
- collect (let* ((hsh (gethash key tensor-table))
- (opt (get-tensor-class-optimization (getf hsh :class))))
- `(type ,(linear-array-type (getf opt :store-type)) ,(getf hsh :store-sym)))))
- (loop
- ,@(loop for key being the hash-keys of tensor-table
- when (> (getf (gethash key tensor-table) :ref-count) 0)
- append (let ((hsh (gethash key tensor-table)))
- `(with ,(getf hsh :offset-sym) of-type index-type = (head ,key))))
- do (locally
- ,@tr-body)
- ;;Optimized for row-order
- while (loop for ,count-sym of-type index-type from (1- ,rank-sym) downto 0
- do (if (= (aref ,idx ,count-sym) (1- (aref ,dims-sym ,count-sym)))
- (progn
- (setf (aref ,idx ,count-sym) 0)
- ,@(loop for key being the hash-keys of tensor-table
- when (> (getf (gethash key tensor-table) :ref-count) 0)
- collect (let ((hsh (gethash key tensor-table)))
- `(decf ,(getf hsh :offset-sym) (* (aref ,(getf hsh :stride-sym) ,count-sym) (1- (aref ,dims-sym ,count-sym)))))))
- (progn
- (incf (aref ,idx ,count-sym))
- ,@(loop for key being the hash-keys of tensor-table
- when (> (getf (gethash key tensor-table) :ref-count) 0)
- collect (let ((hsh (gethash key tensor-table)))
- `(incf ,(getf hsh :offset-sym) (aref ,(getf hsh :stride-sym) ,count-sym))))
- (return t)))
- finally (return nil)))))))))
commit 570ae7eb80324580ee27cfa7ba1d20b11f779e41
Author: Akshay Srinivasan <ak...@cs...>
Date: Mon Aug 19 15:48:54 2013 -0700
Minor fix to copy!
diff --git a/matlisp.asd b/matlisp.asd
index f5290e8..f3c236b 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -161,7 +161,8 @@
:pathname "lapack"
:depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
:components ((:file "lu")
- (:file "chol")))
+ (:file "chol")
+ (:file "eig")))
#+nil
(:module "matlisp-sugar"
:pathname "sugar"
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index 0d71f6d..f7432f7 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -157,7 +157,7 @@
,(recursive-append
(when (subtypep cly 'blas-numeric-tensor)
`(if-let (strd (and (call-fortran? y (t/l1-lb ,cly)) (consecutive-storep y)))
- (t/blas-copy! ,cly x nil y strd)))
+ (t/blas-copy! ,cly (t/coerce ,(field-type cly) x) nil y strd)))
`(t/copy! (t ,cly) x y))))
(copy! x y)))
commit 79a87b8605dba8ae97c8f354e42f1f081b127771
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Aug 19 09:52:54 2013 -0700
Not calling geev for workspace size, does not seem to save much time.
diff --git a/src/lapack/eig.lisp b/src/lapack/eig.lisp
index a907a7b..923aee7 100644
--- a/src/lapack/eig.lisp
+++ b/src/lapack/eig.lisp
@@ -52,6 +52,12 @@
(deft/generic (t/lapack-geev-workspace-inquiry #'subtypep) sym (n jobvl jobvr))
+#+nil
+(deft/method t/lapack-geev-workspace-inquiry (sym blas-numeric-tensor) (n jobvl jobvr)
+ (with-gensyms (n-sym)
+ `(let ((,n-sym ,n))
+ (* 10 ,n-sym))))
+
(deft/method t/lapack-geev-workspace-inquiry (sym blas-numeric-tensor) (n jobvl jobvr)
(using-gensyms (decl (n jobvl jobvr))
(with-gensyms (xxx)
commit fc524fd099c95abfc3af0280e8a200e461cc9493
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Aug 19 01:08:50 2013 -0700
Ported geev
diff --git a/src/base/template.lisp b/src/base/template.lisp
index 8ae4dc3..cd7daf8 100644
--- a/src/base/template.lisp
+++ b/src/base/template.lisp
@@ -39,8 +39,10 @@
(with-gensyms (num-sym)
`(let ((,num-sym ,num))
(cl:conjugate ,num-sym))))
+
(deft/method t/fc (ty real) (num)
num)
+
(defgeneric fc (x)
(:method ((x complex))
(cl:conjugate x))
@@ -55,6 +57,26 @@
(defun field-realp (fil)
(eql (macroexpand-1 `(t/fc ,fil phi)) 'phi))
;;
+(deft/generic (t/frealpart #'subtypep) ty (num))
+
+(deft/method t/frealpart (ty number) (num)
+ (with-gensyms (num-sym)
+ `(let ((,num-sym ,num))
+ (cl:realpart ,num-sym))))
+
+(deft/method t/frealpart (ty real) (num)
+ num)
+;;
+(deft/generic (t/fimagpart #'subtypep) ty (num))
+
+(deft/method t/fimagpart (ty number) (num)
+ (with-gensyms (num-sym)
+ `(let ((,num-sym ,num))
+ (cl:imagpart ,num-sym))))
+
+(deft/method t/fimagpart (ty real) (num)
+ `(t/fid+ ,ty))
+;;
(deft/generic (t/f= #'subtypep) ty (&rest nums))
(deft/method t/f= (ty number) (&rest nums)
(let* ((decl (zipsym nums))
diff --git a/src/foreign-core/lapack.lisp b/src/foreign-core/lapack.lisp
index 0806287..b587e50 100644
--- a/src/foreign-core/lapack.lisp
+++ b/src/foreign-core/lapack.lisp
@@ -717,7 +717,7 @@
(ldvr :integer :input)
(work (* :complex-double-float) :workspace-output)
(lwork :integer :input)
- (rwork (* :double-float) :workspace)
+ (rwork (* :double-float) :workspace-output)
(info :integer :output)
)
diff --git a/src/lapack/eig.lisp b/src/lapack/eig.lisp
index 308204a..a907a7b 100644
--- a/src/lapack/eig.lisp
+++ b/src/lapack/eig.lisp
@@ -50,6 +50,30 @@
(the index-type (head ,A)) (if ,vl (the index-type (head ,vl)) 0) (if ,vr (the index-type (head ,vr)) 0)))))
;;
+(deft/generic (t/lapack-geev-workspace-inquiry #'subtypep) sym (n jobvl jobvr))
+
+(deft/method t/lapack-geev-workspace-inquiry (sym blas-numeric-tensor) (n jobvl jobvr)
+ (using-gensyms (decl (n jobvl jobvr))
+ (with-gensyms (xxx)
+ `(let (,@decl
+ (,xxx (t/store-allocator ,sym 1)))
+ (declare (type index-type ,n)
+ (type character ,jobvl ,jobvr)
+ (type ,(store-type sym) ,xxx))
+ (,(macroexpand-1 `(t/lapack-geev-func ,sym))
+ ,jobvl ,jobvr
+ ,n
+ ,xxx ,n
+ ,xxx ,xxx
+ ,xxx ,n
+ ,xxx ,n
+ ,xxx -1
+ 0)
+ (ceiling (t/frealpart ,(field-type sym) (t/store-ref ,sym ,xxx 0)))))))
+
+#+nil
+(t/lapack-geev-workspace-inquiry complex-tensor 2 #\V #\V)
+
;;
#+nil
(progn
@@ -122,6 +146,49 @@
[3] INFO
where INFO is T if successful, NIL otherwise.
-"))
-
-
+")
+ (:method :before ((a standard-tensor) &optional vl vr)
+ (assert (tensor-squarep a) nil 'tensor-dimension-mismatch)
+ (when vl
+ (assert (and (tensor-squarep vl) (= (nrows vl) (nrows a)) (typep vl (type-of a))) nil 'tensor-dimension-mismatch))
+ (when vr
+ (assert (and (tensor-squarep vr) (= (nrows vr) (nrows a)) (typep vr (type-of a))) nil 'tensor-dimension-mismatch))))
+
+(defmethod geev! ((a blas-numeric-tensor) &optional vl vr)
+ (let ((cla (class-name (class-of A))))
+ (assert (member cla *tensor-type-leaves*)
+ nil 'tensor-abstract-class :tensor-class (list cla))
+ (compile-and-eval
+ `(defmethod geev! ((A ,cla) &optional vl vr)
+ (let* ((n (nrows A))
+ (jobvl (if vl #\V #\N))
+ (jobvr (if vr #\V #\N))
+ (work (t/store-allocator ,cla (t/lapack-geev-workspace-inquiry ,cla n jobvl jobvr)))
+ (wr (t/store-allocator ,cla n))
+ (wi (t/store-allocator ,cla n)))
+ (ecase jobvl
+ ,@(loop :for jvl :in '(#\N #\V)
+ :collect `(,jvl
+ (ecase jobvr
+ ,@(loop :for jvr :in '(#\N #\V)
+ :collect `(,jvr
+ (with-columnification (,cla () (A ,@(when (char= jvl #\V) `(vl)) ,@(when (char= jvr #\V) `(vr))))
+ (multiple-value-bind (osto owr owi ovl ovr owork info)
+ (t/lapack-geev! ,cla
+ A (or (blas-matrix-compatiblep A #\N) 0)
+ ,@(if (char= jvl #\N) `(nil 1) `(vl (or (blas-matrix-compatiblep vl #\N) 0)))
+ ,@(if (char= jvr #\N) `(nil 1) `(vr (or (blas-matrix-compatiblep vr #\N) 0)))
+ wr wi work)
+ (declare (ignore osto owr owi ovl ovr owork))
+ (unless (= info 0)
+ (error "geev returned ~a~%" info))))))))))
+ (values-list (remove-if #'null
+ (list
+ (let ((*check-after-initializing?* nil))
+ (make-instance 'complex-tensor ;',(complexified-type cla)
+ :dimensions (make-index-store (list (nrows A)))
+ :strides (make-index-store (list 1))
+ :head 0
+ :store (t/geev-output-fix ,cla wr wi)))
+ vl vr)))))))
+ (geev! A vl vr))
diff --git a/src/lapack/geev.lisp b/src/lapack/geev.lisp
index 531b004..9789585 100644
--- a/src/lapack/geev.lisp
+++ b/src/lapack/geev.lisp
@@ -372,11 +372,9 @@
(:nn (values "N" "N"))
((:vn t) (values "N" "V"))
(:nv (values "V" "N"))
- (:vv (values "V" "V")))
-
+ (:vv (values "V" "V")))
(let* ((ldvr (if (equal jobvr "V") n 1))
(ldvl (if (equal jobvl "V") n 1)))
-
(multiple-value-bind (store-a store-w store-vl store-vr work info)
(zgeev jobvl
jobvr
commit a3ea5898fd82c16b66bf0ce3d5615e370deb40a8
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Aug 18 23:32:53 2013 -0700
Saving state on geev porting.
diff --git a/src/base/template.lisp b/src/base/template.lisp
index f5d7fb3..8ae4dc3 100644
--- a/src/base/template.lisp
+++ b/src/base/template.lisp
@@ -86,6 +86,16 @@
(defun field-type (clname)
(macroexpand-1 `(t/field-type ,clname)))
+;;Hack? Yes.
+(defun complexified-type (ten)
+ (let ((ty (macroexpand-1 `(t/field-type ,ten))))
+ (if (subtypep ty 'complex) ten
+ (let* ((cty `(complex ,ty))
+ (table-entry (or (gethash 't/field-type matlisp-template::*template-table*) (ERROR "Undefined template : ~a~%" 'T/FIELD-TYPE))))
+ (car (find cty (mapcar #'(lambda (x) (list (cadr x) (funcall (car x) (cadr x))))
+ (getf table-entry :methods))
+ :key #'second :test #'list-eq))))))
+
;;Beware of infinite loops here.
(deft/generic (t/store-element-type #'subtypep) sym ())
(deft/method t/store-element-type (sym standard-tensor) ()
diff --git a/src/lapack/eig.lisp b/src/lapack/eig.lisp
new file mode 100644
index 0000000..308204a
--- /dev/null
+++ b/src/lapack/eig.lisp
@@ -0,0 +1,127 @@
+(in-package #:matlisp)
+
+(deft/generic (t/lapack-geev-func #'subtypep) sym ())
+
+(deft/method t/lapack-geev-func (sym real-tensor) ()
+ 'dgeev)
+;;Make API for real and complex versions similar.
+(definline mzgeev (jobvl jobvr n a lda w rwork vl ldvl vr ldvr work lwork info &optional (head-a 0) (head-vl 0) (head-vr 0))
+ (zgeev jobvl jobvr n a lda w vl ldvl vr ldvr work lwork rwork info head-a head-vl head-vr))
+(deft/method t/lapack-geev-func (sym complex-tensor) ()
+ 'mzgeev)
+;;
+(deft/generic (t/geev-output-fix #'subtypep) sym (wr wi))
+(deft/method t/geev-output-fix (sym real-numeric-tensor) (wr wi)
+ (let ((csym (or (complexified-type sym) (error "No corresponding complex-tensor defined for type ~a." sym))))
+ (using-gensyms (decl (wr wi))
+ (with-gensyms (ret i)
+ `(let* (,@decl
+ (,ret (t/store-allocator ,csym (length ,wr))))
+ (declare (type ,(store-type sym) ,wr ,wi)
+ (type ,(store-type csym) ,ret))
+ (very-quickly
+ (loop :for ,i :from 0 :below (length ,wr)
+ :do (t/store-set ,csym (complex (aref ,wr ,i) (aref ,wi ,i)) ,ret ,i)))
+ ,ret)))))
+
+(deft/method t/geev-output-fix (sym complex-numeric-tensor) (wr wi)
+ (using-gensyms (decl (wr))
+ `(let (,@decl)
+ (declare (type ,(store-type sym) ,wr))
+ ,wr)))
+;;
+(deft/generic (t/lapack-geev! #'subtypep) sym (A lda vl ldvl vr ldvr wr wi work))
+
+(deft/method t/lapack-geev! (sym blas-numeric-tensor) (A lda vl ldvl vr ldvr wr wi work)
+ (using-gensyms (decl (A lda vl ldvl vr ldvr wr wi work))
+ `(let (,@decl)
+ (declare (type ,sym ,A)
+ (type index-type ,lda)
+ (type ,(store-type sym) ,wr ,wi ,work))
+ (,(macroexpand-1 `(t/lapack-geev-func ,sym))
+ (if ,vl #\V #\N) (if ,vr #\V #\N)
+ (nrows ,A)
+ (the ,(store-type sym) (store ,A)) ,lda
+ ,wr ,wi
+ (if ,vl (the ,(store-type sym) (store ,vl)) (cffi:null-pointer)) (if ,vl ,ldvl 1)
+ (if ,vr (the ,(store-type sym) (store ,vr)) (cffi:null-pointer)) (if ,vr ,ldvr 1)
+ ,work (length ,work)
+ 0
+ (the index-type (head ,A)) (if ,vl (the index-type (head ,vl)) 0) (if ,vr (the index-type (head ,vr)) 0)))))
+;;
+
+;;
+#+nil
+(progn
+(let ((*default-tensor-type* 'complex-tensor))
+ (let ((a (copy! #2a((1 2) (3 4)) (zeros '(2 2)))))
+ (t/lapack-geev! complex-tensor a 2 (zeros '(2 2)) 2 (zeros '(2 2)) 2 (t/store-allocator complex-tensor 2) (t/store-allocator complex-tensor 2) (t/store-allocator complex-tensor 100))))
+
+
+(let ((a (copy! #2a((1 2) (3 4)) (zeros '(2 2)))))
+ (t/lapack-geev! real-tensor a 2 nil 1 nil 1 (t/store-allocator real-tensor 2) (t/store-allocator real-tensor 2) (t/store-allocator real-tensor 100)))
+)
+
+(defgeneric geev! (a &optional vl vr)
+ (:documentation "
+ Syntax
+ ======
+ (GEEV a [job])
+
+ Purpose:
+ ========
+ Computes the eigenvalues and left/right eigenvectors of A.
+
+ For an NxN matrix A, its eigenvalues are denoted by:
+
+ lambda(i), j = 1 ,..., N
+
+ The right eigenvectors of A are denoted by v(i) where:
+
+ A * v(i) = lambda(i) * v(i)
+
+ The left eigenvectors of A are denoted by u(i) where:
+
+ H H
+ u(i) * A = lambda(i) * u(i)
+
+ In matrix notation:
+ -1
+ A = V E V
+
+ and
+ -1
+ H H
+ A = U E U
+
+ where lambda(i) is the ith diagonal of the diagonal matrix E,
+ v(i) is the ith column of V and u(i) is the ith column of U.
+
+ The computed eigenvectors are normalized to have Euclidean norm
+ equal to 1 and largest component real.
+
+ Return Values:
+ ==============
+
+ JOB Return Values
+ ------------------------------------------------------------------
+ :NN (default) [1] (DIAG E) An Nx1 vector of eigenvalues
+ [2] INFO
+
+ :VN or T [1] V
+ [2] E
+ [3] INFO
+
+ :NV [1] E
+ [2] U
+ [3] INFO
+
+ :VV [1] V
+ [2] E
+ [3] U
+ [3] INFO
+
+ where INFO is T if successful, NIL otherwise.
+"))
+
+
diff --git a/src/lapack/geev.lisp b/src/lapack/geev.lisp
index 935d4a0..531b004 100644
--- a/src/lapack/geev.lisp
+++ b/src/lapack/geev.lisp
@@ -94,51 +94,6 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package #:matlisp)
-(deft/generic (t/lapack-geev-func #'subtypep) sym ())
-
-(deft/method t/lapack-geev-func (sym real-tensor) ()
- 'dgeev)
-;;Make API for real and complex versions similar.
-(definline mzgeev (jobvl jobvr n a lda w rwork vl ldvl vr ldvr work lwork info)
- (zgeev jobvl jobvr n a lda w vl ldvl vr ldvr work lwork rwork info))
-(deft/method t/lapack-geev-func (sym complex-tensor) ()
- 'mzgeev)
-;;
-
-(deft/generic (t/geev-output-fix #'subtypep) sym (wr wi))
-(deft/method t/geev-output-fix (sym real-numeric-tensor) (wr wi)
- (using-gensyms (decl (wr wi))
- (with-gensyms (ret)
- `(let* (,@decl
- (,ret (t/store-allocator ,sym (* 2 (length ,wr)))))
- (declare (type ,(store-type sym) ,wr ,wi ,ret))
- (very-quickly
- (loop :for i :from 0 :below (length ,wr)
- :do (setf (aref ,ret (* 2 i)) (aref ,wr i)
- (aref ,ret (1+ (* 2 i))) (aref ,wi i))))
- ,ret))))
-
-(deft/method t/geev-output-fix (sym complex-numeric-tensor) (wr wi)
- (using-gensyms (decl (wr))
- `(let (,@decl)
- (declare (type ,(store-type sym) ,wr))
- ,wr)))
-
-(deft/generic (t/lapack-geev! #'subtypep) sym (A lda vl ldvl vr ldvr w st-w))
-
-(deft/method t/lapack-potrf! (sym blas-numeric-tensor) (A lda uplo)
- (using-gensyms (decl (A lda uplo))
- `(let* (,@decl)
- (declare (type ,sym ,A)
- (type index-type ,lda)
- (type character ,uplo))
- (,(macroexpand-1 `(t/lapack-potrf-func ,sym))
- ,uplo
- (nrows ,A)
- (the ,(store-type sym) (store ,A)) ,lda
- 0))))
-
-
(defgeneric geev (a &optional job)
(:documentation
"
commit 69ca54a98c4e4a03e004268297644094b5541cae
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Aug 17 23:30:32 2013 -0700
Made changes to use (complex double-float) arrays in SBCL.
diff --git a/src/base/einstein.lisp b/src/base/einstein.lisp
index 64998a5..264aa92 100644
--- a/src/base/einstein.lisp
+++ b/src/base/einstein.lisp
@@ -121,7 +121,7 @@
(stp (when memp `(aref ,(car (getf plst :strides)) ,(position cidx (car ofst))))))
(get-incs (cdr idxs) (if memp (list `(the index-type (* ,(if tloop 1 stp) (aref ,(car (getf plst :dimensions)) ,(position cidx (car ofst)))))) nil)
(if (or tloop (and (null acc) (not memp))) (cons nil decl)
- (cons
+ (cons
(if memp
`(,dsym ,(if (null acc) stp `(the index-type (- ,stp ,@acc))) :type index-type)
`(,dsym (the index-type (- ,@acc)) :type index-type))
diff --git a/src/base/tweakable.lisp b/src/base/tweakable.lisp
index 6f8b0bb..86dca90 100644
--- a/src/base/tweakable.lisp
+++ b/src/base/tweakable.lisp
@@ -16,6 +16,8 @@
returns a 10x10 matrix with Column major order.
")
+(defparameter *default-tensor-type* 'real-tensor)
+
(defparameter *check-after-initializing?* t
"
If t, then check for invalid values in the field of
diff --git a/src/classes/numeric.lisp b/src/classes/numeric.lisp
index 93893f5..6f3e87e 100644
--- a/src/classes/numeric.lisp
+++ b/src/classes/numeric.lisp
@@ -53,39 +53,42 @@
(deft/method t/l3-lb (sym complex-numeric-tensor) ()
'*complex-l3-fcall-lb*)
-(deft/method t/store-element-type (sym complex-numeric-tensor) ()
- (let ((cplx-type (macroexpand-1 `(t/field-type ,sym))))
- (second cplx-type)))
-
-(deft/method t/compute-store-size (sym complex-numeric-tensor) (size)
- `(* 2 ,size))
-
-(deft/method t/store-ref (sym complex-numeric-tensor) (store idx)
- (let ((store-s (gensym))
- (idx-s (gensym))
- (type (macroexpand-1 `(t/store-element-type ,sym))))
- `(let ((,store-s ,store)
- (,idx-s ,idx))
- (declare (type (simple-array ,type) ,store-s))
- (complex (aref ,store-s (* 2 ,idx-s)) (aref ,store-s (1+ (* 2 ,idx-s)))))))
-
-(deft/method t/store-set (sym complex-numeric-tensor) (value store idx)
- (let ((store-s (gensym))
- (idx-s (gensym))
- (value-s (gensym))
- (type (macroexpand-1 `(t/store-element-type ,sym)))
- (ftype (macroexpand-1 `(t/field-type ,sym))))
- `(let ((,store-s ,store)
- (,idx-s ,idx)
- (,value-s ,value))
- (declare (type (simple-array ,type) ,store-s)
- (type ,ftype ,value-s))
- (setf (aref ,store-s (* 2 ,idx-s)) (cl:realpart ,value-s)
- (aref ,store-s (1+ (* 2 ,idx-s))) (cl:imagpart ,value-s))
- nil)))
-
-(defmethod store-size ((tensor complex-numeric-tensor))
- (floor (/ (length (store tensor)) 2)))
+;;SBCL uses specialized arrays for floating complex arrays.
+#-sbcl
+(progn
+ (deft/method t/store-element-type (sym complex-numeric-tensor) ()
+ (let ((cplx-type (macroexpand-1 `(t/field-type ,sym))))
+ (second cplx-type)))
+
+ (deft/method t/compute-store-size (sym complex-numeric-tensor) (size)
+ `(* 2 ,size))
+
+ (deft/method t/store-ref (sym complex-numeric-tensor) (store idx)
+ (let ((store-s (gensym))
+ (idx-s (gensym))
+ (type (macroexpand-1 `(t/store-element-type ,sym))))
+ `(let ((,store-s ,store)
+ (,idx-s ,idx))
+ (declare (type (simple-array ,type) ,store-s))
+ (complex (aref ,store-s (* 2 ,idx-s)) (aref ,store-s (1+ (* 2 ,idx-s)))))))
+
+ (deft/method t/store-set (sym complex-numeric-tensor) (value store idx)
+ (let ((store-s (gensym))
+ (idx-s (gensym))
+ (value-s (gensym))
+ (type (macroexpand-1 `(t/store-element-type ,sym)))
+ (ftype (macroexpand-1 `(t/field-type ,sym))))
+ `(let ((,store-s ,store)
+ (,idx-s ,idx)
+ (,value-s ,value))
+ (declare (type (simple-array ,type) ,store-s)
+ (type ,ftype ,value-s))
+ (setf (aref ,store-s (* 2 ,idx-s)) (cl:realpart ,value-s)
+ (aref ,store-s (1+ (* 2 ,idx-s))) (cl:imagpart ,value-s))
+ nil)))
+
+ (defmethod store-size ((tensor complex-numeric-tensor))
+ (floor (/ (length (store tensor)) 2))))
(defmethod print-element ((tensor complex-numeric-tensor)
element stream)
@@ -105,4 +108,3 @@
(defleaf scomplex-tensor (complex-numeric-tensor) ())
(deft/method t/store-element-type (sym scomplex-tensor) ()
'single-float))
-
diff --git a/src/ffi/ffi-cffi.lisp b/src/ffi/ffi-cffi.lisp
index 8ac6db3..5175304 100644
--- a/src/ffi/ffi-cffi.lisp
+++ b/src/ffi/ffi-cffi.lisp
@@ -146,6 +146,10 @@
(deftype matlisp-specialized-array ()
`(or (simple-array double-float (*))
(simple-array single-float (*))
+ #+sbcl
+ (simple-array (complex double-float) (*))
+ #+sbcl
+ (simple-array (complex single-float) (*))
;;
(simple-array (signed-byte 64) *)
(simple-array (signed-byte 32) *)
diff --git a/src/level-1/maker.lisp b/src/level-1/maker.lisp
index cf41512..f101dba 100644
--- a/src/level-1/maker.lisp
+++ b/src/level-1/maker.lisp
@@ -22,7 +22,7 @@
(t/zeros ,dtype dims)))
(zeros-generic dims dtype)))
-(definline zeros (dims &optional (type 'real-tensor))
+(definline zeros (dims &optional (type *default-tensor-type*))
(let ((*check-after-initializing?* nil))
(let ((type (etypecase type (standard-class (class-name type)) (symbol type))))
(etypecase dims
commit 6c6f96e88fab82f42a2cd563c53e90c48eb8da24
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Aug 17 17:34:10 2013 -0700
Removed redundant files.
diff --git a/src/classes/complex-tensor.lisp b/src/classes/complex-tensor.lisp
deleted file mode 100644
index 272285a..0000000
--- a/src/classes/complex-tensor.lisp
+++ /dev/null
@@ -1,44 +0,0 @@
-(in-package #:matlisp)
-
-(defmacro-method field-type ((sym (eql 'complex-tensor)))
- '(complex double-float))
-
-(defmacro-method store-element-type ((sym (eql 'complex-tensor)))
- 'double-float)
-
-(defmacro-method compute-store-size ((sym (eql 'complex-tensor)) size)
- '(* 2 size))
-
-(defmacro-method store-ref ((sym (eql 'complex-tensor)) store idx)
- (let ((store-s (gensym))
- (idx-s (gensym))
- (type (macroexpand-1 `(store-element-type ,sym))))
- `(let ((,store-s ,store)
- (,idx-s ,idx))
- (declare (type (simple-array ,type) ,store-s))
- (complex (aref ,store-s (* 2 ,idx-s...
[truncated message content] |
