Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-289-g1f45e5c
|
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 standard-tensor) (size &optional initial-element)
(let ((size-sym (gensym))
diff --git a/src/base/tweakable.lisp b/src/base/tweakable.lisp
index 86dca90..cab998f 100644
--- a/src/base/tweakable.lisp
+++ b/src/base/tweakable.lisp
@@ -28,6 +28,11 @@
which case, may lead to memory error. Use at your own risk.
")
+(defparameter *default-rcond* 1d-15
+ "
+ The default value of condition number to be used for
+ determining the rank of a matrix (used in gelsy).
+")
;;Level 1--------------------------------------------------------;;
(defparameter *real-l1-fcall-lb* 5000
"If the size of the array is less than this parameter, the
diff --git a/src/classes/numeric.lisp b/src/classes/numeric.lisp
index f360bf9..632a85a 100644
--- a/src/classes/numeric.lisp
+++ b/src/classes/numeric.lisp
@@ -64,6 +64,9 @@
(deft/method t/compute-store-size (sym complex-numeric-tensor) (size)
`(* 2 ,size))
+ (deft/method t/store-size (sym complex-numeric-tensor) (vec)
+ `(/ (length ,vec) 2))
+
(deft/method t/store-ref (sym complex-numeric-tensor) (store idx)
(let ((store-s (gensym))
(idx-s (gensym))
diff --git a/src/foreign-core/lapack.lisp b/src/foreign-core/lapack.lisp
index 1c94e74..1dc05ed 100644
--- a/src/foreign-core/lapack.lisp
+++ b/src/foreign-core/lapack.lisp
@@ -1578,9 +1578,9 @@
(m :integer :input)
(n :integer :input)
(nrhs :integer :input)
- (a (* :double-float) :input-output)
+ (a (* :double-float :inc head-a) :input-output)
(lda :integer :input)
- (b (* :double-float) :input-output)
+ (b (* :double-float :inc head-b) :input-output)
(ldb :integer :input)
(jpvt (* :integer) :input-output)
(rcond :double-float :input)
@@ -1865,3 +1865,140 @@
(b (* :complex-double-float :inc head-b) :input-output)
(ldb :integer :input)
(info :integer :output))
+
+(def-fortran-routine zgelsy :void
+ "
+ SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
+ WORK, LWORK, RWORK, INFO )
+ Purpose
+ =======
+
+ ZGELSY computes the minimum-norm solution to a complex 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 callthey 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 unitary 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 permutation of matrix B (the right hand side) is faster and
+ more simple.
+ 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.
+
+ Arguments
+ =========
+
+ M (input) INTEGER
+ The number of rows of the matrix A. M >= 0.
+
+ N (input) INTEGER
+ The number of columns of the matrix A. N >= 0.
+
+ NRHS (input) INTEGER
+ The number of right hand sides, i.e., the number of
+ columns of matrices B and X. NRHS >= 0.
+
+ A (input/output) COMPLEX*16 array, dimension (LDA,N)
+ On entry, the M-by-N matrix A.
+ On exit, A has been overwritten by details of its
+ complete orthogonal factorization.
+
+ LDA (input) INTEGER
+ The leading dimension of the array A. LDA >= max(1,M).
+
+ B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
+ On entry, the M-by-NRHS right hand side matrix B.
+ On exit, the N-by-NRHS solution matrix X.
+
+ LDB (input) INTEGER
+ The leading dimension of the array B. LDB >= max(1,M,N).
+
+ JPVT (input/output) INTEGER array, dimension (N)
+ On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
+ to the front of AP, otherwise column i is a free column.
+ On exit, if JPVT(i) = k, then the i-th column of A*P
+ was the k-th column of A.
+
+ RCOND (input) DOUBLE PRECISION
+ RCOND is used to determine the effective rank of A, which
+ is defined as the order of the largest leading triangular
+ submatrix R11 in the QR factorization with pivoting of A,
+ whose estimated condition number < 1/RCOND.
+
+ RANK (output) INTEGER
+ The effective rank of A, i.e., the order of the submatrix
+ R11. This is the same as the order of the submatrix T11
+ in the complete orthogonal factorization of A.
+
+ WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
+ On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+
+ LWORK (input) INTEGER
+ The dimension of the array WORK.
+ The unblocked strategy requires that:
+ LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
+ where MN = min(M,N).
+ The block algorithm requires that:
+ LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
+ where NB is an upper bound on the blocksize returned
+ by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
+ and ZUNMRZ.
+
+ If LWORK = -1, then a workspace query is assumedthe routine
+ only calculates the optimal size of the WORK array, returns
+ this value as the first entry of the WORK array, and no error
+ message related to LWORK is issued by XERBLA.
+
+ RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
+
+ INFO (output) INTEGER
+ = 0: successful exit
+ < 0: if INFO = -i, the i-th argument had an illegal value
+
+ 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
+"
+ (m :integer :input)
+ (n :integer :input)
+ (nrhs :integer :input)
+ (a (* :complex-double-float :inc head-a) :input-output)
+ (lda :integer :input)
+ (b (* :complex-double-float :inc head-b) :input-output)
+ (ldb :integer :input)
+ (jpvt (* :integer) :input-output)
+ (rcond :complex-double-float :input)
+ (rank :integer :output)
+ (work (* :complex-double-float) :workspace-output)
+ (lwork :integer :input)
+ (rwork (* :double-float) :workspace-output)
+ (info :integer :output))
diff --git a/src/lapack/eig.lisp b/src/lapack/eig.lisp
index f0a92da..ed67faf 100644
--- a/src/lapack/eig.lisp
+++ b/src/lapack/eig.lisp
@@ -45,7 +45,7 @@
,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)
+ ,work (t/store-size ,sym ,work)
0
(the index-type (head ,A)) (if ,vl (the index-type (head ,vl)) 0) (if ,vr (the index-type (head ,vr)) 0)))))
;;
-----------------------------------------------------------------------
Summary of changes:
matlisp.asd | 9 ++-
packages.lisp | 2 +-
src/base/template.lisp | 4 +
src/base/tweakable.lisp | 5 ++
src/classes/numeric.lisp | 3 +
src/foreign-core/lapack.lisp | 141 +++++++++++++++++++++++++++++++++++++-
src/lapack/eig.lisp | 2 +-
src/lapack/gels.lisp | 145 ---------------------------------------
src/lapack/least-squares.lisp | 151 +++++++++++++++++++++++++++++++++++++++++
src/lapack/lu.lisp | 76 ++-------------------
src/level-1/axpy.lisp | 4 +
src/level-1/copy.lisp | 60 ++++++++++++++++
src/reader/loadsave.lisp | 85 +++++++++++++++++++++++
src/special/map.lisp | 14 +++-
14 files changed, 474 insertions(+), 227 deletions(-)
delete mode 100644 src/lapack/gels.lisp
create mode 100644 src/lapack/least-squares.lisp
create mode 100644 src/reader/loadsave.lisp
hooks/post-receive
--
matlisp
|
