Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-147-g8ccded8
|
From: Akshay S. <ak...@us...> - 2012-07-26 13:31:37
|
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 8ccded8d5db3d1918b7af875f4dbddd16dc75f28 (commit)
via 284a1e7bcb18ff7bc25e53d2b636d4fe5c963205 (commit)
via 855c687f17ce0468bd189e8c6f9942ad5cec2999 (commit)
via 9980ae3686cf6361c2e8d8dec95d85f355b3a5d8 (commit)
via 0fc0b662754ddb98367d6add3aeb42f71a9301aa (commit)
via 83ad581a242f7fd2c6416dc115192692a7447c35 (commit)
from 9a94775cd4eb5593fea88f5cf665bcadc198fb6f (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 8ccded8d5db3d1918b7af875f4dbddd16dc75f28
Author: Akshay Srinivasan <aks...@gm...>
Date: Thu Jul 26 18:56:33 2012 +0530
o Made copy and swap a generic method. Added support for permute! to handle tensors.
diff --git a/matlisp.asd b/matlisp.asd
index 3bbf437..d0e7c0f 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -40,10 +40,10 @@
((:file "packages")))
(asdf:defsystem matlisp-config
- :pathname #.(translate-logical-pathname "matlisp:builddir;")
- :depends-on ("matlisp-packages")
- :components
- ((:file "config")))
+ :pathname #.(translate-logical-pathname "matlisp:builddir;")
+ :depends-on ("matlisp-packages")
+ :components
+ ((:file "config")))
(asdf:defsystem matlisp-conditions
:depends-on ("matlisp-packages" "matlisp-config")
@@ -101,8 +101,12 @@
;;
(:file "loopy"
:depends-on ("standard-tensor"))
+ (:file "generic-copy"
+ :depends-on ("standard-tensor" "loopy"))
+ (:file "generic-swap"
+ :depends-on ("standard-tensor" "loopy"))
(:file "permutation"
- :depends-on ("standard-tensor"))
+ :depends-on ("standard-tensor" "generic-copy" "generic-swap"))
(:file "blas-helpers"
:depends-on ("standard-tensor" "permutation"))
(:file "print"
diff --git a/packages.lisp b/packages.lisp
index ca65b2c..e71bae6 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -33,9 +33,10 @@
;;<conditon {accessors*}>
;;Generic errors
#:generic-error #:message
+ #:dimension-mismatch
#:assumption-violated
#:invalid-type #:given #:expected
- #:invalid-arguments
+ #:invalid-arguments #:argnum
#:invalid-value #:given #:expected
#:unknown-token #:token
#:parser-error
diff --git a/src/base/blas-helpers.lisp b/src/base/blas-helpers.lisp
index 1e3fcb6..969c3ea 100644
--- a/src/base/blas-helpers.lisp
+++ b/src/base/blas-helpers.lisp
@@ -38,8 +38,8 @@
(defun blas-matrix-compatible-p (matrix op)
(declare (type standard-matrix matrix))
- (let ((rs (aref (strides matrix) 0))
- (cs (aref (strides matrix) 1)))
+ (let ((rs (row-stride matrix))
+ (cs (col-stride matrix)))
(declare (type index-type rs cs))
(cond
((= cs 1) (values :row-major rs (fortran-nop op)))
diff --git a/src/base/generic-copy.lisp b/src/base/generic-copy.lisp
new file mode 100644
index 0000000..4445153
--- /dev/null
+++ b/src/base/generic-copy.lisp
@@ -0,0 +1,72 @@
+(in-package #:matlisp)
+
+(defgeneric copy! (from-tensor to-tensor)
+ (:documentation
+ "
+ Syntax
+ ======
+ (COPY! x y)
+
+ Purpose
+ =======
+ Copies the contents of X into Y. Returns Y.
+
+ X,Y must have the same dimensions, and
+ ergo the same number of elements.
+
+ Furthermore, X may be a scalar, in which
+ case Y is filled with X.
+")
+ (:method :before ((x cons) (y cons))
+ (assert (= (length x) (length y))))
+ (:method :before ((x array) (y array))
+ (assert (subtypep (array-element-type x) (array-element-type y))
+ nil 'invalid-type
+ :given (array-element-type y) :expected (array-element-type x))
+ (assert (and
+ (= (array-rank x) (array-rank y))
+ (reduce #'(lambda (x y) (and x y))
+ (mapcar #'= (array-dimensions x) (array-dimensions y))))
+ nil 'dimension-mismatch))
+ (:method :before (x (y array))
+ (assert (subtypep (type-of x) (array-element-type y))
+ nil 'invalid-type
+ :given (type-of x) :expected (array-element-type x))))
+
+(defmethod copy! ((from cons) (to cons))
+ (let-rec cdr-writer ((flst from) (tlst to))
+ (if (null flst) to
+ (progn
+ (rplaca tlst (car flst))
+ (cdr-writer (cdr flst) (cdr tlst))))))
+
+(defmethod copy! ((from t) (to cons))
+ (mapl #'(lambda (lst) (rplaca lst from)) to)
+ to)
+
+(defmethod copy! ((from array) (to array))
+ (let ((lst (make-list (array-rank to))))
+ (mod-dotimes (idx (make-index-store (array-dimensions to)))
+ do (progn
+ (idx->list! idx lst)
+ (setf (apply #'aref to lst) (apply #'aref from lst)))))
+ to)
+
+;;
+(defgeneric copy (object)
+ (:documentation
+ "
+ Syntax
+ ======
+ (COPY x)
+
+ Purpose
+ =======
+ Return a copy of X"))
+
+(defmethod copy ((lst cons))
+ (copy-list lst))
+
+(defmethod copy ((arr array))
+ (let ((ret (make-array (array-dimensions arr) :element-type (array-element-type arr))))
+ (copy! arr ret)))
diff --git a/src/base/generic-swap.lisp b/src/base/generic-swap.lisp
new file mode 100644
index 0000000..058f3c1
--- /dev/null
+++ b/src/base/generic-swap.lisp
@@ -0,0 +1,19 @@
+(in-package #:matlisp)
+
+(defgeneric swap! (x y)
+ (:documentation
+"
+ Sytnax
+ ======
+ (SWAP! x y)
+
+ Purpose
+ =======
+ Given tensors X,Y, performs:
+
+ X <-> Y
+
+ and returns Y.
+
+ X, Y must have the same dimensions.
+"))
diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp
index 2e866ec..51a405f 100644
--- a/src/base/permutation.lisp
+++ b/src/base/permutation.lisp
@@ -158,22 +158,28 @@
(make-instance 'permutation-pivot-flip :repr pact))
;;Generic permute! method.
-(defgeneric permute! (seq perm)
+(defgeneric permute! (thing permutation &optional argument)
(:documentation "
- (permute! seq perm)
+ (permute! thing permutation [argument 0])
- Applies the permutation on the sequence.
-")
- (:method :before ((seq sequence) (perm permutation))
+ Permutes the ARGUMENT index of the the array-like object THING, by
+ applying PERMUTATION on it.")
+ (:method :before ((seq sequence) (perm permutation) &optional (arg 0))
+ (declare (ignore arg))
(let ((len (length seq)))
(assert (>= len (group-rank perm)) nil
+ 'permutation-permute-error :seq-len len :group-rank (group-rank perm))))
+ (:method :before ((ten standard-tensor) (perm permutation) &optional (arg 0))
+ (let ((len (aref (dimensions ten) arg)))
+ (assert (>= len (group-rank perm)) nil
'permutation-permute-error :seq-len len :group-rank (group-rank perm)))))
-
-(definline permute (seq perm)
- (permute! (copy-seq seq) perm))
+
+(definline permute (thing perm &optional (arg 0))
+ (permute! (copy thing) perm arg))
;;Action
-(defmethod permute! ((seq cons) (perm permutation-action))
+(defmethod permute! ((seq cons) (perm permutation-action) &optional arg)
+ (declare (ignore arg))
(let ((cseq (make-array (length seq) :initial-contents seq))
(act (repr perm))
(glen (group-rank perm)))
@@ -185,7 +191,8 @@
(rplaca x (aref cseq (aref act i)))
(incf i)))) seq)))
-(defmethod permute! ((seq vector) (perm permutation-action))
+(defmethod permute! ((seq vector) (perm permutation-action) &optional arg)
+ (declare (ignore arg))
(let ((cseq (make-array (length seq) :initial-contents seq))
(act (repr perm)))
(loop
@@ -205,7 +212,8 @@
(if (= i 0) xl
(aref seq (aref pcyc (1- i))))))))
-(defmethod permute! ((seq cons) (perm permutation-cycle))
+(defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg)
+ (declare (ignore arg))
(let ((cseq (make-array (length seq) :initial-contents seq))
(glen (group-rank perm)))
(dolist (cyc (repr perm))
@@ -218,20 +226,23 @@
(rplaca x (aref cseq i))
(incf i)))) seq)))
-(defmethod permute! ((seq vector) (perm permutation-cycle))
+(defmethod permute! ((seq vector) (perm permutation-cycle) &optional arg)
+ (declare (ignore arg))
(dolist (cyc (repr perm) seq)
(declare (type perrepr-vector cyc))
(apply-cycle! seq cyc)))
;;Pivot idx
-(defmethod permute! ((seq vector) (perm permutation-pivot-flip))
+(defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg)
+ (declare (ignore arg))
(let-typed ((pidx (repr perm) :type perrepr-vector))
(loop for i of-type index-type from 0 below (group-rank perm)
unless (= i (aref pidx i))
do (rotatef (aref seq i) (aref seq (aref pidx i)))
finally (return seq))))
-(defmethod permute! ((seq cons) (perm permutation-pivot-flip))
+(defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg)
+ (declare (ignore arg))
(let ((cseq (make-array (length seq) :initial-contents seq))
(glen (group-rank perm)))
(permute! cseq perm)
@@ -242,6 +253,20 @@
(rplaca x (aref cseq i))
(incf i)))) seq)))
+(defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0))
+ (let ((idiv (repr 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)))
+ (loop for i from 0 below (length idiv)
+ do (progn
+ (unless (= i (aref idiv i))
+ (setf (head ttwo) (* (aref idiv i) argstd))
+ (swap! tone ttwo))
+ (incf (head tone) argstd))))))
+ A)
+
;;Conversions----------------------------------------------------;;
(defun action->cycle (act)
"
@@ -412,6 +437,5 @@
(qsort-bounds todo)))))))
(qsort-bounds `((0 ,len)))
(values seq (action->cycle (make-paction perm))))))
-
;;Add a general sorter, this is a very useful thing to have.
;;Add a function to apply permutations to a matrices, tensors.
diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp
index c27af81..3ff61e3 100644
--- a/src/base/standard-tensor.lisp
+++ b/src/base/standard-tensor.lisp
@@ -57,6 +57,13 @@
(loop for ele across a
collect ele))
+(definline idx->list! (a lst)
+ ;;No error checking!
+ (mapl (let ((i 0))
+ #'(lambda (lst)
+ (rplaca lst (aref a i))
+ (incf i)))
+ lst))
;;
(defclass standard-tensor ()
((rank
@@ -72,6 +79,10 @@
:accessor number-of-elements
:type index-type
:documentation "Total number of elements in the tensor.")
+ (element-type
+ :accessor element-type
+ :type symbol
+ :documentation "Element type of the tensor")
;;
(parent-tensor
:accessor parent-tensor
diff --git a/src/classes/complex-tensor.lisp b/src/classes/complex-tensor.lisp
index 9692b49..5e689f0 100644
--- a/src/classes/complex-tensor.lisp
+++ b/src/classes/complex-tensor.lisp
@@ -39,9 +39,8 @@
;;
(defclass complex-tensor (standard-tensor)
- ((store
- :initform nil
- :type complex-store-vector))
+ ((store :type complex-store-vector)
+ (element-type :initform 'complex-type))
(:documentation "Tensor class with complex elements."))
(defclass complex-matrix (standard-matrix complex-tensor)
diff --git a/src/classes/real-tensor.lisp b/src/classes/real-tensor.lisp
index 6ece4df..77c850b 100644
--- a/src/classes/real-tensor.lisp
+++ b/src/classes/real-tensor.lisp
@@ -22,9 +22,8 @@ Allocates real storage. Default initial-element = 0d0.")
;;
(defclass real-tensor (standard-tensor)
- ((store
- :initform nil
- :type real-store-vector))
+ ((store :type real-store-vector)
+ (element-type :initform 'real-type))
(:documentation "Tensor class with real elements."))
(defclass real-matrix (standard-matrix real-tensor)
diff --git a/src/conditions.lisp b/src/conditions.lisp
index 0cbcb67..8a05062 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -18,12 +18,23 @@
,@rest)
,@(loop for mth in methods
collect `(defmethod ,@(cdr mth)))))))
+
+(defun slots-boundp (obj &rest slots)
+ (dolist (slot slots t)
+ (unless (slot-boundp obj slot)
+ (return nil))))
;;Generic conditions---------------------------------------------;;
(defcondition generic-error (error)
((message :reader message :initarg :message :initform ""))
(:method print-object ((c generic-error) stream)
(format stream (message c))))
+(defcondition dimension-mismatch (generic-error)
+ ()
+ (:method print-object ((c generic-error) stream)
+ (format stream "Dimension mismatch.")
+ (call-next-method)))
+
(defcondition assumption-violated (generic-error)
()
(:method print-object ((c assumption-violated) stream)
@@ -35,26 +46,33 @@
(expected-type :reader expected :initarg :expected))
(:documentation "Given an unexpected type.")
(:method print-object ((c invalid-type) stream)
- (format stream "Given object of type ~A, expected ~A.~%" (given c) (expected c))
+ (when (slots-boundp c 'given-type 'expected-type)
+ (format stream "Given object of type ~A, expected ~A.~%" (given c) (expected c)))
(call-next-method)))
(defcondition invalid-arguments (generic-error)
- ()
- (:documentation "Given invalid arguments to the function."))
+ ((argument-number :reader argnum :initarg :argnum))
+ (:documentation "Given invalid arguments to the function.")
+ (:method print-object ((c invalid-arguments) stream)
+ (when (slot-boundp c 'argument-number)
+ (format stream "The ~a'th argument given to the function is invalid." (argnum c)))
+ (call-next-method)))
(defcondition invalid-value (generic-error)
((given-value :reader given :initarg :given)
(expected-value :reader expected :initarg :expected))
(:documentation "Given an unexpected value.")
(:method print-object ((c invalid-value) stream)
- (format stream "Given object ~A, expected ~A.~%" (given c) (expected c))
+ (when (slots-boundp c 'given-value 'expected-value)
+ (format stream "Given object ~A, expected ~A.~%" (given c) (expected c)))
(call-next-method)))
(defcondition unknown-token (generic-error)
((token :reader token :initarg :token))
(:documentation "Given an unknown token.")
(:method print-object ((c unknown-token) stream)
- (format stream "Given unknown token: ~A.~%" (token c))
+ (when (slot-boundp c 'token)
+ (format stream "Given unknown token: ~A.~%" (token c)))
(call-next-method)))
(defcondition parser-error (generic-error)
@@ -66,7 +84,8 @@
(to :reader to :initarg :to))
(:documentation "Cannot coerce one type into another.")
(:method print-object ((c coercion-error) stream)
- (format stream "Cannot coerce ~a into ~a.~%" (from c) (to c))
+ (when (slots-boundp c 'from 'to)
+ (format stream "Cannot coerce ~a into ~a.~%" (from c) (to c)))
(call-next-method)))
(defcondition out-of-bounds-error (generic-error)
@@ -74,7 +93,8 @@
(bound :reader bound :initarg :bound))
(:documentation "General out-of-bounds error")
(:method print-object ((c out-of-bounds-error) stream)
- (format stream "Out-of-bounds error, requested index : ~a, bound : ~a.~%" (requested c) (bound c))
+ (when (slots-boundp c 'requested 'bound)
+ (format stream "Out-of-bounds error, requested index : ~a, bound : ~a.~%" (requested c) (bound c)))
(call-next-method)))
(defcondition non-uniform-bounds-error (generic-error)
@@ -82,7 +102,8 @@
(found :reader found :initarg :found))
(:documentation "Bounds are not uniform")
(:method print-object ((c non-uniform-bounds-error) stream)
- (format stream "The bounds are not uniform, assumed bound : ~a, now found to be : ~a.~%" (assumed c) (found c))
+ (when (slots-boundp c 'assumed 'found)
+ (format stream "The bounds are not uniform, assumed bound : ~a, now found to be : ~a.~%" (assumed c) (found c)))
(call-next-method)))
;;Permutation conditions-----------------------------------------;;
@@ -100,10 +121,10 @@
((sequence-length :reader seq-len :initarg :seq-len)
(group-rank :reader group-rank :initarg :group-rank))
(:documentation "Cannot permute sequence.")
- (:report (lambda (c stream)
- (format stream "Cannot permute sequence.
-sequence-length : ~a
-group-rank: ~a" (seq-len c) (group-rank c)))))
+ (:report (lambda (c stream)
+ (format stream "Cannot permute sequence.")
+ (when (slots-boundp c 'sequence-length 'group-rank)
+ (format stream "~%sequence-length : ~a group-rank: ~a" (seq-len c) (group-rank c))))))
;;Tensor conditions----------------------------------------------;;
(define-condition tensor-error (error)
@@ -115,26 +136,30 @@ group-rank: ~a" (seq-len c) (group-rank c)))))
(store-size :reader store-size :initarg :store-size))
(:documentation "An out of bounds index error for the one-dimensional store.")
(:report (lambda (c stream)
- (format stream "Requested index ~A, but store is only of size ~A." (index c) (store-size c)))))
+ (when (slots-boundp c 'index 'store-size)
+ (format stream "Requested index ~A, but store is only of size ~A." (index c) (store-size c))))))
(define-condition tensor-insufficient-store (tensor-error)
((store-size :reader store-size :initarg :store-size)
(max-idx :reader max-idx :initarg :max-idx))
(:documentation "Store is too small for the tensor with given dimensions.")
(:report (lambda (c stream)
- (format stream "Store size is ~A, but maximum possible index is ~A." (store-size c) (max-idx c)))))
+ (when (slots-boundp c 'max-idx 'store-size)
+ (format stream "Store size is ~A, but maximum possible index is ~A." (store-size c) (max-idx c))))))
(define-condition tensor-not-matrix (tensor-error)
((tensor-rank :reader rank :initarg :rank))
(:documentation "Given tensor is not a matrix.")
(:report (lambda (c stream)
- (format stream "Given tensor with rank ~A, is not a matrix." (rank c)))))
+ (when (slots-boundp c 'tensor-rank)
+ (format stream "Given tensor with rank ~A, is not a matrix." (rank c))))))
(define-condition tensor-not-vector (tensor-error)
((tensor-rank :reader rank :initarg :rank))
(:documentation "Given tensor is not a vector.")
(:report (lambda (c stream)
- (format stream "Given tensor with rank ~A, is not a vector." (rank c)))))
+ (when (slots-boundp c 'tensor-rank)
+ (format stream "Given tensor with rank ~A, is not a vector." (rank c))))))
(define-condition tensor-index-out-of-bounds (tensor-error)
((argument :reader argument :initarg :argument)
@@ -142,46 +167,53 @@ group-rank: ~a" (seq-len c) (group-rank c)))))
(argument-space-dimension :reader dimension :initarg :dimension))
(:documentation "An out of bounds index error")
(:report (lambda (c stream)
- (format stream "~&Out of bounds for argument ~A: requested ~A, but dimension is only ~A." (argument c) (index c) (dimension c)))))
+ (when (slots-boundp c 'argument 'index 'argument-space-dimension)
+ (format stream "~&Out of bounds for argument ~A: requested ~A, but dimension is only ~A." (argument c) (index c) (dimension c))))))
(define-condition tensor-index-rank-mismatch (tensor-error)
((index-rank :reader index-rank :initarg :index-rank)
(rank :reader rank :initarg :rank))
(:documentation "Incorrect number of subscripts for the tensor.")
(:report (lambda (c stream)
- (format stream "Index is of size ~A, whereas the tensor is of rank ~A." (index-rank c) (rank c)))))
+ (when (slots-boundp c 'index-rank 'rank)
+ (format stream "Index is of size ~A, whereas the tensor is of rank ~A." (index-rank c) (rank c))))))
(define-condition tensor-invalid-head-value (tensor-error)
((head :reader head :initarg :head))
(:documentation "Incorrect value for the head of the tensor storage.")
(:report (lambda (c stream)
- (format stream "Head of the store must be >= 0, initialized with ~A." (head c)))))
+ (when (slots-boundp c 'head)
+ (format stream "Head of the store must be >= 0, initialized with ~A." (head c))))))
(define-condition tensor-invalid-dimension-value (tensor-error)
((argument :reader argument :initarg :argument)
(argument-dimension :reader dimension :initarg :dimension))
(:documentation "Incorrect value for one of the dimensions of the tensor.")
(:report (lambda (c stream)
- (format stream "Dimension of argument ~A must be > 0, initialized with ~A." (argument c) (dimension c)))))
+ (when (slots-boundp c 'argument 'argument-dimension)
+ (format stream "Dimension of argument ~A must be > 0, initialized with ~A." (argument c) (dimension c))))))
(define-condition tensor-invalid-stride-value (tensor-error)
((argument :reader argument :initarg :argument)
(argument-stride :reader stride :initarg :stride))
(:documentation "Incorrect value for one of the strides of the tensor storage.")
(:report (lambda (c stream)
- (format stream "Stride of argument ~A must be >= 0, initialized with ~A." (argument c) (stride c)))))
+ (when (slots-boundp c 'argument 'argument-stride)
+ (format stream "Stride of argument ~A must be >= 0, initialized with ~A." (argument c) (stride c))))))
(define-condition tensor-cannot-find-counter-class (tensor-error)
((tensor-class :reader tensor-class :initarg :tensor-class))
(:documentation "Cannot find the counter-class list of the given tensor class")
(:report (lambda (c stream)
- (format stream "Cannot find the counter-class list of the given tensor class: ~a." (tensor-class c)))))
+ (when (slots-boundp c 'tensor-class)
+ (format stream "Cannot find the counter-class list of the given tensor class: ~a." (tensor-class c))))))
(define-condition tensor-cannot-find-optimization (tensor-error)
((tensor-class :reader tensor-class :initarg :tensor-class))
(:documentation "Cannot find optimization information for the given tensor class")
(:report (lambda (c stream)
- (format stream "Cannot find optimization information for the given tensor class: ~a." (tensor-class c)))))
+ (when (slots-boundp c 'tensor-class)
+ (format stream "Cannot find optimization information for the given tensor class: ~a." (tensor-class c))))))
(define-condition tensor-dimension-mismatch (tensor-error)
()
diff --git a/src/ffi/c-ffi.lisp b/src/ffi/c-ffi.lisp
index 109a756..a7fd066 100644
--- a/src/ffi/c-ffi.lisp
+++ b/src/ffi/c-ffi.lisp
@@ -7,7 +7,7 @@
(real ,base-type)
(imag ,base-type)))
-(defccomplex %c.complex-double :double)R
+(defccomplex %c.complex-double :double)
(defccomplex %c.complex-float :float)
;; Get the equivalent CFFI type.
diff --git a/src/foreign-core/lapack.lisp b/src/foreign-core/lapack.lisp
index 9ce7a66..ad75f1b 100644
--- a/src/foreign-core/lapack.lisp
+++ b/src/foreign-core/lapack.lisp
@@ -1588,7 +1588,7 @@
(work (* :double-float) :workspace-output)
(lwork :integer :input)
(info :integer :output))
-
+
(def-fortran-routine dpotrf :void
"
SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
diff --git a/src/lapack/getrf.lisp b/src/lapack/getrf.lisp
index 0875520..d82ede5 100644
--- a/src/lapack/getrf.lisp
+++ b/src/lapack/getrf.lisp
@@ -33,43 +33,51 @@
(assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class)
`(defun ,func-name (A ipiv)
(declare (type ,matrix-class A)
- (type permutation-action ipiv))
- (mlet*
- (((maj-A ld-A fop-A) (blas-matrix-compatible-p A :n) :type (symbol index-type nil)))
- (assert maj-A nil 'tensor-not-consecutive-store)
- (multiple-value-bind (new-a new-ipiv info)
+ (type permutation-pivot-flip ipiv))
+ (mlet* (((maj-A ld-A fop-A) (blas-matrix-compatible-p A :n) :type (symbol index-type nil)))
+ (assert maj-A nil 'tensor-store-not-consecutive)
+ (multiple-value-bind (new-A new-ipiv info)
(,lapack-func
(nrows A) (ncols A) (store A)
ld-A (repr ipiv) 0)
- (declare (ignore new-a new-ipiv))
+ (declare (ignore new-A new-ipiv))
;;Convert from 1-based indexing to 0-based indexing, and fix
;;other Fortran-ic quirks
- (assert (= info 0) nil 'invalid-arguments)
- (format t "~a~%" ipiv)
- (let-typed ((ipiv-repr (repr ipiv) :type perrepr-vector))
- (loop for i of-type fixnum from 0 below (length ipiv-repr)
- do (decf (aref ipiv-repr i)))
- (loop for i of-type fixnum from 0 below (length ipiv-repr)
- do (let ((val (aref ipiv-repr i)))
- (setf (aref ipiv-repr val) i))))
- (values A ipiv info))))))
+ (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info))
+ (let-typed ((pidv (repr ipiv) :type perrepr-vector))
+ (very-quickly
+ (loop for i from 0 below (length pidv)
+ do (decf (aref pidv i)))))
+ (if (eq maj-A :row-major)
+ ;;Crout's decomposition
+ (values A (list :decomposition-type :|U_ii=1| :column-permutation ipiv))
+ ;;Dolittle's decomposition
+ (values A (list :decomposition-type :|L_ii=1| :row-permutation ipiv))))))))
(generate-typed-getrf! real-typed-getrf! (real-matrix dgetrf))
-
-
+(generate-typed-getrf! complex-typed-getrf! (complex-matrix zgetrf))
+
+#+nil
+(let ((A (make-real-tensor '((1 2)
+ (3 4))))
+ (idiv (make-pidx (perv 0 1))))
+ (real-typed-getrf! A idiv))
+
+
(defgeneric getrf! (A ipiv)
(:documentation
"
Syntax
======
- (GETRF a)
+ (GETRF! a)
Purpose
=======
Given an NxM matrix A, compute its LU factorization using
- partial pivoting, row interchanges:
+ partial pivoting, row or column interchanges:
- A = P * L * U
+ A = P * L * U (if A is row-major ordered)
+ A = L * U * P' (if A is col-major ordered)
where:
@@ -97,33 +105,50 @@
[3] INFO = T: successful
i: U(i,i) is exactly zero.
")
- (:method :before ((A standard-matrix) (ipiv permutation-action))
+ (:method :before ((A standard-matrix) (ipiv permutation-pivot-flip))
(assert (>= (group-rank ipiv) (idx-min (dimensions A))) nil 'invalid-value
:given (group-rank ipiv) :expected '(>= (group-rank ipiv) (idx-min (dimensions A))))))
-(defmethod getrf! ((a real-matrix) (ipiv permutation-action))
- (let* ((n (nrows a))
- (m (ncols a))
- (ipiv #+:pre-allocate-workspaces
- (or ipiv *ipiv*)
- #-:pre-allocate-workspaces
- (or ipiv (make-array (min n m) :element-type '(unsigned-byte 32)))))
-
- (declare (type fixnum n m))
- (multiple-value-bind (new-a new-ipiv info)
- (dgetrf n ;; M
- m ;; N
- (store a) ;; A
- n ;; LDA
- ipiv ;; IPIV
- 0) ;; INFO
- (declare (ignore new-a new-ipiv))
- (values a ipiv (if (zerop info)
- t
- info)))))
-
-
-
+(defmethod getrf! ((A real-matrix) (ipiv permutation-pivot-flip))
+ (let* ((copy? (not (consecutive-store-p A)))
+ (cp-A (if copy? (copy A) A))
+ (ret (multiple-value-list (real-typed-getrf! cp-A ipiv))))
+ (when copy?
+ (copy! (first ret) A)
+ (rplaca ret A))
+ (values-list ret)))
+
+(defmethod getrf! ((A complex-matrix) (ipiv permutation-pivot-flip))
+ (let* ((copy? (not (consecutive-store-p A)))
+ (cp-A (if copy? (copy A) A))
+ (ret (multiple-value-list (complex-typed-getrf! cp-A ipiv))))
+ (when copy?
+ (copy! (first ret) A)
+ (rplaca ret A))
+ (values-list ret)))
+
+(defun permute-idx (A arg perm)
+ (declare (type standard-matrix A)
+ (type permutation-pivot-flip perm))
+ (let* ((idiv (repr 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)))
+ (loop for i from 0 below (length idiv)
+ do (progn
+ (unless (= i (aref idiv i))
+ (setf (head ttwo) (* (aref idiv i) argstd))
+ (swap! tone ttwo))
+ (incf (head tone) argstd))))))
+ A)
+
+(defun split-lu (A op-info)
+ (declare (type standard-matrix A))
+ (destructuring-bind (&key decomposition-type row-permutation col-permutation) op-info
+ (assert (member decomposition-type '(:|U_ii=1| :|L_ii=1|)) nil 'invalid-arguments :message "Bad decomposition-type")
+
+;;
(defgeneric lu (a &key with-l with-u with-p)
(:documentation
"
@@ -147,31 +172,7 @@
"))
-
-(defmethod getrf! ((a complex-matrix) &optional ipiv)
- (let* ((n (nrows a))
- (m (ncols a))
- (ipiv #+:pre-allocate-workspaces
- (or ipiv *ipiv*)
- #-:pre-allocate-workspaces
- (or ipiv (make-array (min n m) :element-type '(unsigned-byte 32)))))
-
- (declare (type fixnum n m))
- (multiple-value-bind (new-a new-ipiv info)
- (zgetrf n ;; M
- m ;; N
- (store a) ;; A
- n ;; LDA
- ipiv ;; IPIV
- 0) ;; INFO
- (declare (ignore new-a new-ipiv))
- (values a ipiv (if (zerop info)
- t
- info)))))
-
-
(defmethod lu ((a standard-matrix) &key (with-l t) (with-u t) (with-p t))
-
(multiple-value-bind (lu ipiv info)
(getrf! (copy a))
(declare (ignore info))
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index 8cdb3d0..24e087f 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -105,39 +105,24 @@
(generate-typed-num-copy! complex-typed-num-copy!
(complex-tensor zcopy *complex-l1-fcall-lb*))
;;---------------------------------------------------------------;;
+;;Generic function defined in src;base;generic-copy.lisp
-(defgeneric copy! (from-tensor to-tensor)
- (:documentation
- "
- Syntax
- ======
- (COPY! x y)
-
- Purpose
- =======
- Copies the contents of the tensor X to
- the tensor Y, returns Y.
-
- X,Y must have the same dimensions, and
- ergo the same number of elements.
-
- Furthermore, X may be a scalar, in which
- case Y is filled with X.
-
+(defmethod copy! :before ((x standard-tensor) (y standard-tensor))
+ "
The contents of X must be coercable to
the type of Y. For example,
a COMPLEX-MATRIX cannot be copied to a
- REAL-MATRIX but the converse is possible.
-")
- (:method :before ((x standard-tensor) (y standard-tensor))
- (unless (idx= (dimensions x) (dimensions y))
- (error 'tensor-dimension-mismatch)))
- (:method ((x standard-tensor) (y standard-tensor))
- (mod-dotimes (idx (dimensions x))
- do (setf (tensor-ref y idx) (tensor-ref x idx)))
- y)
- (:method ((x complex-tensor) (y real-tensor))
- (error 'coercion-error :from 'complex-tensor :to 'real-tensor)))
+ REAL-MATRIX but the converse is possible."
+ (assert (idx= (dimensions x) (dimensions y)) nil
+ 'tensor-dimension-mismatch))
+
+(defmethod copy! ((x standard-tensor) (y standard-tensor))
+ (mod-dotimes (idx (dimensions x))
+ do (setf (tensor-ref y idx) (tensor-ref x idx)))
+ y)
+
+(defmethod copy! ((x complex-tensor) (y real-tensor))
+ (error 'coercion-error :from 'complex-tensor :to 'real-tensor))
(defmethod copy! ((x real-tensor) (y real-tensor))
(real-typed-copy! x y))
@@ -165,32 +150,7 @@
(defmethod copy! ((x number) (y complex-tensor))
(complex-typed-num-copy! (coerce-complex x) y))
-;;
-(defgeneric copy (tensor)
- (:documentation
- "
- Syntax
- ======
- (COPY x)
-
- Purpose
- =======
- Return a copy of the tensor X"))
-
-(defmethod copy ((tensor real-tensor))
- (let* ((ret (apply #'make-real-tensor (idx->list (dimensions tensor)))))
- (declare (type real-tensor ret))
- (copy! tensor ret)))
-
-(defmethod copy ((tensor complex-tensor))
- (let* ((ret (apply #'make-complex-tensor (idx->list (dimensions tensor)))))
- (declare (type complex-tensor ret))
- (copy! tensor ret)))
-
-(defmethod copy ((tensor number))
- tensor)
-
-;;
+;; Copy between a Lisp array and a tensor
(defun convert-to-lisp-array (tensor)
"
Syntax
@@ -208,7 +168,92 @@
:element-type (if-ret (getf (get-tensor-class-optimization (class-name (class-of tensor))) :element-type)
(error 'tensor-cannot-find-optimization :tensor-class (class-name (class-of tensor)))))))
(declare (type index-store-vector dims))
+ (let ((lst (make-list (rank tensor))))
+ (very-quickly
+ (mod-dotimes (idx dims)
+ do (setf (apply #'aref ret (idx->list! idx lst)) (tensor-ref tensor idx))))
+ ret)))
+
+(defmethod copy! :before ((x standard-tensor) (y array))
+ (assert (subtypep (element-type x)
+ (array-element-type y))
+ nil 'invalid-type
+ :given (element-type x)
+ :expected (array-element-type y))
+ (assert (and
+ (= (rank x) (array-rank y))
+ (reduce #'(lambda (x y) (and x y))
+ (mapcar #'= (idx->list (dimensions x)) (array-dimensions y))))
+ nil 'dimension-mismatch))
+
+(defmethod copy! ((x real-tensor) (y array))
+ (let-typed ((sto-x (store x) :type real-store-vector)
+ (lst (make-list (rank x)) :type cons))
+ (very-quickly
+ (mod-dotimes (idx (dimensions x))
+ with (linear-sums
+ (of-x (strides x) (head x)))
+ do (setf (apply #'aref y (idx->list! idx lst))
+ (aref sto-x of-x)))))
+ y)
+
+(defmethod copy! ((x complex-tensor) (y array))
+ (let-typed ((sto-x (store x) :type complex-store-vector)
+ (lst (make-list (rank x)) :type cons))
(very-quickly
- (mod-dotimes (idx dims)
- do (setf (apply #'aref ret (idx->list idx)) (tensor-ref tensor idx))))
- ret))
+ (mod-dotimes (idx (dimensions x))
+ with (linear-sums
+ |
