Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch classy updated. 2012-02-23-280-g5304b72
|
From: Akshay S. <ak...@us...> - 2013-09-23 22:50:15
|
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 5304b7204035eab0b7ac2664a6e1949a0689e741 (commit)
via d112fcce019bbf7c536a4047927cfa248bff6239 (commit)
via 54c32278dd5b119ca1157b022ebe1a1b0f945f8e (commit)
via 06961b98935b57db5dd4d9b56bdd93c647ba6484 (commit)
via 22f5a0bcf4a70d769e8448e05840fc9ce8fe7988 (commit)
via e2cf244082f3b9993eb6d9e4f6051349f80ccbbd (commit)
via cd98eb7ed25c777623ccbacac627bb4968574536 (commit)
via 2247ca7cc973977e061ee894efee10faec823f1d (commit)
via e41ab636a047d01b438e86d24ad4b5169d0edfe2 (commit)
from f8b87a620796e228cadb86996b85f4298409ed75 (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 5304b7204035eab0b7ac2664a6e1949a0689e741
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Sep 23 15:41:49 2013 -0700
o Fixed bug in triangle template.
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index 3bcb166..3399780 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -145,7 +145,7 @@
:do (progn
,(if diag?
`(if (=
- (t/store-set ,sym (t/store-ref ,sym ,sto-a ,of-a) ,sto-b ,of-b))))))
+ (t/store-set ,sym (t/store-ref ,sym ,sto-a ,of-a) ,sto-b ,of-b))))))))))
,b))))
;;
commit d112fcce019bbf7c536a4047927cfa248bff6239
Author: Akshay Srinivasan <aks...@gm...>
Date: Thu Sep 19 22:39:13 2013 -0700
Added splot to gnuplot, added triangle-copy template.
diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp
index 59ba49b..d55f9da 100644
--- a/lib-src/gnuplot/gnuplot.lisp
+++ b/lib-src/gnuplot/gnuplot.lisp
@@ -1,11 +1,12 @@
(in-package :matlisp)
(defvar *current-gnuplot-process* nil)
-(defun open-gnuplot-stream (&optional (gnuplot-binary
- #+darwin
- (pathname "/opt/local/bin/gnuplot")
- #+linux
- (pathname "/usr/bin/gnuplot")))
+(defun open-gnuplot-stream (&key (gnuplot-binary
+ #+darwin
+ (pathname "/opt/local/bin/gnuplot")
+ #+linux
+ (pathname "/usr/bin/gnuplot"))
+ (terminal "wxt"))
(setf *current-gnuplot-process* (#+:sbcl
sb-ext:run-program
#+:ccl
@@ -13,7 +14,8 @@
gnuplot-binary nil :input :stream :wait nil :output t))
(gnuplot-send "
set datafile fortran
-")
+set term ~a
+" terminal)
*current-gnuplot-process*)
(defun close-gnuplot-stream ()
@@ -37,7 +39,7 @@ set datafile fortran
(multiple-value-bind (b2 b1) (floor a 256)
(list b2 b1 b0))))
-(defun plot2d (data &key (lines t) (color nil))
+(defun plot (data &key (lines t) (color nil))
(let ((fname "/tmp/matlisp-gnuplot.out"))
(with-open-file (s fname :direction :output :if-exists :supersede :if-does-not-exist :create)
(loop :for i :from 0 :below (loop :for x :in data :minimizing (size x))
@@ -45,7 +47,7 @@ set datafile fortran
(let ((col (if (listp color) color
(let ((lst (list color)))
(setf (cdr lst) lst)
- lst))))
+ lst))))
(let ((cmd (apply #'string+ (cons "plot " (loop :for x :in (cdr data)
:for i := 2 :then (1+ i)
:for clist := col :then (cdr clist)
@@ -60,6 +62,14 @@ set datafile fortran
(setf (aref cmd (- (length cmd) 2)) #\;
(aref cmd (- (length cmd) 1)) #\Newline)
(gnuplot-send cmd)))))
+
+(defun splot (data)
+ (let ((fname "/tmp/matlisp-gnuplot.out"))
+ (with-open-file (s fname :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 "~%"))))
+ (gnuplot-send (string+ "splot \'" fname "\'
+"))))
;; (defclass gnuplot-plot-info ()
;; ((title
diff --git a/packages.lisp b/packages.lisp
index ad00eb4..b0aa7ca 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -86,7 +86,7 @@
#:compile-and-eval
#:getcons #:mapcons
;;Macros
- #:when-let #:if-let #:if-ret #:with-gensyms #:let-rec #:using-gensyms
+ #:when-let #:if-let #:if-ret #:with-gensyms #:let-rec #:using-gensyms #:with-marking
#:mlet* #:make-array-allocator #:let-typed #:let*-typed
#:nconsc #:define-constant
#:macrofy #:looped-mapcar #:defun-compiler-macro
diff --git a/src/lapack/geqr.lisp b/src/lapack/geqr.lisp
index 57c56d5..2c52eba 100644
--- a/src/lapack/geqr.lisp
+++ b/src/lapack/geqr.lisp
@@ -47,42 +47,6 @@
")
(:method :before ((a standard-tensor))
(assert (tensor-matrixp a) nil 'tensor-dimension-mismatch)))
-
-(defmacro loop-lt ((dims-e &rest mats) &rest body)
- (let ((syms (mapcar #'(lambda (x)
- (let ((mat-sym (gensym)))
- `((,mat-sym ,(cadr x))
- (,(gensym "strd") (strides ,mat-sym))
- (,(car x) (head x)))))
- mats)))
- (with-gensyms (i j dims)
- `(let* (,@(apply #'append syms)
- (,dims ,dims-e))
- (with-marking
- (loop :for ,j :from 0 :below (mark (aref ,dims 1))
- :do (progn
- ,@(mapcar #'(lambda (x) `(incf ,(car (third x)) (mark (aref ,(car (second x)) 1)))))
- (loop :repeat :from 0 :below (mark (aref ,dims 0))
- :do (progn
- ,@body))))))))
-
-(deft/generic t/copy-upper-triangle (sym #'subtypep) (a b)
- (using-gensyms (decl (a b))
- (with-gensyms (sto-a sto-b strd-a strd-b)
- `(let (,@decl
- (,sto-a (store ,a))
- (,strd-a (strides ,a))
- (,sto-b (store ,b))
- (,strd-b (strides ,b)))
- (declare (type ,sym ,a ,b)
- (type ,(store-type sym) ,sto-a ,sto-b)
- (type index-store-vector ,strd-a ,strd-b))
- (very-quickly
- (loop :repeat (nrows ,a)
- :for rof-a :of-type index-type := (head a) :then (+ rof-a (aref strd-a 0))
- :for rof-a :of-type index-type := (head a) :then (+ rof-a (aref strd-a 0))
- :do (loop :repeat (ncols b)
- :do (t/store-set ,sym (t/store-ref ,sym sto-a ..) sto-b ..))))))))))
(defmethod geqr! ((a standard-tensor))
(let ((cla (class-name (class-of A))))
diff --git a/src/lapack/lu.lisp b/src/lapack/lu.lisp
index 4da2596..9a9be99 100644
--- a/src/lapack/lu.lisp
+++ b/src/lapack/lu.lisp
@@ -136,6 +136,9 @@
By default WITH-L,WITH-U,WITH-P.
"))
+(defmethod lu ((a standard-tensor) &optional split-lu?)
+ (let ((lu (getrf! (copy a))
+
#+nil
(defmacro make-lu (tensor-class)
(let* ((opt (if-ret (get-tensor-class-optimization-hashtable tensor-class)
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index f7432f7..3bcb166 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -117,6 +117,36 @@
(,of-y (strides ,y) (head ,y)))
:do (t/store-set ,cly ,cx ,sto-y ,of-y)))
,y))))
+;;
+;;This macro is used for interfacing with lapack
+;;Only to be used with matrices!
+(deft/generic (t/copy-triangle! #'subtypep) sym (a b &optional upper? diag?))
+(deft/method t/copy-triangle! (sym standard-tensor) (a b &optional (upper? nil) (diag? t))
+ (using-gensyms (decl (a b))
+ (with-gensyms (sto-a sto-b strd-a strd-b dof-a dof-b of-a of-b)
+ `(let* (,@decl
+ (,sto-a (store ,a))
+ (,strd-a (strides ,a))
+ (,sto-b (store ,b))
+ (,strd-b (strides ,b)))
+ (declare (type ,sym ,a ,b)
+ (type ,(store-type sym) ,sto-a ,sto-b)
+ (type index-store-vector ,strd-a ,strd-b))
+ (with-marking
+ (very-quickly
+ (:mark* ((ndiags (min (nrows ,a) (ncols ,a))))
+ (loop :for i :from 0 :below ndiags
+ :for ,dof-a :of-type index-type := (head ,a) :then (+ ,dof-a (:mark (lvec-foldr #'+ ,strd-a) :type index-type))
+ :for ,dof-b :of-type index-type := (head ,b) :then (+ ,dof-b (:mark (lvec-foldr #'+ ,strd-b) :type index-type))
+ :do (loop :for j :from 0 :below ,(if upper? `(1+ i) `(- ndiags i))
+ :for ,of-a :of-type index-type := ,dof-a :then (,(if upper? '- '+) ,of-a (:mark (aref ,strd-a 0)))
+ :for ,of-b :of-type index-type := ,dof-b :then (,(if upper? '- '+) ,of-b (:mark (aref ,strd-b 0)))
+ ,@(unless diag? `(:unless (= j 0)))
+ :do (progn
+ ,(if diag?
+ `(if (=
+ (t/store-set ,sym (t/store-ref ,sym ,sto-a ,of-a) ,sto-b ,of-b))))))
+ ,b))))
;;
(defmethod copy! :before ((x standard-tensor) (y standard-tensor))
diff --git a/src/special/map.lisp b/src/special/map.lisp
index 98b8443..9dcb66b 100644
--- a/src/special/map.lisp
+++ b/src/special/map.lisp
@@ -39,6 +39,10 @@
:do (t/store-set ,cly (funcall func (t/store-ref ,clx sto-x of-x)) sto-y of-y))))
y)))
(mapsor! func x y))
+
+(definline mapsor (func x)
+ (let ((ret (zeros (dimensions x) (class-of x))))
+ (mapsor! func x ret)))
;;
(defun mapslice (func x &optional (axis 0))
diff --git a/src/utilities/macros.lisp b/src/utilities/macros.lisp
index 199478c..c105360 100644
--- a/src/utilities/macros.lisp
+++ b/src/utilities/macros.lisp
@@ -40,7 +40,7 @@ Example:
(types nil)
(code (mapcons #'(lambda (mrk)
(ecase (car mrk)
- (mark*
+ (:mark*
`(symbol-macrolet (,@(mapcar #'(lambda (decl) (destructuring-bind (ref code &key type) decl
(let ((rsym (gensym (symbol-name ref))))
(push `(,rsym ,code) decls)
@@ -49,14 +49,14 @@ Example:
`(,ref ,rsym))))
(cadr mrk)))
,@(cddr mrk)))
- (mark
+ (:mark
(destructuring-bind (code &key type) (cdr mrk)
(let ((rsym (gensym)))
(push `(,rsym ,code) decls)
(when type
(push `(type ,type ,rsym) types))
rsym)))))
- body '(mark* mark))))
+ body '(:mark* :mark))))
`(let* (,@decls)
,@(when types `((declare ,@types)))
,@code)))
@@ -531,9 +531,9 @@ Example:
(defmacro slowly (&body forms)
"
Macro which encloses @arg{forms} inside
- (declare (optimize (speed 1)))
+ (declare (optimize (speed 1) (debug 3)))
"
- `(with-optimization (:speed 1)
+ `(with-optimization (:speed 1 :debug 3)
,@forms))
)
commit 54c32278dd5b119ca1157b022ebe1a1b0f945f8e
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Sep 17 15:37:01 2013 -0700
Fixed a bug in gemm.
Moved useful stuff from einstein.lisp into utilities.
Added "with-marking" macro for nicer local/global semantics.
diff --git a/packages.lisp b/packages.lisp
index f3d665d..ad00eb4 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -73,7 +73,7 @@
#:vectorify #:copy-n
#:ensure-args #:repsym #:findsym #:find-tag
#:zip #:zip-eq #:zipsym
- #:list-eq #:setadd #:setrem
+ #:list-eq #:setadd #:setrem #:set-eq
#:cut-cons-chain!
#:slot-values #:remmeth
#:recursive-append #:unquote-args #:flatten
@@ -84,6 +84,7 @@
#:lvec-map-foldl! #:lvec-map-foldr!
#:lvec->list #:lvec->list!
#:compile-and-eval
+ #:getcons #:mapcons
;;Macros
#:when-let #:if-let #:if-ret #:with-gensyms #:let-rec #:using-gensyms
#:mlet* #:make-array-allocator #:let-typed #:let*-typed
diff --git a/src/base/einstein.lisp b/src/base/einstein.lisp
index 264aa92..cf29c9b 100644
--- a/src/base/einstein.lisp
+++ b/src/base/einstein.lisp
@@ -1,23 +1,9 @@
(in-package :matlisp)
-(defun get-cons (lst sym)
- (if (atom lst) nil
- (if (eq (car lst) sym)
- (list lst)
- (append (get-cons (car lst) sym) (get-cons (cdr lst) sym)))))
-
(defun has-sym (lst sym)
(if (atom lst) (eql lst sym)
(or (has-sym (car lst) sym) (has-sym (cdr lst) sym))))
-(defun mapcons (func lst keys)
- (if (atom lst) lst
- (let ((tlst (if (member (car lst) keys)
- (funcall func lst)
- lst)))
- (if (atom tlst) tlst
- (mapcar #'(lambda (x) (mapcons func x keys)) tlst)))))
-
;;Only works for distinct objects
(defun generate-permutations (lst)
(if (null (cdr lst)) (list lst)
@@ -26,18 +12,8 @@
(mapcar #'(lambda (y) (cons x y)) (generate-permutations pop-x))))
lst))))
-(defun set-eq (a b &key (test #'eql))
- (and (loop :for ele :in a
- :do (unless (member ele b :test test)
- (return nil))
- :finally (return t))
- (loop :for ele :in b
- :do (unless (member ele a :test test)
- (return nil))
- :finally (return t))))
-
(defun parse-loopx (type place clause)
- (let* ((refs (let ((tmp (get-cons (list place clause) 'ref))
+ (let* ((refs (let ((tmp (getcons (list place clause) 'ref))
(ret nil))
(loop :for ele :in tmp
:do (setf ret (setadd ret ele #'equal)))
diff --git a/src/lapack/geqr.lisp b/src/lapack/geqr.lisp
index 0971de8..57c56d5 100644
--- a/src/lapack/geqr.lisp
+++ b/src/lapack/geqr.lisp
@@ -47,45 +47,6 @@
")
(:method :before ((a standard-tensor))
(assert (tensor-matrixp a) nil 'tensor-dimension-mismatch)))
-
-
-(defmacro with-marking (&rest body)
- (let* ((decls nil)
- (types nil)
- (code (mapcons #'(lambda (mrk)
- (ecase (car mrk)
- (mark*
- `(symbol-macrolet (,@(mapcar #'(lambda (decl) (destructuring-bind (ref code &key type) decl
- (let ((rsym (gensym (symbol-name ref))))
- (push `(,rsym ,code) decls)
- (when type
- (push `(type ,type ,rsym) types))
- `(,ref ,rsym))))
- (cadr mrk)))
- ,@(cddr mrk)))
- (mark
- (destructuring-bind (code &key type) (cdr mrk)
- (let ((rsym (gensym)))
- (push `(,rsym ,code) decls)
- (when type
- (push `(type ,type ,rsym) types))
- rsym)))))
- body '(mark* mark))))
- `(let* (,@decls)
- ,@(when types `((declare ,@types)))
- ,@code)))
-
-(with-marking
- (loop :for i := 0 :then (1+ i)
- :do (mark* ((xi (* 10 2) :type index-type)
- (sum 0 :type index-type))
- (incf sum (mark (* 10 2)))
- (if (= i 10)
- (return sum)))))
-
-(loop-upper-triangle ((dimensions x)
- (of-a a)
- (of-b b)))
(defmacro loop-lt ((dims-e &rest mats) &rest body)
(let ((syms (mapcar #'(lambda (x)
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index 3dfb285..75ba8a2 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -160,5 +160,5 @@
(defmethod gemv (alpha (A standard-tensor) (x standard-tensor)
(beta (eql nil)) (y (eql nil)) &optional (job :n))
- (let ((ret (zeros (nrows A) (class-of A))))
+ (let ((ret (zeros (ecase job (:n (nrows A)) (:t (ncols A))) (class-of A))))
(gemv! alpha A x 1 ret job)))
diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp
index b2cec50..7c824c0 100644
--- a/src/level-3/gemm.lisp
+++ b/src/level-3/gemm.lisp
@@ -140,7 +140,7 @@
(with-columnification (,cla ((a joba) (b jobb)) (c))
(multiple-value-bind (lda opa) (blas-matrix-compatiblep a joba)
(multiple-value-bind (ldb opb) (blas-matrix-compatiblep b jobb)
- (t/blas-gemm! ,cla alpha A lda B ldb beta C ldb opa opb))))))
+ (t/blas-gemm! ,cla alpha A lda B ldb beta C (or (blas-matrix-compatiblep c #\N) 0) opa opb))))))
`(t/gemm! ,cla alpha A B beta C joba jobb))))
C))
(gemm! alpha A B beta C job))
@@ -175,3 +175,11 @@
(defmethod gemm (alpha (A standard-tensor) (B standard-tensor)
beta (C standard-tensor) &optional (job :nn))
(gemm! alpha A B beta (copy C) job))
+
+(defmethod gemm (alpha (A standard-tensor) (B standard-tensor)
+ (beta (eql nil)) (C (eql nil)) &optional (job :nn))
+ (let ((ret (destructuring-bind (job-a job-b) (split-job job)
+ (zeros (list (ecase job-a (#\N (nrows A)) ((#\C #\T) (ncols A)))
+ (ecase job-b (#\N (ncols B)) ((#\C #\T) (nrows B))))
+ (class-of A)))))
+ (gemm! alpha A B 1 ret job)))
diff --git a/src/utilities/functions.lisp b/src/utilities/functions.lisp
index f790208..6494995 100644
--- a/src/utilities/functions.lisp
+++ b/src/utilities/functions.lisp
@@ -56,6 +56,16 @@
(cdr lst)
(cons (car lst) (setrem (cdr lst) a test)))))
+(defun set-eq (a b &key (test #'eql))
+ (and (loop :for ele :in a
+ :do (unless (member ele b :test test)
+ (return nil))
+ :finally (return t))
+ (loop :for ele :in b
+ :do (unless (member ele a :test test)
+ (return nil))
+ :finally (return t))))
+
(declaim (inline copy-n))
(defun copy-n (vec lst n)
(declare (type vector vec)
@@ -66,6 +76,20 @@
:do (setf (car vlst) (aref vec i)))
lst)
+(defun getcons (lst sym)
+ (if (atom lst) nil
+ (if (eq (car lst) sym)
+ (list lst)
+ (append (getcons (car lst) sym) (getcons (cdr lst) sym)))))
+
+(defun mapcons (func lst keys)
+ (if (atom lst) lst
+ (let ((tlst (if (member (car lst) keys)
+ (funcall func lst)
+ lst)))
+ (if (atom tlst) tlst
+ (mapcar #'(lambda (x) (mapcons func x keys)) tlst)))))
+
(defun find-tag (lst tag)
(let ((car (car lst)))
(if (atom car)
diff --git a/src/utilities/macros.lisp b/src/utilities/macros.lisp
index f19868f..199478c 100644
--- a/src/utilities/macros.lisp
+++ b/src/utilities/macros.lisp
@@ -10,6 +10,57 @@
`(defconstant ,name (if (boundp ',name) (symbol-value ',name) ,value)
,@(when doc (list doc))))
+(defmacro with-marking (&rest body)
+ "
+ This macro basically declares local-variables globally,
+ while keeping semantics and scope local.
+
+Example:
+ > (macroexpand-1
+ `(with-marking
+ (loop :for i := 0 :then (1+ i)
+ :do (mark* ((xi (* 10 2) :type index-type)
+ (sum 0 :type index-type))
+ (incf sum (mark (* 10 2)))
+ (if (= i 10)
+ (return sum))))))
+
+ (LET* ((#:G1083 (* 10 2)) (#:SUM1082 0) (#:XI1081 (* 10 2)))
+ (DECLARE (TYPE INDEX-TYPE #:SUM1082)
+ (TYPE INDEX-TYPE #:XI1081))
+ (LOOP :FOR I := 0 :THEN (1+ I)
+ :DO (SYMBOL-MACROLET ((XI #:XI1081) (SUM #:SUM1082))
+ (INCF SUM #:G1083)
+ (IF (= I 10)
+ (RETURN SUM)))))
+ T
+ >
+"
+ (let* ((decls nil)
+ (types nil)
+ (code (mapcons #'(lambda (mrk)
+ (ecase (car mrk)
+ (mark*
+ `(symbol-macrolet (,@(mapcar #'(lambda (decl) (destructuring-bind (ref code &key type) decl
+ (let ((rsym (gensym (symbol-name ref))))
+ (push `(,rsym ,code) decls)
+ (when type
+ (push `(type ,type ,rsym) types))
+ `(,ref ,rsym))))
+ (cadr mrk)))
+ ,@(cddr mrk)))
+ (mark
+ (destructuring-bind (code &key type) (cdr mrk)
+ (let ((rsym (gensym)))
+ (push `(,rsym ,code) decls)
+ (when type
+ (push `(type ,type ,rsym) types))
+ rsym)))))
+ body '(mark* mark))))
+ `(let* (,@decls)
+ ,@(when types `((declare ,@types)))
+ ,@code)))
+
(defmacro mlet* (vars &rest body)
"
This macro extends the syntax of let* to handle multiple values, it also handles
commit 06961b98935b57db5dd4d9b56bdd93c647ba6484
Merge: 22f5a0b e2cf244
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Sep 17 14:57:54 2013 -0700
Merge branch 'classy' of bicycle.cs.washington.edu:/homes/gws/akshays/git/matlisp into classy
commit 22f5a0bcf4a70d769e8448e05840fc9ce8fe7988
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Sep 17 14:56:33 2013 -0700
Added default gnuplot path for darwin.
diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp
index f9b48dc..59ba49b 100644
--- a/lib-src/gnuplot/gnuplot.lisp
+++ b/lib-src/gnuplot/gnuplot.lisp
@@ -1,7 +1,11 @@
(in-package :matlisp)
(defvar *current-gnuplot-process* nil)
-(defun open-gnuplot-stream (&optional (gnuplot-binary (pathname "/usr/bin/gnuplot")))
+(defun open-gnuplot-stream (&optional (gnuplot-binary
+ #+darwin
+ (pathname "/opt/local/bin/gnuplot")
+ #+linux
+ (pathname "/usr/bin/gnuplot")))
(setf *current-gnuplot-process* (#+:sbcl
sb-ext:run-program
#+:ccl
commit e2cf244082f3b9993eb6d9e4f6051349f80ccbbd
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Sep 17 14:33:14 2013 -0700
Saving state on QR.
diff --git a/src/lapack/geqr.lisp b/src/lapack/geqr.lisp
index e1af058..0971de8 100644
--- a/src/lapack/geqr.lisp
+++ b/src/lapack/geqr.lisp
@@ -48,21 +48,62 @@
(:method :before ((a standard-tensor))
(assert (tensor-matrixp a) nil 'tensor-dimension-mismatch)))
-(defmacro loop-upper-triangle ((dims-e &rest mats) &rest body)
+
+(defmacro with-marking (&rest body)
+ (let* ((decls nil)
+ (types nil)
+ (code (mapcons #'(lambda (mrk)
+ (ecase (car mrk)
+ (mark*
+ `(symbol-macrolet (,@(mapcar #'(lambda (decl) (destructuring-bind (ref code &key type) decl
+ (let ((rsym (gensym (symbol-name ref))))
+ (push `(,rsym ,code) decls)
+ (when type
+ (push `(type ,type ,rsym) types))
+ `(,ref ,rsym))))
+ (cadr mrk)))
+ ,@(cddr mrk)))
+ (mark
+ (destructuring-bind (code &key type) (cdr mrk)
+ (let ((rsym (gensym)))
+ (push `(,rsym ,code) decls)
+ (when type
+ (push `(type ,type ,rsym) types))
+ rsym)))))
+ body '(mark* mark))))
+ `(let* (,@decls)
+ ,@(when types `((declare ,@types)))
+ ,@code)))
+
+(with-marking
+ (loop :for i := 0 :then (1+ i)
+ :do (mark* ((xi (* 10 2) :type index-type)
+ (sum 0 :type index-type))
+ (incf sum (mark (* 10 2)))
+ (if (= i 10)
+ (return sum)))))
+
+(loop-upper-triangle ((dimensions x)
+ (of-a a)
+ (of-b b)))
+
+(defmacro loop-lt ((dims-e &rest mats) &rest body)
(let ((syms (mapcar #'(lambda (x)
- (let ((mat-sym (gensyms)))
- `((,mat-sym ,x)
- (,(gensym "sto") (store ,mat-sym))
- (,(gensym "strides") (strides ,mat-sym))
- (,(gensym "dimensions") (dimensions ,mat-sym)))))
+ (let ((mat-sym (gensym)))
+ `((,mat-sym ,(cadr x))
+ (,(gensym "strd") (strides ,mat-sym))
+ (,(car x) (head x)))))
mats)))
(with-gensyms (i j dims)
- `(let (,@(apply #'append syms)
- (,dims ,dims-e))
- (loop :for ,i :from 0 :below (aref ,dims 0)
- :do (loop :for ,j :from 0 :below (aref ,dims 1)
- :do (progn
- ,@body)))))))
+ `(let* (,@(apply #'append syms)
+ (,dims ,dims-e))
+ (with-marking
+ (loop :for ,j :from 0 :below (mark (aref ,dims 1))
+ :do (progn
+ ,@(mapcar #'(lambda (x) `(incf ,(car (third x)) (mark (aref ,(car (second x)) 1)))))
+ (loop :repeat :from 0 :below (mark (aref ,dims 0))
+ :do (progn
+ ,@body))))))))
(deft/generic t/copy-upper-triangle (sym #'subtypep) (a b)
(using-gensyms (decl (a b))
diff --git a/tests/loopy-tests.lisp b/tests/loopy-tests.lisp
index cceb145..1f8a22e 100644
--- a/tests/loopy-tests.lisp
+++ b/tests/loopy-tests.lisp
@@ -1,15 +1,20 @@
(in-package :matlisp)
-(defun tdcopy (n)
- (let* ((t-a (make-real-tensor-dims n n n))
+
+(defun tcopy (n &optional (rank 2))
+ (let* ((dims (make-list rank :initial-element n))
+ (t-a (zeros dims))
(st-a (store t-a))
- (t-b (make-real-tensor-dims n n n))
+ (t-b (zeros dims))
(st-b (store t-b)))
- (with-optimization (:speed 3 :safety 0 :space 0)
- (mod-dotimes (idx (idxv n n))
- with (linear-sums
- (of (idxv (* n n) n)))
- do (dcopy n st-a 1 st-b 1 of of)))))
+ (declare (type (simple-array double-float (*)) st-a st-b))
+ (time
+ (very-quickly
+ (mod-dotimes (idx (dimensions t-a))
+ :with (linear-sums
+ (of-a (strides t-a) (head t-a))
+ (of-b (strides t-b) (head t-b)))
+ :do (setf (aref st-b of-b) (aref st-a of-a)))))))
(defun tcopy (n)
(let* ((t-a (make-real-tensor-dims n n n))
commit cd98eb7ed25c777623ccbacac627bb4968574536
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Sep 16 12:05:07 2013 -0700
Saving state.
diff --git a/src/lapack/geqr.lisp b/src/lapack/geqr.lisp
index 17c1818..e1af058 100644
--- a/src/lapack/geqr.lisp
+++ b/src/lapack/geqr.lisp
@@ -48,8 +48,22 @@
(:method :before ((a standard-tensor))
(assert (tensor-matrixp a) nil 'tensor-dimension-mismatch)))
-
-(defmacro loop-upper-triangle (cla
+(defmacro loop-upper-triangle ((dims-e &rest mats) &rest body)
+ (let ((syms (mapcar #'(lambda (x)
+ (let ((mat-sym (gensyms)))
+ `((,mat-sym ,x)
+ (,(gensym "sto") (store ,mat-sym))
+ (,(gensym "strides") (strides ,mat-sym))
+ (,(gensym "dimensions") (dimensions ,mat-sym)))))
+ mats)))
+ (with-gensyms (i j dims)
+ `(let (,@(apply #'append syms)
+ (,dims ,dims-e))
+ (loop :for ,i :from 0 :below (aref ,dims 0)
+ :do (loop :for ,j :from 0 :below (aref ,dims 1)
+ :do (progn
+ ,@body)))))))
+
(deft/generic t/copy-upper-triangle (sym #'subtypep) (a b)
(using-gensyms (decl (a b))
(with-gensyms (sto-a sto-b strd-a strd-b)
commit 2247ca7cc973977e061ee894efee10faec823f1d
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Sep 15 11:55:55 2013 -0700
Saving working state on geqr.
diff --git a/src/foreign-core/lapack.lisp b/src/foreign-core/lapack.lisp
index b587e50..1c94e74 100644
--- a/src/foreign-core/lapack.lisp
+++ b/src/foreign-core/lapack.lisp
@@ -1068,7 +1068,7 @@
(m :integer :input)
(n :integer :input)
(k :integer :input)
- (a (* :complex-double-float) :input-output)
+ (a (* :complex-double-float :inc head-a) :input-output)
(lda :integer :input)
(tau (* :complex-double-float) :input)
(work (* :complex-double-float) :workspace-output)
@@ -1236,7 +1236,7 @@
(m :integer :input)
(n :integer :input)
(k :integer :input)
- (a (* :double-float) :input-output)
+ (a (* :double-float :inc head-a) :input-output)
(lda :integer :input)
(tau (* :double-float) :input)
(work (* :double-float) :workspace-output)
diff --git a/src/lapack/geqr.lisp b/src/lapack/geqr.lisp
index 3d129a1..17c1818 100644
--- a/src/lapack/geqr.lisp
+++ b/src/lapack/geqr.lisp
@@ -1,43 +1,26 @@
-;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
-;;;
-;;; $Id: geqr.lisp,v 1.7 2002/01/20 00:42:25 simsek Exp $
-;;;
-;;; $Log: geqr.lisp,v $
-;;; Revision 1.7 2002/01/20 00:42:25 simsek
-;;; o removed a spurious ignore
-;;;
-;;; Revision 1.6 2002/01/08 19:40:45 rtoy
-;;; The functions we use are exported now.
-;;;
-;;; Revision 1.5 2001/10/29 18:00:28 rtoy
-;;; Updates from M. Koerber to support QR routines with column pivoting:
-;;;
-;;; o Add an integer4 type and allocate-integer4-store routine.
-;;; o Add the necessary Fortran routines
-;;; o Add Lisp interface to the Fortran routines
-;;; o Update geqr for the new routines.
-;;;
-;;; Revision 1.4 2001/10/29 17:34:34 rtoy
-;;; I (RLT) stupidly deleted too much from M. Koerber's update. This is
-;;; his latest version.
-;;;
-;;; Revision 1.3 2001/10/26 15:19:25 rtoy
-;;; Renamed optional SKINNY parameter to ECON.
-;;;
-;;; Revision 1.2 2001/10/26 13:37:03 rtoy
-;;; Correctly handle the case when rows > cols and we want the [q1 q2]
-;;; form. Fix from M. Koerber.
-;;;
-;;; Revision 1.1 2001/10/25 21:51:58 rtoy
-;;; Initial revision for QR routines.
-;;;
-
(in-package #:matlisp)
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;; Set up the methods required to handle general matricies of Real
-;; and complex types. There are numerous other special cases, but
-;; they will not be considered for this first release. mak
+(deft/generic (t/lapack-geqrf-func #'subtypep) sym ())
+(deft/method t/lapack-geqrf-func (sym real-tensor) ()
+ 'matlisp-lapack:dgeqrf)
+(deft/method t/lapack-geqrf-func (sym complex-tensor) ()
+ 'matlisp-lapack:zgeqrf)
+;;
+(deft/generic (t/lapack-geqrf-workspace-inquiry #'subtypep) sym (m n))
+(deft/method t/lapack-geqrf-workspace-inquiry (sym blas-numeric-tensor) (m n)
+ (using-gensyms (decl (m n))
+ (with-gensyms (xxx)
+ `(let (,@decl
+ (,xxx (t/store-allocator ,sym 1)))
+ (declare (type index-type ,m ,n)
+ (type ,(store-type sym) ,xxx))
+ (,(macroexpand-1 `(t/lapack-geqrf-func ,sym))
+ ,m ,n
+ ,xxx ,m
+ ,xxx ,xxx -1 0)
+ (ceiling (t/frealpart ,(field-type sym) (t/store-ref ,sym ,xxx 0)))))))
+
+;;
(defgeneric geqr! (a)
(:documentation
"
@@ -61,68 +44,59 @@
[2] R
If the factorization can not be done, Q and R are set to NIL.
-
- NOTE: THIS FUNCTION IS DESTRUCTIVE.
-"))
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;; One could simply use LWORK = (MAX 1 N), but this call might result
-;; in some optimization in performance. For small matricies this is
-;; probably a 'no-net-gain' operation...but I seldom use small matricies
-;; in my work ;-) ... mak
-(let ((xx (allocate-real-store 1))
- (work (allocate-real-store 1)))
-
- (defun dgeqrf-workspace-inquiry (m n)
- (multiple-value-bind (store-a store-tau store-work lwork info)
- (lapack:dgeqrf m n xx m xx work -1 0)
-
- (declare (ignore store-a store-tau store-work lwork info))
-
- (values (ceiling (realpart (aref work 0)))))))
-
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-(let ((xx (allocate-complex-store 1))
- (work (allocate-complex-store 1)))
-
- (defun zgeqrf-workspace-inquiry (m n)
-
- (multiple-value-bind (store-a store-tau store-work lwork info)
- (lapack:zgeqrf m n xx m xx work -1 0)
-
- (declare (ignore store-a store-tau store-work lwork info))
-
- (values (ceiling (realpart (aref work 0)))))))
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;; Okay...now we build up the specific method for real and comples
-(defmethod geqr! ((a real-matrix))
-
- (let* ((m (nrows a))
- (n (ncols a))
- (k (min m n)) ; THESE ROUTINES ONLY RETURN A MINIMUM Q!
- (tau (allocate-real-store k)) ; reflection factors
- (lwork (dgeqrf-workspace-inquiry m n)) ; optimum work array size
- (work (allocate-real-store lwork))) ; and the work area
-
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ;; Do the Householder portion of the decomposition
- (multiple-value-bind (q-r new-tau new-work info)
- (lapack:dgeqrf m n (store a) m tau work lwork 0)
-
- (declare (ignore new-work))
- ;; Q-R and NEW-TAU aren't needed either since the (STORE A) and WORK
- ;; get modified
-
- (if (not (zerop info))
- ;; If INFO is not zero, then an error occured. Return Nil
- ;; for the Q and R and print a warning
- (progn (warn "QR Decomp failed: Argument ~d in call to DGEQRF is bad" (- info))
- (values nil nil))
-
- ;; If we are here, then INFO == 0 and all is well...
- (let ((r (make-real-matrix k n)))
+")
+ (:method :before ((a standard-tensor))
+ (assert (tensor-matrixp a) nil 'tensor-dimension-mismatch)))
+
+
+(defmacro loop-upper-triangle (cla
+(deft/generic t/copy-upper-triangle (sym #'subtypep) (a b)
+ (using-gensyms (decl (a b))
+ (with-gensyms (sto-a sto-b strd-a strd-b)
+ `(let (,@decl
+ (,sto-a (store ,a))
+ (,strd-a (strides ,a))
+ (,sto-b (store ,b))
+ (,strd-b (strides ,b)))
+ (declare (type ,sym ,a ,b)
+ (type ,(store-type sym) ,sto-a ,sto-b)
+ (type index-store-vector ,strd-a ,strd-b))
+ (very-quickly
+ (loop :repeat (nrows ,a)
+ :for rof-a :of-type index-type := (head a) :then (+ rof-a (aref strd-a 0))
+ :for rof-a :of-type index-type := (head a) :then (+ rof-a (aref strd-a 0))
+ :do (loop :repeat (ncols b)
+ :do (t/store-set ,sym (t/store-ref ,sym sto-a ..) sto-b ..))))))))))
+
+(defmethod geqr! ((a standard-tensor))
+ (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 geqr! ((a ,cla))
+ (let* ((m (nrows a))
+ (n (ncols a))
+ (k (min m n)) ; THESE ROUTINES ONLY RETURN A MINIMUM Q!
+ (tau (t/store-allocator ,cla k)) ; reflection factors
+ (lwork (t/lapack-geqrf-workspace-inquiry m n)) ; optimum work array size
+ (work (t/store-allocator ,cla lwork))) ; and the work area
+ (declare (type index-type lwork m n k)
+ (type ,(store-type cla) tau work))
+ ;; Do the Householder portion of the decomposition
+ (with-columnification (,cla () (A))
+ (multiple-value-bind (q-r new-tau new-work info)
+ (,(macroexpand-1 `(t/lapack-geqrf-func ,cla))
+ m n
+ (the ,(store-type cla) (store A)) (or (blas-matrix-compatiblep A #\N) 0)
+ tau work lwork 0 (the index-type (head A)))
+ (declare (ignore q-r new-tau new-work))
+ (unless (= info 0)
+ (error "geqrf returned ~a~%" info))
+
+
+
+ ;; If we are here, then INFO == 0 and all is well...
+ (let ((r (make-real-matrix k n)))
;; Extract the matrix R from Q-R
(dotimes (row k)
(loop for col from row below n do
commit e41ab636a047d01b438e86d24ad4b5169d0edfe2
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Sep 14 12:03:29 2013 -0700
Made real-version of eig more generic.
diff --git a/src/lapack/eig.lisp b/src/lapack/eig.lisp
index c603342..f0a92da 100644
--- a/src/lapack/eig.lisp
+++ b/src/lapack/eig.lisp
@@ -199,8 +199,6 @@
vl vr)))))))
(geev! A vl vr))
;;
-
-
(defgeneric eig (matrix &optional job)
(:method :before ((matrix standard-tensor) &optional (job :nn))
(assert (tensor-matrixp matrix) nil 'tensor-dimension-mismatch)
@@ -213,7 +211,6 @@
(vr (when revec? (zeros (list n n) (class-of matrix)))))
(geev! (copy matrix) vl vr)))
-
(defun geev-fix-up-eigvec (n eigval eigvec)
(let* ((evec (copy! eigvec (zeros (list n n) (complexified-type (class-of eigvec)))))
(tmp (zeros n (complexified-type (class-of eigvec))))
@@ -236,19 +233,17 @@
(return nil)))
evec))
-(defmethod eig ((matrix real-tensor) &optional (job :nn))
+(defmethod eig ((matrix real-numeric-tensor) &optional (job :nn))
(mlet* ((n (nrows matrix))
((levec? revec?) (values-list (mapcar #'(lambda (x) (char= x #\V)) (split-job job))))
(ret (multiple-value-list
(geev! (copy matrix)
- (when levec? (zeros (list n n) 'real-tensor))
- (when revec? (zeros (list n n) 'real-tensor)))))
+ (when levec? (zeros (list n n) (class-of matrix)))
+ (when revec? (zeros (list n n) (class-of matrix))))))
(eig (car ret)))
- (if (let ((stoe (store eig)))
- (loop :for i :from 0 :below n
- :do (unless (zerop (t/fc (t/field-type complex-tensor) (t/store-ref complex-tensor stoe i)))
- (return nil))
- :finally (return t)))
- (values-list ret)
- (values-list (cons eig (mapcar #'(lambda (mat) (geev-fix-up-eigvec n eig mat)) (cdr ret)))))))
-
+ (if (loop :for i :from 0 :below n
+ :do (unless (zerop (imagpart (ref eig i)))
+ (return nil))
+ :finally (return t))
+ (values-list ret)
+ (values-list (cons eig (mapcar #'(lambda (mat) (geev-fix-up-eigvec n eig mat)) (cdr ret)))))))
diff --git a/src/lapack/geev.lisp b/src/lapack/geev.lisp
deleted file mode 100644
index 9789585..0000000
--- a/src/lapack/geev.lisp
+++ /dev/null
@@ -1,524 +0,0 @@
-;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Copyright (c) 2000 The Regents of the University of California.
-;;; All rights reserved.
-;;;
-;;; Permission is hereby granted, without written agreement and without
-;;; license or royalty fees, to use, copy, modify, and distribute this
-;;; software and its documentation for any purpose, provided that the
-;;; above copyright notice and the following two paragraphs appear in all
-;;; copies of this software.
-;;;
-;;; IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
-;;; FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
-;;; ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
-;;; THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
-;;; SUCH DAMAGE.
-;;;
-;;; THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
-;;; INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-;;; MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
-;;; PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
-;;; CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
-;;; ENHANCEMENTS, OR MODIFICATIONS.
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Originally written by Raymond Toy
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; $Id: geev.lisp,v 1.11 2005/08/20 13:50:27 rtoy Exp $
-;;;
-;;; $Log: geev.lisp,v $
-;;; Revision 1.11 2005/08/20 13:50:27 rtoy
-;;; Fix problem with geev-workspace-inquiry functions when given a job
-;;; of :vv. This is a slightly modified version of the patch by Paul
-;;; Ledbetter III, on matlisp-user, 2005-08-16.
-;;;
-;;; We also removed the extra xxx array and used work instead.
-;;;
-;;; Revision 1.10 2001/10/26 15:24:16 rtoy
-;;; From M. Koerber:
-;;;
-;;; When determining LWORK, if JOBVR is V, the LDVR must be >= N.
-;;;
-;;; Revision 1.9 2001/06/22 12:52:41 rtoy
-;;; Use ALLOCATE-REAL-STORE and ALLOCATE-COMPLEX-STORE to allocate space
-;;; instead of using the error-prone make-array.
-;;;
-;;; Revision 1.8 2001/03/07 00:14:56 rtoy
-;;; Asking dgeev for the desired workspace size now works. (Didn't have
-;;; the return values matched up correctly!) (Needs a fix to dgeev.f for
-;;; this to work, though.)
-;;;
-;;; Revision 1.7 2001/03/06 21:58:06 rtoy
-;;; o The workspace inquiry function doesn't seem to work for DGEEV.
-;;; Don't use it in geev.
-;;; o The workspace was too small when inquiring the workspace for ZGEEV.
-;;; Make it larger.
-;;;
-;;; Revision 1.6 2001/02/23 14:04:20 rtoy
-;;; The Fortran geev routines allow the user to inquire about the optimum
-;;; size of the work array. Use that to allocate the appropriate amount
-;;; of space.
-;;;
-;;; Revision 1.5 2001/02/23 13:13:56 rtoy
-;;; The length of the work array was half-sized! (Despite the name,
-;;; complex-matrix-element-type is not a complex number. It's just a real
-;;; number)
-;;;
-;;; Revision 1.4 2000/07/11 18:02:03 simsek
-;;; o Added credits
-;;;
-;;; Revision 1.3 2000/07/11 02:11:56 simsek
-;;; o Added support for Allegro CL
-;;;
-;;; Revision 1.2 2000/05/08 17:19:18 rtoy
-;;; Changes to the STANDARD-MATRIX class:
-;;; o The slots N, M, and NXM have changed names.
-;;; o The accessors of these slots have changed:
-;;; NROWS, NCOLS, NUMBER-OF-ELEMENTS
-;;; The old names aren't available anymore.
-;;; o The initargs of these slots have changed:
-;;; :nrows, :ncols, :nels
-;;;
-;;; Revision 1.1 2000/04/14 00:11:12 simsek
-;;; o This file is adapted from obsolete files 'matrix-float.lisp'
-;;; 'matrix-complex.lisp' and 'matrix-extra.lisp'
-;;; o Initial revision.
-;;;
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-(in-package #:matlisp)
-
-(defgeneric geev (a &optional job)
- (: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.
-"))
-
-
-(defmethod geev :before ((a standard-matrix) &optional (job :NN))
- (if (not (square-matrix-p a))
- (error "argument A given to GEEV is not a square matrix")
- (if (not (member job '(:nn nn :vn vn :nv nv :vv vv t)))
- (error "argument JOB given to GEEV is not recognized"))))
-
-(defun geev-fix-up-eigvec (n real-eig-p eigval eigvec)
- (if real-eig-p
- (make-instance 'real-matrix :nrows n :ncols n :store eigvec)
- ;; We have to carefully handle complex-valued eigenvectors and eigenvalues
- (let ((evec (make-complex-matrix n n)))
- (do ((col 0 (incf col))
- (posn 0))
- ((>= col n) evec)
- (cond ((zerop (aref eigval col))
- (dotimes (row n)
- (setf (matrix-ref evec row col) (aref eigvec posn))
- (incf posn)))
- (t
- (dotimes (row n)
- (let ((next-posn (+ posn n)))
- (setf (matrix-ref evec row col)
- (complex (aref eigvec posn) (aref eigvec next-posn)))
- (setf (matrix-ref evec row (1+ col))
- (complex (aref eigvec posn) (- (aref eigvec next-posn))))
- (incf posn)))
- ;; Skip over the next column, which we've already used
- (incf col)
- (incf posn n)))))))
-
-(defun geev-fix-up-eigen (n wr wi vr vl left-eig right-eig)
- (let ((res nil)
- ;; Eigenvalues are real unless the max of wi is not zero.
- (real-eig (zerop (aref wi (1- (blas::idamax n wi 1))))))
-
- (when right-eig
- (push (geev-fix-up-eigvec n real-eig wi vr) res))
-
- (if real-eig
- (if (or right-eig left-eig)
- (let ((eigval (make-real-matrix n n)))
- (dotimes (k n)
- (setf (matrix-ref eigval k k) (aref wr k)))
- (push eigval res))
- (let ((eigval (make-real-matrix n 1)))
- (dotimes (k n)
- (setf (matrix-ref eigval k) (aref wr k)))
- (push eigval res)))
- (if (or right-eig left-eig)
- (let ((eigval (make-complex-matrix n n)))
- (dotimes (k n)
- (setf (matrix-ref eigval k k) (complex (aref wr k) (aref wi k))))
- (push eigval res))
- (let ((eigval (make-complex-matrix n 1)))
- (dotimes (k n)
- (setf (matrix-ref eigval k) (complex (aref wr k) (aref wi k))))
- (push eigval res))))
-
- (when left-eig
- (push (geev-fix-up-eigvec n real-eig wi vl) res))
-
- (push t res)
- (values-list (nreverse res))))
-
-
-(let ((work (allocate-real-store 1)))
- (defun dgeev-workspace-inquiry (n job)
- ;; Ask geev how much space it wants for the work array
- (multiple-value-bind (jobvl jobvr)
- (case job
- (:nn (values "N" "N"))
- ((:vn t) (values "N" "V"))
- (:nv (values "V" "N"))
- (: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-wr store-wi store-vl store-vr
- work info)
- (dgeev jobvl
- jobvr
- n ; N
- work ; A
- n ; LDA
- work ; WR
- work ; WI
- work ; VL
- ldvl ; LDVL
- work ; VR
- ldvr ; LDVR
- work ; WORK
- -1 ; LWORK
- 0 ) ; INFO
- (declare (ignore store-a store-wr store-wi store-vl store-vr))
- (assert (zerop info))
- (ceiling (realpart (aref work 0))))))))
-
-
-(defmethod geev ((a real-matrix) &optional (job :NN))
- (let* ((n (nrows a))
- (a (copy a))
- (xxx (allocate-real-store 1))
- (wr (allocate-real-store n))
- (wi (allocate-real-store n))
- (lwork (dgeev-workspace-inquiry n job))
- (work (allocate-real-store lwork)))
-
- (declare (type fixnum n)
- (type (simple-array real-matrix-element-type (*)) xxx wr wi))
-
- (case job
- (:nn
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "N" ; JOBVL
- "N" ; JOBVR
- n ; N
- (store a) ; A
- n ; LDA
- wr ; WR
- wi ; WI
- xxx ; VL
- 1 ; LDVL
- xxx ; VR
- 1 ; LDVR
- work ; WORK
- lwork ; LWORK
- 0 ) ; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl nil nil)
- (values nil nil))))
-
- ((:vn t)
- (let* ((vr (allocate-real-store (* n n))))
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "N" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- wr ;; WR
- wi ;; WI
- xxx ;; VL
- 1 ;; LDVL
- vr ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- 0 ) ;; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl nil t)
- (values nil nil)))))
-
- (:nv
- (let* ((vl (allocate-real-store (* n n))))
-
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "V" ;; JOBVL
- "N" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- wr ;; WR
- wi ;; WI
- vl ;; VL
- n ;; LDVL
- xxx ;; VR
- 1 ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- 0 ) ;; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl t nil)
- (values nil nil)))))
-
- (:vv
- (let* ((vl (allocate-real-store (* n n)))
- (vr (allocate-real-store (* n n))))
-
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "V" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- wr ;; WR
- wi ;; WI
- vl ;; VL
- n ;; LDVL
- vr ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- 0 ) ;; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl t t)
- (values nil nil)))))
-
-
- )))
-
-
-(let ((work (allocate-complex-store 1)))
- (defun zgeev-workspace-inquiry (n job)
- ;; Ask geev how much space it wants for the work array
- (multiple-value-bind (jobvl jobvr)
- (case job
- (:nn (values "N" "N"))
- ((:vn t) (values "N" "V"))
- (:nv (values "V" "N"))
- (: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
- n ; N
- work ; A
- n ; LDA
- work ; W
- work ; VL
- 1 ; LDVL
- work ; VR
- ldvr ; LDVR
- work ; WORK
- -1 ; LWORK
- work ; RWORK
- 0 ) ; INFO
- (declare (ignore store-a store-w store-vl store-vr info))
- ;; The desired size in in work[0], which we convert to an
- ;; integer.
- (ceiling (aref work 0)))))))
-
-;; Hmm, should this really be 4 (5) different methods, one for each
-;; possible value of job?
-
-(defmethod geev ((a complex-matrix) &optional (job :NN))
- (let* ((n (nrows a))
- (a (copy a))
- (w (make-complex-matrix-dim n 1))
- (xxx (allocate-complex-store 1))
- (lwork (zgeev-workspace-inquiry n job))
- (work (allocate-complex-store lwork))
- (rwork (allocate-complex-store n)))
-
- (declare (type fixnum lwork n)
- (type (simple-array complex-matrix-element-type (*)) xxx work)
- (type (simple-array real-matrix-element-type (*)) rwork))
-
- (case job
- (:nn
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "N" ; JOBVL
- "N" ; JOBVR
- n ; N
- (store a) ; A
- n ; LDA
- (store w) ; W
- xxx ; VL
- 1 ; LDVL
- xxx ; VR
- 1 ; LDVR
- work ; WORK
- lwork ; LWORK
- rwork ; RWORK
- 0 ) ; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values w t)
- (values nil nil))))
-
- ((:vn t)
- (let* ((vr (make-complex-matrix-dim n n)))
-
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "N" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- (store w) ;; W
- xxx ;; VL
- 1 ;; LDVL
- (store vr) ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- rwork ;; RWORK
- 0 ) ;; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values vr (diag w) t)
- (values nil nil)))))
-
- (:nv
- (let* ((vl (make-complex-matrix-dim n n)))
-
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "V" ;; JOBVL
- "N" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- (store w) ;; W
- (store vl) ;; VL
- n ;; LDVL
- xxx ;; VR
- 1 ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- rwork ;; RWORK
- 0 ) ;; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values (diag w) vl t)
- (values nil nil)))))
-
- (:vv
- (let* ((vr (make-complex-matrix-dim n n))
- (vl (make-complex-matrix-dim n n)))
-
-
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "V" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- (store w) ;; W
- (store vl) ;; VL
- n ;; LDVL
- (store vr) ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- rwork ;; RWORK
- 0 ) ;; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values vr (diag w) vl t)
- (values nil nil)))))
-
-
- )))
-
-
-
-(defun eig (a &optional (job :nn))
- "
- Syntax
- ======
- (EIG a [job])
-
- Purpose
- =======
- Computes the eigenvalues and left/right eigenvector of A.
-
- EIG is an alias for GEEV, for more help see GEEV.
-"
- (geev a job))
-----------------------------------------------------------------------
Summary of changes:
lib-src/gnuplot/gnuplot.lisp | 22 ++-
packages.lisp | 5 +-
src/base/einstein.lisp | 26 +--
src/foreign-core/lapack.lisp | 4 +-
src/lapack/eig.lisp | 23 +--
src/lapack/geev.lisp | 524 ------------------------------------------
src/lapack/geqr.lisp | 154 +++++--------
src/lapack/lu.lisp | 3 +
src/level-1/copy.lisp | 30 +++
src/level-2/gemv.lisp | 2 +-
src/level-3/gemm.lisp | 10 +-
src/special/map.lisp | 4 +
src/utilities/functions.lisp | 24 ++
src/utilities/macros.lisp | 55 +++++-
tests/loopy-tests.lisp | 21 +-
15 files changed, 224 insertions(+), 683 deletions(-)
delete mode 100644 src/lapack/geev.lisp
hooks/post-receive
--
matlisp
|
