Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-140-gcd30ca8
|
From: Akshay S. <ak...@us...> - 2012-07-24 05:48:35
|
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 cd30ca81e687388cf532e30e08f79b68cf56c325 (commit)
from 3d2b1c49901f857eff0b30ebecaeb251d35e1755 (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 cd30ca81e687388cf532e30e08f79b68cf56c325
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Jul 24 11:13:16 2012 +0530
o Added (:c, :h) to gemv!
o Moved all the BLAS tweakable parameters to a file
diff --git a/matlisp.asd b/matlisp.asd
index be38acc..fefa76d 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -100,7 +100,9 @@
(:file "blas-helpers"
:depends-on ("standard-tensor" "permutation"))
(:file "print"
- :depends-on ("standard-tensor"))))
+ :depends-on ("standard-tensor"))
+ ;;Probably not the right place, but should do.
+ (:file "tweakable")))
(:module "matlisp-classes"
:pathname "classes"
:depends-on ("matlisp-base")
@@ -122,7 +124,9 @@
(:file "dot"
:depends-on ("realimag"))
(:file "axpy"
- :depends-on ("copy" "scal"))))
+ :depends-on ("copy" "scal"))
+ (:file "trans"
+ :depends-on ("scal" "copy"))))
(:module "matlisp-level-2"
:pathname "level-2"
:depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1")
diff --git a/packages.lisp b/packages.lisp
index b295a79..ca65b2c 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -33,6 +33,7 @@
;;<conditon {accessors*}>
;;Generic errors
#:generic-error #:message
+ #:assumption-violated
#:invalid-type #:given #:expected
#:invalid-arguments
#:invalid-value #:given #:expected
diff --git a/src/classes/complex-tensor.lisp b/src/classes/complex-tensor.lisp
index ee4dd08..9692b49 100644
--- a/src/classes/complex-tensor.lisp
+++ b/src/classes/complex-tensor.lisp
@@ -23,12 +23,20 @@
(make-array (* 2 size) :element-type 'complex-base-type
:initial-element (coerce 0 'complex-base-type)))
-(definline coerce-complex (x)
+(definline coerce-complex-unforgiving (x)
(coerce x 'complex-type))
-(definline coerce-complex-base (x)
+(defun coerce-complex (x)
+ (restart-case (coerce-complex-unforgiving x)
+ (use-value (value) (coerce-complex value))))
+
+(definline coerce-complex-base-unforgiving (x)
(coerce x 'complex-base-type))
+(defun coerce-complex-base (x)
+ (restart-case (coerce-complex-base-unforgiving x)
+ (use-value (value) (coerce-complex-base value))))
+
;;
(defclass complex-tensor (standard-tensor)
((store
@@ -62,7 +70,7 @@ Cannot hold complex numbers."))
(tensor-store-defs (complex-tensor complex-type complex-base-type)
:store-allocator allocate-complex-store
- :coercer coerce-complex
+ :coercer coerce-complex-unforgiving
:reader
(lambda (tstore idx)
(complex (aref tstore (* 2 idx))
diff --git a/src/classes/real-tensor.lisp b/src/classes/real-tensor.lisp
index 7d360a6..6ece4df 100644
--- a/src/classes/real-tensor.lisp
+++ b/src/classes/real-tensor.lisp
@@ -13,9 +13,13 @@
"(allocate-real-store size [initial-element])
Allocates real storage. Default initial-element = 0d0.")
-(definline coerce-real (x)
+(definline coerce-real-unforgiving (x)
(coerce x 'real-type))
+(defun coerce-real (x)
+ (restart-case (coerce-real-unforgiving x)
+ (use-value (value) (coerce-real value))))
+
;;
(defclass real-tensor (standard-tensor)
((store
@@ -43,7 +47,7 @@ Allocates real storage. Default initial-element = 0d0.")
;;
(tensor-store-defs (real-tensor real-type real-type)
:store-allocator allocate-real-store
- :coercer coerce-real
+ :coercer coerce-real-unforgiving
:reader
(lambda (tstore idx)
(aref tstore idx))
@@ -68,5 +72,3 @@ Allocates real storage. Default initial-element = 0d0.")
(defmethod print-element ((tensor real-tensor)
element stream)
(format stream "~11,5,,,,,'Eg" element))
-
-
diff --git a/src/conditions.lisp b/src/conditions.lisp
index b384847..0cbcb67 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -24,6 +24,12 @@
(:method print-object ((c generic-error) stream)
(format stream (message c))))
+(defcondition assumption-violated (generic-error)
+ ()
+ (:method print-object ((c assumption-violated) stream)
+ (format stream "An assumption assumed when writing the software has been violated. Proceed with caution.")
+ (call-next-method)))
+
(defcondition invalid-type (generic-error)
((given-type :reader given :initarg :given)
(expected-type :reader expected :initarg :expected))
diff --git a/src/ffi/c-ffi.lisp b/src/ffi/c-ffi.lisp
index 7a6fba9..109a756 100644
--- a/src/ffi/c-ffi.lisp
+++ b/src/ffi/c-ffi.lisp
@@ -7,37 +7,47 @@
(real ,base-type)
(imag ,base-type)))
-(defccomplex c-complex-double :double)
-(defccomplex c-complex-float :float)
+(defccomplex %c.complex-double :double)R
+(defccomplex %c.complex-float :float)
;; Get the equivalent CFFI type.
;; If the type is an array, get the type of the array element type.
-(defun c->cffi-type (type)
+(defun %c.cffi-type (type)
"Convert the given Fortran FFI type into a type understood by CFFI."
(cond
- ((and (listp type) (eq (first type) '*))
- `(:pointer ,(c->cffi-type
- (case (second type)
- ;;CDR coding ?
- (:complex-single-float :single-float)
- (:complex-double-float :double-float)
- (t (second type))))))
+ ;; '* means arrays of a type, which isn't necessarily the same as pointer-to-type
+ ;; (* :complex-single-float) expands to (:pointer :float)
+ ;; (:pointer :complex-single-float) expands to (:pointer (:struct %c.complex-float))
+ ((consp type)
+ (cond
+ ((eq (first type) '*)
+ `(:pointer ,(%c.cffi-type
+ (case (second type)
+ ;;CDR coding ?
+ (:complex-single-float :single-float)
+ (:complex-double-float :double-float)
+ (t (second type))))))
+ ;;We assume you what you're doing, and
+ ;;CFFI knows the type.
+ ((eq (first type) :pointer)
+ type)
((callback-type-p type)
- `(:pointer ,(c->cffi-type :callback)))
+ `(:pointer ,(%c.cffi-type :callback)))
((eq type :complex-single-float)
- `(:struct c-complex-float))
+ `(:struct %c.complex-float))
((eq type :complex-double-float)
- `(:struct c-complex-double))
+ `(:struct %c.complex-double))
(t (case type
(:void :void)
(:integer :int)
(:long :long)
- (:single-float 'c-complex-float)
- (:double-float 'c-complex-double)
+ (:single-float :float)
+ (:double-float :double)
(:string :string)
+ (:character :char)
;; Pass a pointer to the function.
(:callback :void)
- ;;We assume the type is known to CFFI.
+ ;;We assume that the type is known to CFFI.
(t type)))))
;; Check if given type is a string
diff --git a/src/ffi/f77-ffi.lisp b/src/ffi/f77-ffi.lisp
index d68b04b..602eba0 100644
--- a/src/ffi/f77-ffi.lisp
+++ b/src/ffi/f77-ffi.lisp
@@ -50,9 +50,9 @@
`(:pointer ,(%f77.cffi-type :double-float)))
(t (ecase type
(:void :void)
- (:integer :int)
+ (:integer :int32)
(:character :char)
- (:long :long)
+ (:long :int64)
(:single-float :float)
(:double-float :double)
(:string :string)
diff --git a/src/ffi/ffi-cffi.lisp b/src/ffi/ffi-cffi.lisp
index 9597bc2..8429984 100644
--- a/src/ffi/ffi-cffi.lisp
+++ b/src/ffi/ffi-cffi.lisp
@@ -22,6 +22,10 @@
:string :character
:callback))
+(define-constant +ffi-array-types+
+ '(:single-float :double-float
+ :integer :long))
+
;; Separte the body of code into documentation and parameter lists.
(defun parse-doc-&-parameters (body &optional header footer)
(if (stringp (first body))
@@ -144,13 +148,9 @@
;;
(simple-array (signed-byte 64) *)
(simple-array (signed-byte 32) *)
- (simple-array (signed-byte 16) *)
- (simple-array (signed-byte 8) *)
;;
(simple-array (unsigned-byte 64) *)
(simple-array (unsigned-byte 32) *)
- (simple-array (unsigned-byte 16) *)
- (simple-array (unsigned-byte 8) *)
;;
cffi:foreign-pointer))
diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp
index e789d7a..ce870da 100644
--- a/src/level-1/axpy.lisp
+++ b/src/level-1/axpy.lisp
@@ -99,26 +99,19 @@
,(funcall (getf opt :value-writer) '(+ num-from val) 't-sto 't-of))))))))
to)))
-;;Tweakable
-(defparameter *real-axpy-fortran-call-lower-bound* 20000
- "If the size of the array is less than this parameter, the
- lisp version of axpy is called in order to avoid FFI overheads")
-(generate-typed-axpy! real-typed-axpy! (real-tensor
- daxpy
- *real-axpy-fortran-call-lower-bound*))
-(generate-typed-num-axpy! real-typed-num-axpy! (real-tensor
- daxpy
- *real-axpy-fortran-call-lower-bound*))
-;;Tweakable
-(defparameter *complex-axpy-fortran-call-lower-bound* 10000
- "If the size of the array is less than this parameter, the
- lisp version of axpy is called in order to avoid FFI overheads")
-(generate-typed-axpy! complex-typed-axpy! (complex-tensor
- zaxpy
- *complex-axpy-fortran-call-lower-bound*))
-(generate-typed-num-axpy! complex-typed-num-axpy! (complex-tensor
- zaxpy
- *complex-axpy-fortran-call-lower-bound*))
+;;Real
+(generate-typed-axpy! real-typed-axpy!
+ (real-tensor daxpy *real-l1-fcall-lb*))
+
+(generate-typed-num-axpy! real-typed-num-axpy!
+ (real-tensor daxpy *real-l1-fcall-lb*))
+
+;;Complex
+(generate-typed-axpy! complex-typed-axpy!
+ (complex-tensor zaxpy *complex-l1-fcall-lb*))
+
+(generate-typed-num-axpy! complex-typed-num-axpy!
+ (complex-tensor zaxpy *complex-l1-fcall-lb*))
;;---------------------------------------------------------------;;
(defgeneric axpy! (alpha x y)
diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp
index ba062de..8cdb3d0 100644
--- a/src/level-1/copy.lisp
+++ b/src/level-1/copy.lisp
@@ -91,40 +91,20 @@
to)))
-;;Tweakable
-(defparameter *real-copy-fortran-call-lower-bound* 20000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL running
- on x86-64 linux. A reasonable value would be something above 1000.")
-(generate-typed-copy! real-typed-copy! (real-tensor
- dcopy
- *real-copy-fortran-call-lower-bound*))
-(generate-typed-num-copy! real-typed-num-copy! (real-tensor
- dcopy
- *real-copy-fortran-call-lower-bound*))
-
-;;Tweakable
-(defparameter *complex-copy-fortran-call-lower-bound* 10000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL
- running on x86-64 linux. A reasonable value would be something
- above 1000.")
-
-(generate-typed-copy! complex-typed-copy! (complex-tensor
- zcopy
- *complex-copy-fortran-call-lower-bound*))
-(generate-typed-num-copy! complex-typed-num-copy! (complex-tensor
- zcopy
- *complex-copy-fortran-call-lower-bound*))
-;;---------------------------------------------------------------;;
+;;Real
+(generate-typed-copy! real-typed-copy!
+ (real-tensor dcopy *real-l1-fcall-lb*))
+
+(generate-typed-num-copy! real-typed-num-copy!
+ (real-tensor dcopy *real-l1-fcall-lb*))
-(defun test-copy (n r)
- (let ((x (make-real-tensor n)))
- (time (dotimes (i r)
- (copy! pi x)))
- t))
+;;Complex
+(generate-typed-copy! complex-typed-copy!
+ (complex-tensor zcopy *complex-l1-fcall-lb*))
+
+(generate-typed-num-copy! complex-typed-num-copy!
+ (complex-tensor zcopy *complex-l1-fcall-lb*))
+;;---------------------------------------------------------------;;
(defgeneric copy! (from-tensor to-tensor)
(:documentation
diff --git a/src/level-1/dot.lisp b/src/level-1/dot.lisp
index 60ad32a..15bf751 100644
--- a/src/level-1/dot.lisp
+++ b/src/level-1/dot.lisp
@@ -27,16 +27,11 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package #:matlisp)
-(defparameter *real-dot-fortran-call-lower-bound* 20000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL running
- on x86-64 linux. A reasonable value would be something above 1000.")
(defun real-typed-dot (x y conjugate-p)
(declare (type real-vector x y)
(ignore conjugate-p))
(let ((call-fortran? (> (number-of-elements x)
- *real-dot-fortran-call-lower-bound*)))
+ *real-l1-fcall-lb*)))
(cond
(call-fortran?
(ddot (number-of-elements x)
@@ -57,16 +52,10 @@
summing (* (aref sto-x of-x) (aref sto-y of-y)) into dot of-type real-type
finally (return dot))))))))
-
-(defparameter *complex-dot-fortran-call-lower-bound* 10000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL running
- on x86-64 linux. A reasonable value would be something above 1000.")
(defun complex-typed-dot (x y conjugate-p)
(declare (type complex-vector x y))
(let ((call-fortran? (> (number-of-elements x)
- *complex-dot-fortran-call-lower-bound*)))
+ *complex-l1-fcall-lb*)))
(cond
(call-fortran?
(if conjugate-p
diff --git a/src/level-1/scal.lisp b/src/level-1/scal.lisp
index c7604b1..5bdb664 100644
--- a/src/level-1/scal.lisp
+++ b/src/level-1/scal.lisp
@@ -50,21 +50,18 @@
,(funcall (getf opt :value-writer) 'scal-val 't-sto 't-of))))))))
to)))
-;; TODO: Maybe add zdscal support ? Don't think the difference between
-;; zdscal and zscal is significant, except for very large arrays.
-(defparameter *real-scal-fortran-call-lower-bound* 20000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL running
- on x86-64 linux. A reasonable value would be something above 1000.")
-(generate-typed-scal! real-typed-scal! (real-tensor dscal *real-scal-fortran-call-lower-bound*))
-
-(defparameter *complex-scal-fortran-call-lower-bound* 10000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL running
- on x86-64 linux. A reasonable value would be something above 1000.")
-(generate-typed-scal! complex-typed-scal! (complex-tensor zscal *complex-scal-fortran-call-lower-bound*))
+;;Real
+(generate-typed-scal! real-typed-scal!
+ (real-tensor dscal *real-l1-fcall-lb*))
+
+;;Complex
+(definline zordscal (nele alpha x incx &optional hd-x)
+ (if (zerop (imagpart alpha))
+ (zdscal nele (realpart alpha) x incx hd-x)
+ (zscal nele alpha x incx hd-x)))
+
+(generate-typed-scal! complex-typed-scal!
+ (complex-tensor zordscal *complex-l1-fcall-lb*))
;;---------------------------------------------------------------;;
(defgeneric scal! (alpha x)
diff --git a/src/level-1/swap.lisp b/src/level-1/swap.lisp
index 7f80c78..7204214 100644
--- a/src/level-1/swap.lisp
+++ b/src/level-1/swap.lisp
@@ -53,19 +53,11 @@
do ,(funcall (getf opt :swapper) 'f-sto 'f-of 't-sto 't-of)))))))
y)))
-(defparameter *real-swap-fortran-call-lower-bound* 20000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL running
- on x86-64 linux. A reasonable value would be something above 1000.")
-(generate-typed-swap! real-typed-swap! (real-tensor dswap *real-swap-fortran-call-lower-bound*))
+(generate-typed-swap! real-typed-swap!
+ (real-tensor dswap *real-l1-fcall-lb*))
-(defparameter *complex-scal-fortran-call-lower-bound* 10000
- "
- If the dimension of the arguments is less than this parameter,
- then the Lisp version of copy is used. Default set with SBCL running
- on x86-64 linux. A reasonable value would be something above 1000.")
-(generate-typed-swap! complex-typed-swap! (complex-tensor zswap *complex-scal-fortran-call-lower-bound*))
+(generate-typed-swap! complex-typed-swap!
+ (complex-tensor zswap *complex-l1-fcall-lb*))
;;---------------------------------------------------------------;;
(defgeneric swap! (x y)
diff --git a/src/level-1/trans.lisp b/src/level-1/trans.lisp
new file mode 100644
index 0000000..56bad3e
--- /dev/null
+++ b/src/level-1/trans.lisp
@@ -0,0 +1,210 @@
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
+;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;;
+;;; Copyright (c) 2000 The Regents of the University of California.
+;;; All rights reserved.
+;;;
+;;; Permission is hereby granted, without written agreement and without
+;;; license or royalty fees, to use, copy, modify, and distribute this
+;;; software and its documentation for any purpose, provided that the
+;;; above copyright notice and the following two paragraphs appear in all
+;;; copies of this software.
+;;;
+;;; IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+;;; FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
+;;; ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
+;;; THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
+;;; SUCH DAMAGE.
+;;;
+;;; THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
+;;; INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+;;; MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
+;;; PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
+;;; CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
+;;; ENHANCEMENTS, OR MODIFICATIONS.
+;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+(in-package #:matlisp)
+
+(defun transpose! (A &optional permutation)
+ "
+ Syntax
+ ======
+ (TRANSPOSE! a [permutation])
+
+ Purpose
+ =======
+ Exchange the arguments of the tensor in place. The default
+ is to swap the first and last arguments of the tensor.
+
+ Settable
+ ========
+ (setf (TRANSPOSE! tensor permutation) value)
+
+ is basically the same as
+ (copy! value (TRANSPOSE! tensor permutation)).
+
+ NOTE: This will have side-effects even if copy! doesn't succeed."
+ (declare (type standard-tensor a))
+ (if permutation
+ (progn
+ (permute! (strides A) permutation)
+ (permute! (dimensions A) permutation))
+ (let-typed ((rnk (rank A) :type index-type)
+ (dim-A (dimensions A) :type index-store-vector)
+ (strd-A (strides A) :type index-store-vector))
+ (rotatef (aref dim-A (1- rnk)) (aref dim-A 0))
+ (rotatef (aref strd-A (1- rnk)) (aref strd-A 0))))
+ A)
+
+(defun (setf transpose!) (value A &optional permutation)
+ (copy! value (transpose! A permutation)))
+
+(defun transpose~ (A &optional permutation)
+ "
+ Syntax
+ ======
+ (TRANSPOSE~ a permutation)
+
+ Purpose
+ =======
+ Like TRANSPOSE!, but the permuted strides and dimensions are part of
+ a new tensor object instead, the store being shared with the given
+ tensor.
+
+ Settable
+ ========
+ (setf (TRANSPOSE~ tensor permutation) value)
+
+ is basically the same as
+ (copy! value (TRANSPOSE~ tensor permutation))"
+ (declare (type standard-tensor A))
+ (let ((displaced (make-instance (class-of A) :store (store A)
+ :dimensions (copy-seq (dimensions A))
+ :strides (copy-seq (strides A))
+ :parent-tensor A)))
+ (transpose! displaced permutation)))
+
+(defun (setf transpose~) (value A &optional permutation)
+ (declare (type standard-tensor A))
+ (copy! value (transpose~ A permutation)))
+
+(definline transpose (A &optional permutation)
+ "
+ Syntax
+ ======
+ (TRANSPOSE~ a permutation)
+
+ Purpose
+ =======
+ Like TRANSPOSE!, but the permutation is applied on a copy of
+ the given tensor.
+
+ Settable
+ ========
+ (setf (TRANSPOSE tensor permutation) value)
+
+ is the same as (setf (transpose~ ..) ..)"
+ (declare (type standard-tensor A))
+ (copy (transpose~ A permutation)))
+
+(defun (setf transpose) (value A &optional permutation)
+ (declare (type standard-tensor A))
+ (copy! value (transpose~ A permutation)))
+
+;;This is a bit more complicated, now that we are no longer in S_2
+;;Computing the inverse permutation is trivial in the cyclic representation,
+;;but probably not worth the trouble for this ugly macro.
+#+nil
+(defmacro with-transpose! (matlst &rest body)
+ `(progn
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
+ ,@body
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
+
+
+;;
+(defun mconjugate! (A)
+ "
+ Syntax
+ ======
+ (mconjugate! A)
+
+ Purpose
+ =======
+ Destructively modifies A into its complex conjugate (not hermitian conjugate).
+
+ (tensor-imagpart~ A) <- (- (tensor-imagpart~ A)) "
+ (etypecase A
+ (real-tensor A)
+ (complex-tensor
+ (scal! -1 (tensor-imagpart~ A))
+ A)
+ (number (conjugate A))))
+
+(definline mconjugate (A)
+ "
+ Syntax
+ ======
+ (mconjugate A)
+
+ Purpose
+ =======
+ Like mconjugate!, but non-destructive."
+ (etypecase A
+ (standard-tensor (copy (mconjugate! A)))
+ (number (conjugate A))))
+
+;;
+(defun htranspose! (A &optional permutation)
+"
+ Syntax
+ ======
+ (HTRANSPOSE! A [permutation])
+
+ Purpose
+ =======
+ Hermitian transpose of A (destructive).
+"
+ (declare (type standard-tensor A))
+ (transpose! A permutation)
+ (when (typep A 'complex-tensor)
+ (mconjugate! A))
+ A)
+
+(definline ctranspose! (A &optional permutation)
+ "
+ Syntax
+ ======
+ (CTRANSPOSE! A [permutation])
+
+ Purpose
+ =======
+ Conjugate transpose of A (destructive).
+"
+ (htranspose! A permutation))
+
+(definline htranspose (A &optional permutation)
+"
+ Syntax
+ ======
+ (HTRANSPOSE A [permutation])
+
+ Purpose
+ =======
+ Like HTRANSPOSE!, but non-destructive."
+ (declare (type standard-tensor A))
+ (let ((result (copy A)))
+ (htranspose! result permutation)))
+
+(definline ctranspose (A &optional permutation)
+ "
+ Syntax
+ ======
+ (CTRANSPOSE A [permutation])
+
+ Purpose
+ =======
+ Like CTRANSPOSE!, but non-destructive."
+ (htranspose A permutation))
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index 7fffe5e..d4dfadd 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -59,34 +59,32 @@
`(+ (* alpha dotp) val) 'sto-y 'of-y))))))))))
y)))
-;;Tweakable
-(defparameter *real-gemv-fortran-call-lower-bound* 1000
- "
- If the maximum dimension in the MV is lower than this
- parameter, then the lisp code is used by default, instead of
- calling BLAS. Used to avoid the FFI overhead when calling
- MM with small matrices.
- Default set with SBCL on x86-64 linux. A reasonable value
- is something between 800 and 2000.")
-(generate-typed-gemv! real-typed-gemv! (real-matrix real-vector
- dgemv
- *real-gemv-fortran-call-lower-bound*))
-
-;;Tweakable
-(defparameter *complex-gemv-fortran-call-lower-bound* 600
- "
- If the maximum dimension in the MV is lower than this
- parameter, then the lisp code is used by default, instead of
- calling BLAS. Used to avoid the FFI overhead when calling
- MM with small matrices.
- Default set with SBCL on x86-64 linux. A reasonable value
- is something between 400 and 1000.")
-(generate-typed-gemv! complex-typed-gemv! (complex-matrix complex-vector
- zgemv
- *complex-gemv-fortran-call-lower-bound*))
+;;Real
+(generate-typed-gemv! real-base-typed-gemv!
+ (real-matrix real-vector dgemv *real-l2-fcall-lb*))
+
+(definline real-typed-gemv! (alpha A x beta y job)
+ (real-base-typed-gemv! alpha A x beta y (ecase job ((:n :t) job) (:h :t) (:c :n))))
+
+;;Complex
+(generate-typed-gemv! complex-base-typed-gemv!
+ (complex-matrix complex-vector zgemv *complex-l2-fcall-lb*))
+
+(defun complex-typed-gemv! (alpha A x beta y job)
+ (declare (type complex-matrix A)
+ (type complex-vector x y))
+ (if (member job '(:n :t))
+ (complex-base-typed-gemv! alpha A x beta y job)
+ ;;The CBLAS way.
+ (let ((cx (mconjugate x)))
+ (complex-base-typed-gemv! (cl:conjugate alpha) A cx
+ (cl:conjugate beta) (mconjugate! y) (ecase job (:h :t) (:c :n)))
+ (mconjugate! y))))
+
;;---------------------------------------------------------------;;
-;;Can't support "C" because the dual isn't supported by BLAS.
+;;Can't support "C" because its dual (complex-conjugate without transpose)
+;;isn't supported by BLAS (which'd be needed for row-major matrices).
(defgeneric gemv! (alpha A x beta y &optional job)
(:documentation
"
@@ -111,13 +109,15 @@
JOB Operation
---------------------------------------------------
:N (default) alpha * A * x + beta * y
- :T alpha * A'* x + beta * y
+ :T alpha * transpose(A)* x + beta * y
+ :C alpha * conjugate(A) * x + beta * y
+ :H alpha * transpose o conjugate(A) + beta * y
")
(:method :before ((alpha number) (A standard-matrix) (x standard-vector)
(beta number) (y standard-vector)
&optional (job :n))
- (assert (member job '(:n :t)) nil 'invalid-value
- :given job :expected `(member job '(:n :t))
+ (assert (member job '(:n :t :c :h)) nil 'invalid-value
+ :given job :expected `(member job '(:n :t :c :h))
:message "Inside gemv!")
(assert (not (eq x y)) nil 'invalid-arguments
:message "GEMV!: x and y cannot be the same vector")
@@ -143,7 +143,7 @@
(unless (= beta 1)
(complex-typed-scal! (coerce-complex beta) y))
(unless (= alpha 0)
- (if (complexp alpha)
+ (if (not (zerop (imagpart alpha)))
(let ((A.x (make-real-tensor (aref (dimensions y) 0)))
(vw-y (tensor-realpart~ y)))
(real-typed-gemv! (coerce-real 1) A x (coerce-real 0) A.x job)
@@ -204,6 +204,8 @@
---------------------------------------------------
:N (default) alpha * A * x + beta * y
:T alpha * A'* x + beta * y
+ :C alpha * conjugate(A) * x + beta * y
+ :H alpha * transpose o conjugate(A) + beta * y
"))
(defmethod gemv ((alpha number) (A standard-matrix) (x standard-vector)
@@ -227,5 +229,5 @@
(typep A 'complex-matrix) (typep x 'complex-vector))
#'make-complex-tensor
#'make-real-tensor)
- (list (ecase job (:n (nrows A)) (:t (ncols A)))))))
+ (list (ecase job ((:n :c) (nrows A)) ((:t :h) (ncols A)))))))
(gemv! alpha A x 0 result job)))
diff --git a/src/old/trans.lisp b/src/old/trans.lisp
deleted file mode 100644
index f20b2de..0000000
--- a/src/old/trans.lisp
+++ /dev/null
@@ -1,208 +0,0 @@
-;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Copyright (c) 2000 The Regents of the University of California.
-;;; All rights reserved.
-;;;
-;;; Permission is hereby granted, without written agreement and without
-;;; license or royalty fees, to use, copy, modify, and distribute this
-;;; software and its documentation for any purpose, provided that the
-;;; above copyright notice and the following two paragraphs appear in all
-;;; copies of this software.
-;;;
-;;; IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
-;;; FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
-;;; ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
-;;; THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
-;;; SUCH DAMAGE.
-;;;
-;;; THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
-;;; INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-;;; MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
-;;; PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
-;;; CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
-;;; ENHANCEMENTS, OR MODIFICATIONS.
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Originally written by Raymond Toy.
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; $Id: trans.lisp,v 1.6 2011/01/25 18:36:56 rtoy Exp $
-;;;
-;;; $Log: trans.lisp,v $
-;;; Revision 1.6 2011/01/25 18:36:56 rtoy
-;;; Merge changes from automake-snapshot-2011-01-25-1327 to get the new
-;;; automake build infrastructure.
-;;;
-;;; Revision 1.5.2.1 2011/01/25 18:16:53 rtoy
-;;; Use cl:real instead of real.
-;;;
-;;; Revision 1.5 2001/06/22 12:52:41 rtoy
-;;; Use ALLOCATE-REAL-STORE and ALLOCATE-COMPLEX-STORE to allocate space
-;;; instead of using the error-prone make-array.
-;;;
-;;; Revision 1.4 2000/07/11 18:02:03 simsek
-;;; o Added credits
-;;;
-;;; Revision 1.3 2000/07/11 02:11:56 simsek
-;;; o Added support for Allegro CL
-;;;
-;;; Revision 1.2 2000/05/08 17:19:18 rtoy
-;;; Changes to the STANDARD-MATRIX class:
-;;; o The slots N, M, and NXM have changed names.
-;;; o The accessors of these slots have changed:
-;;; NROWS, NCOLS, NUMBER-OF-ELEMENTS
-;;; The old names aren't available anymore.
-;;; o The initargs of these slots have changed:
-;;; :nrows, :ncols, :nels
-;;;
-;;; Revision 1.1 2000/04/14 00:11:12 simsek
-;;; o This file is adapted from obsolete files 'matrix-float.lisp'
-;;; 'matrix-complex.lisp' and 'matrix-extra.lisp'
-;;; o Initial revision.
-;;;
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-;; notes
-;; =====
-;; transposition is usually a redundant operation. For example, all BLAS/LAPACK
-;; operators take an extra argument asking whether the matrix is transposed,
-;; for small operations transposition doesn't make a difference, for repeated
-;; small operations -- it may, so you need to use this feature of the
-;; interfaced BLAS/LAPACK functions.
-;;
-;; also, the intent that TRANSPOSE creates a new matrix should be made clear,
-;; for example, taking the transpose of a row/column vector is easy, due to
-;; representation, but this will not create a new matrix.
-
-(in-package #:matlisp)
-
-(defun transpose! (matrix)
-"
- Syntax
- ======
- (TRANSPOSE! matrix)
-
- Purpose
- =======
- Exchange row and column strides so that effectively
- the matrix is destructively transposed in place
- (without much effort).
-"
- (typecase matrix
- (standard-matrix
- (progn
- (rotatef (nrows matrix) (ncols matrix))
- (rotatef (row-stride matrix) (col-stride matrix))
- matrix))
- (number matrix)
- (t (error "Don't know how to take the transpose of ~A." matrix))))
-
-(defmacro with-transpose! (matlst &rest body)
- `(progn
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
- ,@body
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
-
-;;
-(defgeneric transpose~ (matrix)
- (:documentation
-"
- Syntax
- ======
- (TRANSPOSE~ matrix)
-
- Purpose
- =======
- Create a new matrix object which represents the transpose of the
- the given matrix.
-
- Store is shared with \"matrix\".
-
- Settable
- ========
- (setf (TRANSPOSE~ matrix) value)
-
- is basically the same as
-
- (copy! value (TRANSPOSE~ matrix))
-"))
-
-(defun (setf transpose~) (value matrix)
- (copy! value (transpose~ matrix)))
-
-;;
-(defmethod transpose~ ((matrix number))
- matrix)
-
-(defmethod transpose~ ((matrix real-matrix))
- (mlet* (((hd nr nc rs cs st) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride store))
- :type (fixnum fixnum fixnum fixnum fixnum (real-matrix-store-type *))))
- (make-instance 'sub-real-matrix
- :nrows nc :ncols nr
- :store st
- :head hd
- :row-stride cs :col-stride rs
- :parent matrix)))
-
-(defmethod transpose~ ((matrix complex-matrix))
- (mlet* (((hd nr nc rs cs st) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride store))
- :type (fixnum fixnum fixnum fixnum fixnum (complex-matrix-store-type *))))
- (make-instance 'sub-complex-matrix
- :nrows nc :ncols nr
- :store st
- :head hd
- :row-stride cs :col-stride rs
- :parent matrix)))
-
-;;
-(declaim (inline transpose))
-(defun transpose (matrix)
-"
- Syntax
- ======
- (TRANSPOSE matrix)
-
- Purpose
- =======
- Creates a new matrix which is the transpose of MATRIX.
-"
- (copy (transpose~ matrix)))
-
-;;
-(defun ctranspose! (matrix)
-"
- Syntax
- ======
- (CTRANSPOSE! matrix)
-
- Purpose
- =======
- Exchange row and column strides so that effectively
- the matrix is destructively transposed in place
- (without much effort). Also scale the imagpart with -1,
- so that the end result is the Hermitian conjugate.
-"
- (transpose! matrix)
- (when (typep matrix 'complex-matrix)
- (scal! -1d0 (mimagpart~ matrix)))
- matrix)
-
-;;
-(defun ctranspose (matrix)
-"
- Syntax
- ======
- (CTRANSPOSE matrix)
-
- Purpose
- =======
- Returns a new matrix which is the conjugate transpose
- of MATRIX.
-"
- (let ((result (copy matrix)))
- (ctranspose! result)))
\ No newline at end of file
-----------------------------------------------------------------------
Summary of changes:
matlisp.asd | 8 +-
packages.lisp | 1 +
src/classes/complex-tensor.lisp | 14 ++-
src/classes/real-tensor.lisp | 10 +-
src/conditions.lisp | 6 +
src/ffi/c-ffi.lisp | 42 +++++---
src/ffi/f77-ffi.lisp | 4 +-
src/ffi/ffi-cffi.lisp | 8 +-
src/level-1/axpy.lisp | 33 +++----
src/level-1/copy.lisp | 44 +++------
src/level-1/dot.lisp | 15 +---
src/level-1/scal.lisp | 25 ++---
src/level-1/swap.lisp | 16 +---
src/level-1/trans.lisp | 210 +++++++++++++++++++++++++++++++++++++++
src/level-2/gemv.lisp | 64 ++++++------
src/old/trans.lisp | 208 --------------------------------------
16 files changed, 347 insertions(+), 361 deletions(-)
create mode 100644 src/level-1/trans.lisp
delete mode 100644 src/old/trans.lisp
hooks/post-receive
--
matlisp
|
