Menu
▾
▴
matlisp-commit
|
From: Akshay S. <ak...@us...> - 2013-11-15 10:50:20
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, tensor has been updated
via 1f45e5ca07fb6ec6e83117fdb4a3ded5fa3e2b4f (commit)
via e357e8266b9d8c1590fff4177057582577277846 (commit)
via 98719db2d568c5b022b93d2520db95e89f210d77 (commit)
via 705275944ef0cbd6caad409f5a3b3148641fca32 (commit)
via 767a08754c5f93918b9bba5e3503f0286191f179 (commit)
via 2b4808f24ea2cb5413b5af069f044a2f1ac1eef2 (commit)
via 66de9b29137420f14bf37f794dbec1129664676f (commit)
from 5481334c9a288f9ced6967f1995d4d30d0e39e2f (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 1f45e5ca07fb6ec6e83117fdb4a3ded5fa3e2b4f
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Oct 27 20:46:26 2013 -0700
Made changes to loadsave.
diff --git a/src/reader/loadsave.lisp b/src/reader/loadsave.lisp
index 2437607..5d51b9c 100644
--- a/src/reader/loadsave.lisp
+++ b/src/reader/loadsave.lisp
@@ -46,18 +46,40 @@ returning two values: the string and the number of bytes read."
(multiple-value-bind (lns nrows) (split-seq #'(lambda (x) (member x newlines)) f-string)
(unless (null lns)
(let* ((ncols (second (multiple-value-list (split-seq #'(lambda (x) (member x delimiters)) (car lns)))))
- (ret (zeros (list nrows ncols) 'real-tensor)))
- (loop :for line :in lns
- :for i := 0 :then (1+ i)
- :do (loop :for num :in (split-seq #'(lambda (x) (member x delimiters)) line)
- :for j := 0 :then (1+ j)
- :do (setf (ref ret i j) (t/coerce (t/field-type real-tensor) (read-from-string num)))))
+ (ret (zeros (if (> ncols 1) (list nrows ncols) (list nrows)) 'real-tensor)))
+ (if (> ncols 1)
+ (loop :for line :in lns
+ :for i := 0 :then (1+ i)
+ :do (loop :for num :in (split-seq #'(lambda (x) (member x delimiters)) line)
+ :for j := 0 :then (1+ j)
+ :do (setf (ref ret i j) (t/coerce (t/field-type real-tensor) (read-from-string num)))))
+ (loop :for line :in lns
+ :for i := 0 :then (1+ i)
+ :do (setf (ref ret i) (t/coerce (t/field-type real-tensor) (read-from-string (car (split-seq #'(lambda (x) (member x delimiters)) line)))))))
ret)))))
(defun savetxt (fname mat &key (delimiter #\Tab) (newline #\Newline))
(with-open-file (out fname :direction :output :if-exists :overwrite :if-does-not-exist :create)
- (let ((ncols (ncols mat)))
- (loop :for i :from 0 :below (nrows mat)
- :do (loop :for j :from 0 :below ncols
- :do (format out "~a~a" (ref mat i j) (if (= j (1- ncols)) newline delimiter)))))
+ (cond
+ ((tensor-matrixp mat)
+ (let ((ncols (ncols mat)))
+ (loop :for i :from 0 :below (nrows mat)
+ :do (loop :for j :from 0 :below ncols
+ :do (format out "~a~a" (ref mat i j) (if (= j (1- ncols)) newline delimiter))))))
+ ((tensor-vectorp mat)
+ (loop :for i :from 0 :below (aref (dimensions mat) 0)
+ :do (format out "~a~a" (ref mat i) newline)))
+ (t
+ (let ((dims (dimensions mat))
+ (strd (strides mat)))
+ (format out ":head~a~a~a" delimiter (head mat) newline)
+ (format out ":dimensions~a" delimiter)
+ (loop :for i :from 0 :below (length dims)
+ :do (format out "~a~a" (aref dims i) (if (= i (1- (length dims))) newline delimiter)))
+ (format out ":strides~a" delimiter)
+ (loop :for i :from 0 :below (length dims)
+ :do (format out "~a~a" (aref strd i) (if (= i (1- (length dims))) newline delimiter)))
+ (let ((sto (store mat)))
+ (loop :for i :from 0 :below (length sto)
+ :do (format out "~a~a" (aref sto i) (if (= i (1- (length sto))) newline delimiter)))))))
nil))
diff --git a/src/special/map.lisp b/src/special/map.lisp
index 9dcb66b..1edc1d2 100644
--- a/src/special/map.lisp
+++ b/src/special/map.lisp
@@ -67,15 +67,18 @@
:collect (prog1 (funcall func (sub-tensor~ v-x nil))
(incf (slot-value v-x 'head) st-x))))))
-(defmacro tensor-foldl (type func ten init &optional (init-type (field-type type)))
+(defmacro tensor-foldl (type func ten init &key (init-type (field-type type)) (key nil))
(using-gensyms (decl (ten init))
- (with-gensyms (sto idx of funcsym)
+ (with-gensyms (sto idx of funcsym keysym)
`(let* (,@decl
,@(unless (symbolp func)
`((,funcsym ,func)))
+ ,@(unless (symbolp key)
+ `((,keysym ,key)))
(,sto (store ,ten)))
(declare (type ,type ,ten)
,@(unless (symbolp func) `((type function ,funcsym)))
+ ,@(unless (symbolp key) `((type function ,keysym)))
(type ,(store-type type) ,sto)
,@(when init-type
`((type ,init-type ,init))))
@@ -85,5 +88,10 @@
(,of (strides ,ten)))
:do (setf ,init (,@(if (symbolp func)
`(,func)
- `(funcall ,funcsym)) ,init (t/store-ref ,type ,sto ,of)))))
+ `(funcall ,funcsym)) ,init ,(recursive-append
+ (when key
+ (if (symbolp key)
+ `(,key)
+ `(funcall ,keysym)))
+ `(t/store-ref ,type ,sto ,of))))))
,init))))
commit e357e8266b9d8c1590fff4177057582577277846
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Oct 27 03:09:35 2013 -0700
Added reader/loadsave.lisp.
diff --git a/matlisp.asd b/matlisp.asd
index 0afc593..8b57116 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -181,10 +181,10 @@
(:file "mplusminus")
#+nil
(:file "mtimesdivide")))
- #+nil
(:module "matlisp-reader"
- :pathname "reader"
- :components ((:file "infix")))))
+ :pathname "reader"
+ :components (#+nil(:file "infix")
+ (:file "loadsave")))))
;; (defclass f2cl-cl-source-file (asdf:cl-source-file)
diff --git a/src/reader/loadsave.lisp b/src/reader/loadsave.lisp
new file mode 100644
index 0000000..2437607
--- /dev/null
+++ b/src/reader/loadsave.lisp
@@ -0,0 +1,63 @@
+(in-package #:matlisp)
+
+#+(not :sbcl)
+(defun file->string (path)
+ "Sucks up an entire file from PATH into a freshly-allocated string,
+returning two values: the string and the number of bytes read."
+ (declare (optimize (safety 0) (speed 3)))
+ (with-open-file (s path :external-format :iso8859-1)
+ (let* ((len (file-length s))
+ (data (make-array len :element-type 'standard-char)))
+ (values data (read-sequence data s)))))
+
+#+sbcl
+(defun file->string (path)
+"Sucks up an entire file from PATH into a freshly-allocated string,
+returning two values: the string and the number of bytes read."
+ (let* ((fsize (with-open-file (s path)
+ (file-length s)))
+ (data (make-array fsize :element-type 'standard-char))
+ (fd (sb-posix:open path 0)))
+ (unwind-protect (sb-posix:read fd (sb-sys:vector-sap data) fsize)
+ (sb-posix:close fd))
+ (values data fsize)))
+;;
+(definline split-seq (test seq &optional (filter-empty? t))
+ "Split a string, wherever the given character occurs."
+ (loop :for i :from (1- (length seq)) :downto -1
+ :with prev := (length seq)
+ :with split-list := nil
+ :with split-count := 0
+ :do (when (or (< i 0) (funcall test (aref seq i)))
+ (let ((str (subseq seq (1+ i) prev)))
+ (when (or (< (1+ i) prev) (not filter-empty?))
+ (incf split-count)
+ (push str split-list))
+ (setf prev i)))
+ :finally (return (values split-list split-count))))
+;;
+(defun splitlines (string)
+ "Split the given string wherever the Carriage-return occurs."
+ (split-seq #'(lambda (x) (or (char= x #\Newline) (char= x #\Return))) string))
+
+;;
+(defun loadtxt (fname &key (delimiters '(#\Space #\Tab #\,)) (newlines '(#\Newline #\;)))
+ (let* ((f-string (file->string fname)))
+ (multiple-value-bind (lns nrows) (split-seq #'(lambda (x) (member x newlines)) f-string)
+ (unless (null lns)
+ (let* ((ncols (second (multiple-value-list (split-seq #'(lambda (x) (member x delimiters)) (car lns)))))
+ (ret (zeros (list nrows ncols) 'real-tensor)))
+ (loop :for line :in lns
+ :for i := 0 :then (1+ i)
+ :do (loop :for num :in (split-seq #'(lambda (x) (member x delimiters)) line)
+ :for j := 0 :then (1+ j)
+ :do (setf (ref ret i j) (t/coerce (t/field-type real-tensor) (read-from-string num)))))
+ ret)))))
+
+(defun savetxt (fname mat &key (delimiter #\Tab) (newline #\Newline))
+ (with-open-file (out fname :direction :output :if-exists :overwrite :if-does-not-exist :create)
+ (let ((ncols (ncols mat)))
+ (loop :for i :from 0 :below (nrows mat)
+ :do (loop :for j :from 0 :below ncols
+ :do (format out "~a~a" (ref mat i j) (if (= j (1- ncols)) newline delimiter)))))
+ nil))
commit 98719db2d568c5b022b93d2520db95e89f210d77
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Oct 27 01:41:17 2013 -0700
Saving state.
diff --git a/src/lapack/gels.lisp b/src/lapack/gels.lisp
deleted file mode 100644
index 0d8008d..0000000
--- a/src/lapack/gels.lisp
+++ /dev/null
@@ -1,145 +0,0 @@
-;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
-
-(in-package "MATLISP")
-
-(defgeneric gelsy! (a b rcond)
- (:documentation "Destructive version of GELSY. See GELSY."))
-
-(defgeneric gelsy (a b rcond)
- (:documentation
- "
- Syntax
- =======
-
- (GELSY A B &key TOL)
-
- INPUT
- -----
- A A Matlisp matrix of size M x N
- B A Matlisp matrix of size M x P
- RCOND A condition number
-
- OUTPUT
- ------
- X A Matlisp matrix of size N x NRHS
- RANK An integer
-
- Purpose
- =======
-
- Compute the minimum-norm solution to a real linear least
- squares problem:
- minimize || A * X - B ||
- using a complete orthogonal factorization of A. A is an M-by-N
- matrix which may be rank-deficient.
-
- Several right hand side vectors b and solution vectors x can be
- handled in a single call; they are stored as the columns of the
- M-by-NRHS right hand side matrix B and the N-by-NRHS solution
- matrix X.
-
- The routine first computes a QR factorization with column pivoting:
- A * P = Q * [ R11 R12 ]
- [ 0 R22 ]
- with R11 defined as the largest leading submatrix whose estimated
- condition number is less than 1/RCOND. The order of R11, RANK,
- is the effective rank of A.
-
- Then, R22 is considered to be negligible, and R12 is annihilated
- by orthogonal transformations from the right, arriving at the
- complete orthogonal factorization:
- A * P = Q * [ T11 0 ] * Z
- [ 0 0 ]
- The minimum-norm solution is then
- X = P * Z' [ inv(T11)*Q1'*B ]
- [ 0 ]
- where Q1 consists of the first RANK columns of Q.
-
- This routine is basically identical to the original xGELSX except
- three differences:
- o The call to the subroutine xGEQPF has been substituted by the
- the call to the subroutine xGEQP3. This subroutine is a Blas-3
- version of the QR factorization with column pivoting.
- o Matrix B (the right hand side) is updated with Blas-3.
- o The permutation of matrix B (the right hand side) is faster and
- more simple.
-
- Further Details
- ===============
-
- Based on contributions by
- A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
- E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
- G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
-
- =====================================================================
-"))
-
-(defun check-info (info function-name)
- (unless (= info 0)
- (error "~a: error in argument ~d" function-name (- info))))
-
-(defun dgelsy-workspace-inquiry (m n nrhs a lda b ldb jpvt rcond rank)
- (let ((work (allocate-real-store 1)))
- (multiple-value-bind
- (store-a store-b store-jpvt rank store-work info)
- (lapack::dgelsy m n nrhs (store a) lda (store b) ldb jpvt rcond rank
- work -1 0)
- (declare (ignore store-a store-b store-jpvt rank store-work))
- (check-info info "dgelsy"))
- (values (ceiling (realpart (aref work 0))))))
-
-(defmethod gelsy! ((a real-matrix) (b real-matrix) rcond)
- (let* ((m (nrows a))
- (n (ncols a))
- (nrhs (ncols b))
- (jpvt (allocate-integer4-store n 0))
- (lda m)
- (ldb (max n m))
- (b-arg b))
- (when (and (< m n))
- ;; In this case we need to extend the matrix which stores B
- ;; since it will be used to store the computation result
- (setq b-arg (make-real-matrix n nrhs))
- (dotimes (i m)
- (dotimes (j nrhs)
- (setf (matrix-ref b-arg i j)
- (matrix-ref b i j)))))
- (let* ((lwork (dgelsy-workspace-inquiry m n nrhs a lda b-arg ldb jpvt
- rcond 0))
- (work (allocate-real-store lwork)))
- (assert (= m (nrows b)))
- (multiple-value-bind
- (store-a store-b store-jpvt rank store-work info)
- (lapack::dgelsy m n nrhs (store a) lda (store b-arg) ldb jpvt rcond 0
- work lwork 0)
- (declare (ignore store-a store-jpvt store-work))
- (check-info info "dgelsy")
- (let ((x (make-real-matrix n nrhs)))
- ;; extract the matrix X from B
- (dotimes (i n)
- (dotimes (j nrhs)
- (setf (matrix-ref x i j)
- (aref store-b (fortran-matrix-indexing i j n)))))
- (values x rank))))))
-
-(defmethod gelsy ((a real-matrix) (b real-matrix) rcond)
- (gelsy! (copy a) (copy b) rcond))
-
-;; Example
-
-#|
-
-(let* ((m 100)
- (n 100)
- (a (rand m n))
- (x (rand n 1))
- (b (m* a x))
- (eps (coerce (expt 2 -52) 'double-float))
- (rcond (* eps (max m n))))
- (multiple-value-bind (r1 rank)
- (matlisp::gelsy a b rcond)
- (list rank
- (norm (m- x r1)))))
-
-|#
diff --git a/src/lapack/lu.lisp b/src/lapack/lu.lisp
index fcd9d92..1c5e7f7 100644
--- a/src/lapack/lu.lisp
+++ b/src/lapack/lu.lisp
@@ -136,77 +136,11 @@
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)
- (error 'tensor-cannot-find-optimization :tensor-class tensor-class)))
- (matrix-class (getf opt :matrix)))
- `(defmethod lu ((A ,matrix-class) &optional (split-lu? t))
- (multiple-value-bind (lu ipiv info)
- (getrf! (with-order :col-major
- (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A)))))
- (declare (ignore info))
- (let* ((n (nrows a))
- (m (ncols a))
- (p (min n m)))
- (declare (type fixnum n m p))
- ;; Extract the lower triangular part, if requested
- (if split-lu?
- (if (= p m)
- (let*-typed ((umat (,(getf opt :zero-maker) (list p m)) :type ,matrix-class)
- ;;
- (u.rstd (row-stride umat) :type index-type)
- (u.cstd (col-stride umat) :type index-type)
- (u.of (head umat) :type index-type)
- (u.sto (store umat) :type ,(linear-array-type (getf opt :store-type)))
- ;;
- (lu.rstd (row-stride lu) :type index-type)
- (lu.cstd (col-stride lu) :type index-type)
- (lu.of (head lu) :type index-type)
- (lu.sto (store lu) :type ,(linear-array-type (getf opt :store-type))))
- (very-quickly
- (loop :for i :of-type index-type :from 0 :below p
- :do (let-typed ((lu.of-ii lu.of :type index-type))
- (loop :repeat (- m i)
- :do (progn
- (,(getf opt :reader-writer) lu.sto lu.of u.sto u.of)
- (,(getf opt :value-writer) (,(getf opt :fid+)) lu.sto lu.of)
- (incf lu.of lu.cstd)
- (incf u.of u.cstd)))
- (,(getf opt :value-writer) (,(getf opt :fid*)) lu.sto lu.of-ii)
- (incf lu.of (- lu.rstd (the index-type (* (- m i 1) lu.cstd))))
- (incf u.of (- u.rstd (the index-type (* (- m i 1) u.cstd)))))))
- (values lu umat ipiv))
- (let*-typed ((lmat (,(getf opt :zero-maker) (list n p)) :type ,matrix-class)
- ;;
- (l.rstd (row-stride lmat) :type index-type)
- (l.cstd (col-stride lmat) :type index-type)
- (l.of (head lmat) :type index-type)
- (l.sto (store lmat) :type ,(linear-array-type (getf opt :store-type)))
- ;;
- (lu.rstd (row-stride lu) :type index-type)
- (lu.cstd (col-stride lu) :type index-type)
- (lu.of (head lu) :type index-type)
- (lu.sto (store lu) :type ,(linear-array-type (getf opt :store-type))))
- (very-quickly
- (loop :for j :of-type index-type :from 0 :below p
- :do (progn
- (,(getf opt :value-writer) (,(getf opt :fid*)) l.sto l.of)
- (loop :repeat (- n j 1)
- :do (progn
- (incf lu.of lu.rstd)
- (incf l.of l.rstd)
- (,(getf opt :reader-writer) lu.sto lu.of l.sto l.of)
- (,(getf opt :value-writer) (,(getf opt :fid+)) lu.sto lu.of)))
- (incf lu.of (- lu.cstd (the index-type (* (- n j 2) lu.rstd))))
- (incf l.of (- l.cstd (the index-type (* (- n j 2) l.rstd)))))))
- (values lmat lu ipiv)))
- (values lu ipiv)))))))
-
-
+;; (defmethod lu ((a standard-tensor) &optional (split-lu? t))
+;; (multiple-value-bind (lu perm) (getrf! (copy a))
+;; (if (not split-lu?) (values lu perm)
+;; (let* ((min (lvec-min (dimensions lu)))
+;; (
;;
(deft/generic (t/lapack-getrs-func #'subtypep) sym ())
(deft/method t/lapack-getrs-func (sym real-tensor) ()
diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp
index d6dd42a..78cf827 100644
--- a/src/level-1/axpy.lisp
+++ b/src/level-1/axpy.lisp
@@ -165,3 +165,7 @@
")
(:method (alpha x (y standard-tensor))
(axpy! alpha x (copy y))))
+
+(defmethod axpy (alpha (x standard-tensor) (y (eql nil)))
+ (let ((tmp (zeros (dimensions x) (class-of x))))
+ (axpy! alpha x tmp)))
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index c06dcd2..c90bbff 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -120,6 +120,35 @@
;;
;;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?))
+(deft/method t/copy-triangle! (sym standard-tensor) (a b &optional (upper? t))
+ (using-gensyms (decl (diag 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))
+ (ecase ,diag
+ (t
+ (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)))
+ :do (t/store-set ,sym (t/store-ref ,sym ,sto-a ,of-a) ,sto-b ,of-b)))))))
+
+ ,b))))
+;;
(deft/generic (t/copy-triangle! #'subtypep) sym (a b &optional upper?))
(deft/method t/copy-triangle! (sym standard-tensor) (a b &optional (upper? t))
(using-gensyms (decl (a b))
@@ -162,6 +191,37 @@
:for ,of-b :of-type index-type := (head ,b) :then (+ ,of-b (:mark (lvec-foldr #'+ (strides ,b)) :type index-type))
:do (t/store-set ,sym ,@(if num? `(,a) `((t/store-ref ,sym ,sto-a ,of-a))) ,sto-b ,of-b)))))
,b))))
+
+;;
+(defgeneric copy-triangle! (x y &key upper? diag?)
+ (:method :before ((x standard-tensor) (y standard-tensor) &key upper? diag?)
+ (assert (and (tensor-matrixp x) (tensor-matrixp y)
+ (= (lvec-min (dimensions x)) (lvec-min (dimensions y))))
+ nil 'tensor-dimension-mismatch)))
+
+
+(defmethod copy-triangle! ((x standard-tensor) (y standard-tensor) &key (upper? t) (diag? t))
+ (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*)
+ (eql clx cly))
+ nil 'tensor-abstract-class :tensor-class (list clx cly))
+ (compile-and-eval
+ (let ((expr `())))
+ `(defmethod copy-triangle! ((x ,clx) (y ,cly) &key (upper? t) (diag? t))
+ (ecase diag?
+ (t ;;copy diagonal
+ (if upper? (t/copy-triangle! ,clx x y t) (t/copy-triangle! ,clx x y nil)))
+ (number
+ (let ((num (t/coerce ,(t/field-type clx) diag?)))
+ (if upper? (t/copy-triangle! ,clx x y t) (t/copy-triangle! ,clx x y nil))
+ (t/copy-diagonal! ,clx num y t)))
+ (nil
+ (let ((num
+
+|#
+
;;
(defmethod copy! :before ((x standard-tensor) (y standard-tensor))
(assert (very-quickly (lvec-eq (the index-store-vector (dimensions x)) (the index-store-vector (dimensions y)) #'=)) nil
commit 705275944ef0cbd6caad409f5a3b3148641fca32
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Oct 12 03:10:51 2013 -0700
Added least-squares to asd file.
diff --git a/matlisp.asd b/matlisp.asd
index b4b9ac0..0afc593 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -164,7 +164,8 @@
:depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
:components ((:file "lu")
(:file "chol")
- (:file "eig")))
+ (:file "eig")
+ (:file "least-squares")))
(:module "matlisp-special"
:pathname "special"
:depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
commit 767a08754c5f93918b9bba5e3503f0286191f179
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Oct 12 03:09:32 2013 -0700
Moved gels.lisp to least-squares.lisp, replaced its contents with new interface.
diff --git a/src/lapack/least-squares.lisp b/src/lapack/least-squares.lisp
new file mode 100644
index 0000000..7816998
--- /dev/null
+++ b/src/lapack/least-squares.lisp
@@ -0,0 +1,151 @@
+(in-package :matlisp)
+
+(deft/generic (t/lapack-gelsy-func #'subtypep) sym ())
+(deft/method t/lapack-gelsy-func (sym real-tensor) ()
+ 'dgelsy)
+
+(definline mzgelsy (m n nrhs a lda b ldb jpvt rcond rank work lwork info &optional (head-a 0) (head-b 0))
+ (zgelsy m n nrhs a lda b ldb jpvt rcond rank work lwork (t/store-allocator complex-tensor n) info head-a head-b))
+(deft/method t/lapack-gelsy-func (sym complex-tensor) ()
+ 'mzgelsy)
+;;
+(deft/generic (t/lapack-gelsy! #'subtypep) sym (A lda B ldb rcond work))
+(deft/method t/lapack-gelsy! (sym blas-numeric-tensor) (A lda B ldb rcond work)
+ (using-gensyms (decl (A lda B ldb rcond work))
+ (with-gensyms (jpvt)
+ `(let* (,@decl
+ (,jpvt (make-array (ncols ,A) :element-type '(unsigned-byte 32) :initial-element 0)))
+ (declare (type ,sym ,A ,B)
+ (type index-type ,lda ,ldb)
+ ;;BEWARE: This will throw an error, if you use (simple-array (complex double-float) (*)) for store.
+ (type ,(store-element-type sym) ,rcond)
+ (type ,(store-type sym) work)
+ (type (simple-array (unsigned-byte 32) (*)) ,jpvt))
+ (,(macroexpand-1 `(t/lapack-gelsy-func ,sym))
+ (nrows ,A) (ncols ,A) (ncols ,B)
+ (the ,(store-type sym) (store ,A)) ,lda
+ (the ,(store-type sym) (store ,B)) ,ldb
+ ,jpvt ,rcond 0
+ ,work (t/store-size ,sym ,work)
+ 0
+ (the index-type (head ,A)) (the index-type (head ,B)))))))
+
+(deft/generic (t/lapack-gelsy-workspace-inquiry #'subtypep) sym (m n nrhs))
+(deft/method t/lapack-gelsy-workspace-inquiry (sym blas-numeric-tensor) (m n nrhs)
+ (using-gensyms (decl (m n nrhs))
+ (with-gensyms (xxx)
+ `(let* (,@decl
+ (,xxx (t/store-allocator ,sym 1)))
+ (declare (type index-type ,m ,n ,nrhs)
+ (type ,(store-type sym) ,xxx))
+ (,(macroexpand-1 `(t/lapack-gelsy-func ,sym))
+ ,m ,n ,nrhs
+ ,xxx ,m
+ ,xxx ,m
+ ,xxx (t/coerce (t/store-element-type ,sym) 0) 0
+ ,xxx -1
+ 0)
+ (ceiling (t/frealpart ,(field-type sym) (t/store-ref ,sym ,xxx 0)))))))
+;;
+
+(defgeneric gelsy (A B &optional rcond)
+ (:documentation "
+ Syntax
+ =======
+
+ (GELSY A B &key TOL)
+
+ INPUT
+ -----
+ A A Matlisp matrix of size M x N
+ B A Matlisp matrix of size M x P
+ RCOND A condition number
+
+ OUTPUT
+ ------
+ X A Matlisp matrix of size N x NRHS
+ RANK An integer
+
+ Purpose
+ =======
+
+ Compute the minimum-norm solution to a real linear least
+ squares problem:
+ minimize || A * X - B ||
+ using a complete orthogonal factorization of A. A is an M-by-N
+ matrix which may be rank-deficient.
+
+ Several right hand side vectors b and solution vectors x can be
+ handled in a single call; they are stored as the columns of the
+ M-by-NRHS right hand side matrix B and the N-by-NRHS solution
+ matrix X.
+
+ The routine first computes a QR factorization with column pivoting:
+ A * P = Q * [ R11 R12 ]
+ [ 0 R22 ]
+ with R11 defined as the largest leading submatrix whose estimated
+ condition number is less than 1/RCOND. The order of R11, RANK,
+ is the effective rank of A.
+
+ Then, R22 is considered to be negligible, and R12 is annihilated
+ by orthogonal transformations from the right, arriving at the
+ complete orthogonal factorization:
+ A * P = Q * [ T11 0 ] * Z
+ [ 0 0 ]
+ The minimum-norm solution is then
+ X = P * Z' [ inv(T11)*Q1'*B ]
+ [ 0 ]
+ where Q1 consists of the first RANK columns of Q.
+
+ This routine is basically identical to the original xGELSX except
+ three differences:
+ o The call to the subroutine xGEQPF has been substituted by the
+ the call to the subroutine xGEQP3. This subroutine is a Blas-3
+ version of the QR factorization with column pivoting.
+ o Matrix B (the right hand side) is updated with Blas-3.
+ o The permutation of matrix B (the right hand side) is faster and
+ more simple.
+
+ Further Details
+ ===============
+
+ Based on contributions by
+ A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
+ E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
+ G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
+
+ =====================================================================
+")
+ (:method :before ((A standard-tensor) (B standard-tensor) &optional rcond)
+ (assert (and (tensor-matrixp A) (tensor-matrixp B) (= (nrows A) (nrows B))) nil 'tensor-dimension-mismatch)
+ (assert (or (null rcond) (> rcond 0)) nil 'invalid-value :expected '(> rcond 0) :given rcond :message "Invalid rcond.")))
+
+(defmethod gelsy ((A standard-tensor) (B standard-tensor) &optional (rcond *default-rcond*))
+ (let ((cla (class-name (class-of A)))
+ (clb (class-name (class-of B))))
+ (assert (and (member cla *tensor-type-leaves*) (member clb *tensor-type-leaves*))
+ nil 'tensor-abstract-class :tensor-class (list cla clb))
+ (cond
+ ((eql cla clb)
+ (compile-and-eval
+ `(defmethod gelsy ((oA ,cla) (B ,clb) &optional (rcond *default-rcond*))
+ (let* ((A (let ((*default-stride-ordering* :col-major)) (copy oA)))
+ (lwork (max (t/lapack-gelsy-workspace-inquiry ,cla (nrows A) (ncols A) (ncols B)) 1))
+ (work (t/store-allocator ,cla lwork)))
+ (declare (type index-type lwork)
+ (type ,(store-type cla) work)
+ (type ,cla A))
+ (let* ((rank-A 0)
+ (mn (max (nrows A) (ncols A)))
+ (X (let ((*default-stride-ordering* :col-major)) (zeros (list mn (ncols B)) ',cla))))
+ (copy! B (sub-tensor~ X `((0 * ,(nrows A)) (* * *)) t))
+ (multiple-value-bind (sto-a sto-b jpvt rank work-out info) (t/lapack-gelsy! ,cla A (or (blas-matrix-compatiblep A #\N) 0) X (or (blas-matrix-compatiblep X #\N) 0) rcond work)
+ ;;TODO: Implement inverse permutation-action, and return jpvt.
+ (declare (ignore sto-a sto-b work-out jpvt))
+ (setf rank-a rank)
+ (unless (= info 0)
+ (error "gelsy returned ~a." info)))
+ (values (copy (sub-tensor~ X `((0 * ,(ncols A)) (* * *)) t)) rank-a)))))
+ (gelsy A B rcond))
+ (t
+ (error "Don't know how to apply getrs! to classes ~a." (list cla clb))))))
commit 2b4808f24ea2cb5413b5af069f044a2f1ac1eef2
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Oct 12 03:01:41 2013 -0700
Fixed silly bug, related to output spacing.
diff --git a/src/foreign-core/lapack.lisp b/src/foreign-core/lapack.lisp
index 1dc05ed..5941cc4 100644
--- a/src/foreign-core/lapack.lisp
+++ b/src/foreign-core/lapack.lisp
@@ -2000,5 +2000,5 @@
(rank :integer :output)
(work (* :complex-double-float) :workspace-output)
(lwork :integer :input)
- (rwork (* :double-float) :workspace-output)
+ (rwork (* :double-float) :workspace) ;;Not workspace-output (!)
(info :integer :output))
commit 66de9b29137420f14bf37f794dbec1129664676f
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Oct 12 02:52:07 2013 -0700
o Added interface for gelsy.
diff --git a/packages.lisp b/packages.lisp
index b0aa7ca..6d8868c 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -154,7 +154,7 @@
#:dgeqrf #:zgeqrf #:dgeqp3 #:zgeqp3
#:dorgqr #:zungqr
#:dpotrs #:zpotrs #:dpotrf #:zpotrf
- #:dgelsy)
+ #:dgelsy #:zgelsy)
(:documentation "LAPACK routines"))
(defpackage "MATLISP-DFFTPACK"
diff --git a/src/base/template.lisp b/src/base/template.lisp
index cd7daf8..65a296c 100644
--- a/src/base/template.lisp
+++ b/src/base/template.lisp
@@ -130,6 +130,10 @@
(deft/method t/compute-store-size (sym standard-tensor) (size)
size)
+(deft/generic (t/store-size #'subtypep) sym (ele))
+(deft/method t/store-size (sym standard-tensor) (ele)
+ `(length ,ele))
+
(deft/generic (t/store-allocator #'subtypep) sym (size &optional initial-element))
(deft/method t/store-allocator (sym sta...
[truncated message content] |
Thanks for helping keep SourceForge clean.
X