Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-204-g468ad7e
|
From: Akshay S. <ak...@us...> - 2013-02-03 09:15:30
|
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 468ad7e97a5c82ae3ed1c39d47085811b719e8b0 (commit)
from d43798f53f1a103b043af6fa6742ea347c777a2f (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 468ad7e97a5c82ae3ed1c39d47085811b719e8b0
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Feb 3 01:10:31 2013 -0800
Fixed some bugs. Finally ported getrf(LU) of LAPACK.
diff --git a/matlisp.asd b/matlisp.asd
index 520d66c..e83b2ee 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -153,10 +153,10 @@
:pathname "level-3"
:depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1" "matlisp-level-2")
:components ((:file "gemm")))
- #+nil(:module "matlisp-lapack"
+ (:module "matlisp-lapack"
:pathname "lapack"
:depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
- :components ((:file "gesv")))
+ :components ((:file "getrf")))
#+nil
(:module "matlisp-sugar"
:pathname "sugar"
diff --git a/packages.lisp b/packages.lisp
index 5f02a07..64e09af 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -62,6 +62,7 @@
#:tensor-cannot-find-optimization #:tensor-class
#:tensor-dimension-mismatch
#:tensor-store-not-consecutive
+ #:tensor-method-does-not-exist
))
(defpackage "MATLISP-UTILITIES"
diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp
index a0d1eb0..778ceba 100644
--- a/src/base/permutation.lisp
+++ b/src/base/permutation.lisp
@@ -106,12 +106,10 @@
(declare (ignore arg))
(let ((len (length seq)))
(assert (>= len (permutation-size perm)) nil
- 'permutation-permute-error :seq-len len :per-size (permutation-size perm)))))
-
-;; (:method :before ((ten standard-tensor) (perm permutation) &optional (arg 0))
-;; (let ((len (aref (dimensions ten) arg)))
-;; (assert (>= len (permutation-size perm)) nil
-;; 'permutation-permute-error :seq-len len :permutation-size (permutation-size perm)))))
+ 'permutation-permute-error :seq-len len :per-size (permutation-size perm))))
+ (:method :before ((ten standard-tensor) (perm permutation) &optional (arg 0))
+ (assert (>= (aref (dimensions ten) arg) (permutation-size perm)) nil
+ 'permutation-permute-error :seq-len (aref (dimensions ten) arg) :permutation-size (permutation-size perm))))
(definline permute (thing perm &optional (arg 0))
(permute! (copy thing) perm arg))
@@ -136,9 +134,8 @@
:do (setf (aref seq i) (aref cseq (aref act i)))
:finally (return seq))))
-;; (defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0))
-;; (let ((cyc (action->cycle perm)))
-;; (permute! ten cyc arg)))
+(defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0))
+ (permute! ten (action->pivot-flip perm) arg))
;;Cycle
(definline apply-cycle! (seq pcyc)
@@ -168,29 +165,8 @@
:do (apply-cycle! seq cyc)))
seq)
-;; (defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0))
-;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:)))
-;; (rplaca (nthcdr arg slst) 0)
-;; (values (sub-tensor~ A slst) (sub-tensor~ A slst)))
-;; (let-typed ((cyclst (store perm) :type cons)
-;; (cp-ten (make-instance (class-of tone)
-;; :dimensions (copy-seq (dimensions tone))))
-;; (std-arg (aref (strides A) arg) :type index-type)
-;; (hd-sl (head ttwo) :type index-type))
-;; (dolist (cyc cyclst)
-;; (declare (type pindex-store-vector cyc))
-;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc))))))
-;; (copy! tone cp-ten)
-;; (loop :for i :of-type index-type :downfrom (1- (length cyc)) :to 0
-;; :do (progn
-;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc i))))
-;; (copy!
-;; (if (= i 0) cp-ten
-;; (progn
-;; (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i)))))
-;; ttwo))
-;; tone))))))
-;; A)
+(defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0))
+ (permute! A (action->pivot-flip (cycle->action perm)) arg))
;Pivot idx
(definline apply-flips! (seq pflip)
@@ -213,20 +189,21 @@
(copy-n cseq seq size))
seq)
-;; (defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0))
-;; (let ((idiv (store perm)))
-;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:)))
-;; (rplaca (nthcdr arg slst) 0)
-;; (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil)))
-;; (let ((argstd (aref (strides A) arg))
-;; (hd-sl (head ttwo)))
-;; (loop :for i :from 0 :below (length idiv)
-;; :do (progn
-;; (unless (= i (aref idiv i))
-;; (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd)))
-;; (swap! tone ttwo))
-;; (incf (head tone) argstd))))))
-;; A)
+(defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0))
+ (multiple-value-bind (t1 t2) (let ((slst (make-list (rank A) :initial-element '(* * *))))
+ (rplaca (nthcdr arg slst) (list 0 '* 1))
+ (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil)))
+ (let-typed ((argstd (aref (strides A) arg) :type index-type)
+ (hd-sl (head t2) :type index-type)
+ (idiv (store perm) :type pindex-store-vector))
+ (very-quickly
+ (loop :for i :from 0 :below (length idiv)
+ :do (progn
+ (unless (= i (aref idiv i))
+ (setf (head t2) (the index-type (+ hd-sl (the index-type (* (aref idiv i) argstd)))))
+ (swap! t1 t2))
+ (setf (head t1) (the index-type (+ argstd (the index-type (head t1))))))))))
+ A)
;;Conversions----------------------------------------------------;;
(defun action->cycle (act)
diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp
index 12eb8d2..4e8e80a 100644
--- a/src/base/standard-tensor.lisp
+++ b/src/base/standard-tensor.lisp
@@ -374,7 +374,7 @@
(declare (type index-type lidx))
(let-typed ((sto-x (store tensor) :type ,(linear-array-type store-element-type)))
(,reader sto-x lidx)))
- (defmethod (setf tensor-ref) (value (tensor ,tensor-class) lidx)
+ (defmethod (setf tensor-store-ref) (value (tensor ,tensor-class) lidx)
(declare (type index-type lidx))
(let-typed ((sto-x (store tensor) :type ,(linear-array-type store-element-type)))
(,value-writer (,coercer-unforgiving value) sto-x lidx)))
diff --git a/src/conditions.lisp b/src/conditions.lisp
index 862dab0..43bdacc 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -224,3 +224,11 @@
(:report (lambda (c stream)
(declare (ignore c))
(format stream "The strides of the store, of the given tensor are not conscutive."))))
+
+(define-condition tensor-method-does-not-exist (tensor-error)
+ ((method-name :reader method-name :initarg :method-name)
+ (tensor-class :reader tensor-class :initarg :tensor-class))
+ (:documentation "The method is not defined for the given tensor class.")
+ (:report (lambda (c stream)
+ (when (slots-boundp c 'method 'tensor-class)
+ (format stream "The tensor class: ~a, does not have the method: ~a defined as a specialization." (method-name c) (tensor-class c))))))
diff --git a/src/lapack/getrf.lisp b/src/lapack/getrf.lisp
index da4bea0..d5a17ae 100644
--- a/src/lapack/getrf.lisp
+++ b/src/lapack/getrf.lisp
@@ -3,14 +3,14 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Copyright (c) 2000 The Regents of the University of California.
-;;; All rights reserved.
-;;;
+;;; 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
@@ -44,11 +44,11 @@
(if (eq maj-A :col-major)
(multiple-value-bind (n-A n-ipiv info) (,lapack-func
(nrows A) (ncols A) (store A)
- ld-A (repr ipiv) 0)
+ ld-A (store ipiv) 0)
(declare (ignore n-A n-ipiv))
(assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info))
;;Convert back to 0-based indexing.
- (let-typed ((pidv (repr ipiv) :type perrepr-vector))
+ (let-typed ((pidv (store ipiv) :type pindex-store-vector))
(very-quickly
(loop :for i :of-type index-type :from 0 :below (length pidv)
:do (decf (aref pidv i)))))
@@ -74,54 +74,58 @@
Given an NxM matrix A, compute its LU factorization using
partial pivoting, row or column interchanges:
- A = P * L * U (if A is row-major ordered)
- A = L * U * P' (if A is col-major ordered)
+ A = P * L * U (if A is row-major ordered)
+ A = L * U * P' (if A is col-major ordered)
where:
- P: permutation matrix
- L: lower triangular with unit diagonal elements
- (lower trapezoidal when N>M)
- U: upper triangular
- (upper trapezoidal when N<M)
+ P: permutation matrix
+ L: lower triangular with unit diagonal elements
+ (lower trapezoidal when N>M)
+ U: upper triangular
+ (upper trapezoidal when N<M)
Return Values
=============
- [1] The factors L and U from the factorization A = P*L*U where the
+ [1] The factors L and U from the factorization A = P*L*U where the
unit diagonal elements of L are not stored. (overwriting A)
[2] IPIV
[3] INFO = T: successful
- i: U(i,i) is exactly zero.
+ i: U(i,i) is exactly zero.
"))
(defmethod getrf! ((A real-matrix))
- (let ((ipiv (make-instance 'permutation-pivot-flip
- :repr (perrepr-id-action (lvec-min (dimensions A))))))
+ (let ((ipiv (let* ((*check-after-initializing?* nil)
+ (ret (make-instance 'permutation-pivot-flip :store (pindex-id (lvec-min (dimensions A))))))
+ (setf (permutation-size ret) (length (store ret)))
+ ret)))
(real-typed-getrf! A ipiv)))
(defmethod getrf! ((A complex-matrix))
- (let ((ipiv (make-instance 'permutation-pivot-flip
- :repr (perrepr-id-action (lvec-min (dimensions A))))))
+ (let ((ipiv (let* ((*check-after-initializing?* nil)
+ (ret (make-instance 'permutation-pivot-flip :store (pindex-id (lvec-min (dimensions A))))))
+ (setf (permutation-size ret) (length (store ret)))
+ ret)))
(complex-typed-getrf! A ipiv)))
;;
-(defgeneric lu (a &key with-l with-u with-p)
+(defgeneric lu (a &optional split-lu?)
(:documentation
"
Syntax
======
(LU a [:WITH-P with-p] [:WITH-L with-l] [:WITH-U with-u])
-
+
Purpose
=======
- Computes the LU decomposition of A.
+ Computes the LU decomposition of A.
This functions is an interface to GETRF!
Return Values
=============
[1] the factors L,U from the factorization in a single matrix,
- where the unit diagonal elements of L are not stored
+ where the unit diagonal elements of L are not stored
[2]-[4] If WITH-X then X, in the order L,U,P
By default WITH-L,WITH-U,WITH-P.
@@ -130,81 +134,70 @@
;;Sure I can do this with a defmethod (slowly), but where's the fun in that ? :)
(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)))
+ (error 'tensor-cannot-find-optimization :tensor-class tensor-class)))
+ (matrix-class (getf opt :matrix)))
`(eval-when (:compile-toplevel :load-toplevel :execute)
- (defmethod lu ((A ,matrix-class) &key with-l with-u with-p)
+ (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* ((ret (list lu))
- (n (nrows a))
+ (let* ((n (nrows a))
(m (ncols a))
(p (min n m)))
(declare (type fixnum n m p))
;; Extract the lower triangular part, if requested
- (when with-l
- (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 :element-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 :element-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)))
- (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))))))
- (push lmat ret))))
- ;; Extract the upper triangular part, if requested
- (when with-u
- (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 :element-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 :element-type))))
- (progn
- (loop :for i :of-type index-type :from 0 :below p
- :do (progn
- (loop :repeat (- m i)
- :do (progn
- (,(getf opt :reader-writer) lu.sto lu.of u.sto u.of)
- (incf lu.of lu.cstd)
- (incf u.of u.cstd)))
- (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))))))
- (push umat ret))))
- ;; Extract the permutation matrix, if requested
- (if with-p
- (let* ((npiv (length ipiv))
- (pmat (make-real-matrix-dim n n))
- (pidx (make-array n :element-type '(unsigned-byte 32))))
- ;; Compute the P matrix from the pivot vector
- (dotimes (k n)
- (setf (aref pidx k) k))
- (dotimes (k npiv)
- (rotatef (aref pidx k) (aref pidx (1- (aref ipiv k)))))
- (dotimes (k n)
- (setf (matrix-ref pmat (aref pidx k) k) 1))
- (push pmat result)))
- ;; Return the final result
- (values-list (nreverse ret))))))))
+ (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))))))))
(make-lu real-tensor)
+(make-lu complex-tensor)
-----------------------------------------------------------------------
Summary of changes:
matlisp.asd | 4 +-
packages.lisp | 1 +
src/base/permutation.lisp | 69 ++++++-----------
src/base/standard-tensor.lisp | 2 +-
src/conditions.lisp | 8 ++
src/lapack/getrf.lisp | 171 ++++++++++++++++++++---------------------
6 files changed, 117 insertions(+), 138 deletions(-)
hooks/post-receive
--
matlisp
|
