Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch matlisp-cffi updated. 2012-02-23-69-g9bfeec0
|
From: Akshay S. <ak...@us...> - 2012-03-20 18:35:45
|
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, matlisp-cffi has been updated
via 9bfeec0a8b2e5604b2ce6b7ad6be62c3fd3f09c4 (commit)
from fd41f88aefad9d87a8c9183f946ac14c3b564de8 (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 9bfeec0a8b2e5604b2ce6b7ad6be62c3fd3f09c4
Author: Akshay Srinivasan <aks...@gm...>
Date: Wed Mar 21 00:02:07 2012 +0530
-> Added classes for handling sub-matrices
diff --git a/matlisp.asd b/matlisp.asd
index 7898782..4a84895 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -145,7 +145,6 @@
"matlisp-lapack-wrappers")
:components ((:file "compat")
(:file "help")
- (:file "diag")
(:file "special")
(:file "reader")
(:file "trans")
diff --git a/packages.lisp b/packages.lisp
index 7b53ddc..e4fe2da 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -295,41 +295,44 @@
(defpackage :matlisp
(:use :common-lisp :fortran-ffi-accessors :blas :lapack :dfftpack :quadpack :matlisp-lib :utilities)
- (:shadow "REAL")
- (:export
- "*PRINT-MATRIX*"
- "AXPY!"
- "AXPY"
- "BLAS-COPYABLE-P"
- "COL-VECTOR-P"
- "COMPLEX-COERCE"
- "COMPLEX-MATRIX"
- "COMPLEX-MATRIX-ARRAY-TYPE"
- "COMPLEX-MATRIX-ELEMENT-TYPE"
- "COMPLEX-MATRIX-STORE-TYPE"
- "COPY!"
- "COPY"
- "CONVERT-TO-LISP-ARRAY"
- "CTRANSPOSE"
- "DIAG"
+ (:shadow #:real)
+ (:export #:*print-matrix*
+ ;;Level 1 BLAS
+ #:axpy! #:axpy
+ #:copy! #:copy
+ #:scal! #:scal
+ ;;Level 2 BLAS
+ #:gemv! #:gemv
+ ;;Level 3 BLAS
+ #:gemm! #:gemm
+ ;;Fortran stuff
+ #:blas-copyable-p #:blas-matrix-compatible-p
+ #:fortran-op #:fortran-nop #:fortran-snop
+ ;;Standard-matrix
+ #:standard-matrix
+ #:nrows #:ncols #:number-of-elements
+ #:head #:row-stride #:col-stride
+ #:store #:store-size
+ ;;Generic functions on standard-matrix
+ #:fill-matrix
+ #:ctranspose! #:ctranspose #:transpose #:transpose!
+ #:row-or-col-vector-p #:row-vector-p #:col-vector-p
+ #:row #:col #:diag #:sub-matrix
+ ;;Real-double-matrix
+ #:real-matrix #:real-matrix-element-type #:real-matrix-store-type
+ ;;Complex-double-matrix
+ #:complex-matrix #:complex-matrix-element-type #:complex-matrix-store-type #:complex-coerce #:complex-double-float
+ #:mrealpart #:mimagpart
+ ;;
+ "CONVERT-TO-LISP-ARRAY"
"DOT"
"EIG"
"EYE"
"FFT"
"FFT"
- "FILL-MATRIX"
- "FLOAT-MATRIX"
- "FLOAT-MATRIX-ARRAY-TYPE"
- "FLOAT-MATRIX-ELEMENT-TYPE"
- "FORTRAN-COMPLEX-MATRIX-INDEXING"
- "FORTRAN-MATRIX-INDEXING"
"GEEV"
"GELSY!"
"GELSY"
- #:gemv!
- #:gemv
- "GEMM!"
- "GEMM"
"GESV!"
"GESV"
"GETRF!"
@@ -337,7 +340,6 @@
"GETRS!"
"HELP"
"IFFT"
- "IMAG"
"JOIN"
"LOAD-BLAS-&-LAPACK-BINARIES"
"LOAD-BLAS-&-LAPACK-LIBRARIES"
@@ -375,8 +377,6 @@
"MASINH"
"MATAN"
"MATANH"
- "MATLISP-HERALD"
- "MATLISP-VERSION"
"MATRIX-REF"
"MCOS"
"MCOSH"
@@ -391,11 +391,6 @@
"MTANH"
"NCOLS"
"NORM"
- "NROWS"
- "NUMBER-OF-COLS"
- "NUMBER-OF-ELEMENTS"
- "NUMBER-OF-ELEMS"
- "NUMBER-OF-ROWS"
"ONES"
"PRINT-ELEMENT"
"QR"
@@ -405,28 +400,18 @@
"POTRS!"
"RAND"
"REAL"
- "REAL-MATRIX"
- "REAL-MATRIX-ELEMENT-TYPE"
- "REAL-MATRIX-STORE-TYPE"
"RESHAPE!"
"RESHAPE"
- "ROW-OR-COL-VECTOR-P"
- "ROW-VECTOR-P"
"SAVE-MATLISP"
- "SCAL!"
- "SCAL"
"SEQ"
"SET-M*!-SWAP-SIZE"
"SIZE"
"SQUARE-MATRIX-P"
- "STANDARD-MATRIX"
+ "STORE-INDEXING"
"SUM"
"SVD"
"SWAP!"
"TR"
- "TRANSPOSE"
- "TRANSPOSE!"
- "VEC"
"UNLOAD-BLAS-&-LAPACK-LIBRARIES"
"ZEROS"
;; From Quadpack
diff --git a/src/axpy.lisp b/src/axpy.lisp
index 69a7005..662cfd6 100644
--- a/src/axpy.lisp
+++ b/src/axpy.lisp
@@ -134,11 +134,11 @@ don't know how to coerce COMPLEX to REAL"))
(generate-typed-axpy!-func complex-double-axpy!-typed complex-double-float complex-matrix-store-type complex-matrix blas:zaxpy)
(defmethod axpy! ((alpha cl:real) (x real-matrix) (y complex-matrix))
- (real-double-axpy!-typed (coerce alpha 'double-float) x (realpart! y)))
+ (real-double-axpy!-typed (coerce alpha 'double-float) x (mrealpart y)))
(defmethod axpy! ((alpha complex) (x real-matrix) (y complex-matrix))
- (real-double-axpy!-typed (coerce (realpart alpha) 'double-float) x (realpart! y))
- (real-double-axpy!-typed (coerce (imagpart alpha) 'double-float) x (imagpart! y)))
+ (real-double-axpy!-typed (coerce (realpart alpha) 'double-float) x (mrealpart y))
+ (real-double-axpy!-typed (coerce (imagpart alpha) 'double-float) x (mimagpart y)))
(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix))
(complex-double-axpy!-typed (complex-coerce alpha) x y))
diff --git a/src/blas.lisp b/src/blas.lisp
index e7f76d9..1e34345 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -56,23 +56,9 @@
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-#+nil
-(defpackage "BLAS"
-#+:cmu (:use "COMMON-LISP" "ALIEN" "C-CALL" "FORTRAN-FFI-ACCESSORS")
-#+:sbcl (:use "COMMON-LISP" "SB-ALIEN" "SB-C" "FORTRAN-FFI-ACCESSORS")
-#+:allegro (:use "COMMON-LISP" "FOREIGN-FUNCTIONS" "FORTRAN-FFI-ACCESSORS")
- (:export
-"IDAMAX" "DASUM" "DDOT" "DNRM2"
-"DROT" "DSCAL" "DSWAP" "DCOPY" "DAXPY"
-"DCABS1" "DZASUM" "DZNRM2" "IZAMAX"
-"ZDSCAL" "ZSCAL" "ZSWAP" "ZCOPY" "ZAXPY" "ZDOTC" "ZDOTU"
-"DGEMV" "DSYMV" "DTRMV" "DTRSV" "DGER" "DSYR" "DSYR2"
-"ZGEMV" "ZHEMV" "ZTRMV" "ZTRSV" "ZGERC" "ZGERU" "ZHER2"
-"DGEMM" "DSYRK" "DSYR2K" "DTRMM" "DTRSM"
-"ZGEMM" "ZTRMM" "ZTRSM" "ZHERK" "ZHER2K" ))
-
(in-package "BLAS")
+;;
(def-fortran-routine daxpy :void
"
Syntax
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 704f823..3ad8280 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -43,18 +43,18 @@
:type (complex-matrix-store-type *)))
(:documentation "A class of matrices with complex elements."))
+(defclass sub-complex-matrix (complex-matrix)
+ ((parent-matrix
+ :initarg :parent
+ :accessor parent
+ :type complex-matrix))
+ (:documentation "A class of matrices with complex elements."))
+
;;
(defmethod initialize-instance ((matrix complex-matrix) &rest initargs)
- (setf (store-size matrix) (length (get-arg :store initargs)))
+ (setf (store-size matrix) (/ (length (get-arg :store initargs)) 2))
(call-next-method))
-(defmethod initialize-instance :after ((matrix complex-matrix) &rest initargs)
- (declare (ignore initargs))
- (let ((ss (store-size matrix)))
- (declare (type fixnum ss))
- (unless (>= ss (* 2 (number-of-elements matrix)))
- (error "Store is not large enough to hold the matrix."))))
-
;;
(defmethod matrix-ref-1d ((matrix complex-matrix) (idx fixnum))
(let ((store (store matrix)))
@@ -69,6 +69,75 @@
(aref store (+ 1 (* 2 idx))) (imagpart coerced-value))))
;;
+(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)))
+
+;;
+(defmethod sub-matrix ((matrix complex-matrix) (origin list) (dim list))
+ (destructuring-bind (o-i o-j) origin
+ (destructuring-bind (nr-s nc-s) dim
+ (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 *))))
+ (unless (and (< -1 o-i (+ o-j nr-s) nr) (< -1 o-j (+ o-j nc-s) nc))
+ (error "Bad index and/or size.
+Cannot create a sub-matrix of size (~a ~a) starting at (~a ~a)" nr-s nc-s o-i o-j))
+ (make-instance 'sub-complex-matrix
+ :nrows nr-s :ncols nc-s
+ :store st
+ :head (store-indexing o-i o-j hd rs cs)
+ :row-stride rs :col-stride cs)))))
+
+;;
+(defmethod row ((matrix complex-matrix) (i fixnum))
+ (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 *))))
+ (unless (< -1 i nr)
+ (error "Index ~a is outside the valid range for the given matrix." i))
+ (make-instance 'sub-complex-matrix
+ :nrows 1 :ncols nc
+ :store st
+ :head (store-indexing i 0 hd rs cs)
+ :row-stride rs :col-stride cs)))
+
+;;
+(defmethod col ((matrix complex-matrix) (j fixnum))
+ (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 *))))
+ (unless (< -1 j nc)
+ (error "Index ~a is outside the valid range for the given matrix." j))
+ (make-instance 'sub-complex-matrix
+ :nrows nr :ncols 1
+ :store st
+ :head (store-indexing 0 j hd rs cs)
+ :row-stride rs :col-stride cs)))
+
+;;
+(defmethod diag ((matrix complex-matrix) &optional (d 0))
+ (declare (type fixnum d))
+ (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 *)))
+ ((f-i f-j) (if (< d 0)
+ (values (- d) 0)
+ (values 0 d))
+ :type (fixnum fixnum)))
+ (unless (and (< -1 f-i nr) (< -1 f-j nc))
+ (error "Index ~a is outside the valid range for the given matrix." d))
+ (let ((d-s (min (- nr f-i) (- nc f-j))))
+ (declare (type fixnum d-s))
+ (make-instance 'sub-complex-matrix
+ :nrows 1 :ncols d-s
+ :store st
+ :head (store-indexing f-i f-j hd rs cs)
+ :row-stride 1 :col-stride (+ rs cs)))))
+
+;;
(declaim (inline allocate-complex-store))
(defun allocate-complex-store (size)
(make-array (* 2 size) :element-type 'complex-matrix-element-type
@@ -291,32 +360,28 @@
;;
-(defun realpart! (mat)
- (cond
- ((typep mat 'real-matrix) mat)
- ((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat)
- :nrows (nrows mat) :ncols (ncols mat)
- :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
- :head (* 2 (head mat))))))
-
-(defun imagpart! (mat)
- (cond
- ((typep mat 'real-matrix) nil)
- ((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat)
- :nrows (nrows mat) :ncols (ncols mat)
- :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
- :head (+ 1 (* 2 (head mat)))))))
-
-
-(defun conjugate! (mat)
- (cond
- ((typep mat 'real-matrix) nil)
- ((typep mat 'complex-matrix) (progn
- (transpose! (scal! -1d0 (imagpart! mat)))
- mat))))
-
-(defmacro with-conjugate! (matlst &rest body)
- `(progn
- ,@(mapcar #'(lambda (mat) `(conjugate! ,mat)) matlst)
- ,@body
- ,@(mapcar #'(lambda (mat) `(conjugate! ,mat)) matlst)))
\ No newline at end of file
+(defun mrealpart (mat)
+ (typecase mat
+ (real-matrix mat)
+ (complex-matrix (make-instance 'sub-real-matrix
+ :parent mat :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (* 2 (head mat))))
+ (number (cl:realpart mat))))
+
+(defun mimagpart (mat)
+ (typecase mat
+ (real-matrix nil)
+ (complex-matrix (make-instance 'sub-real-matrix
+ :parent mat :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (+ 1 (* 2 (head mat)))))
+ (number (cl:imagpart mat))))
+
+(defun mconjugate! (mat)
+ (typecase mat
+ (real-matrix mat)
+ (complex-matrix (scal! -1d0 (mimagpart mat))))
+ mat)
\ No newline at end of file
diff --git a/src/copy.lisp b/src/copy.lisp
index 8ce687d..74b1ca0 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -219,8 +219,8 @@ don't know how to coerce a COMPLEX to a REAL"))
(complex-double-copy!-typed x y))
(defmethod copy! ((x real-matrix) (y complex-matrix))
- (real-double-copy!-typed x (realpart! y))
- (scal! 0d0 (imagpart! y))
+ (real-double-copy!-typed x (mrealpart y))
+ (scal! 0d0 (mimagpart y))
y)
(defmethod copy! ((x number) (y complex-matrix))
diff --git a/src/diag.lisp b/src/diag.lisp
deleted file mode 100644
index 23ee370..0000000
--- a/src/diag.lisp
+++ /dev/null
@@ -1,235 +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: diag.lisp,v 1.8 2011/01/25 18:36:56 rtoy Exp $
-;;;
-;;; $Log: diag.lisp,v $
-;;; Revision 1.8 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.7.2.1 2011/01/25 18:16:53 rtoy
-;;; Use cl:real instead of real.
-;;;
-;;; Revision 1.7 2004/05/24 16:34:22 rtoy
-;;; More SBCL support from Robert Sedgewick. The previous SBCL support
-;;; was incomplete.
-;;;
-;;; Revision 1.6 2002/07/29 00:29:36 rtoy
-;;; Don't use *1x1-real-array*
-;;;
-;;; Revision 1.5 2001/10/29 16:21:02 rtoy
-;;; (setf diag) was broken on CMUCL. Use the Allegro version.
-;;; From M. Koerber.
-;;;
-;;; 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.
-;;;
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(in-package "MATLISP")
-
-(defgeneric diag (matrix)
- (:documentation
- "
- Syntax
- ======
- (DIAG x)
-
- Purpose
- =======
- Given the matrix X, returns the diagonal elements of X as a
- column vector if X is a matrix. If X is a vector
- returns a new matrix whose diagonal is X.
-
- Settable
- ========
- (SETF (DIAG x) y) works as follows.
-
- If Y is a scalar then the diagonal elements of X are assigned to Y.
- If Y is a vector then the diagonal elements of X are assigned to
- the elements of Y.
- If Y is a matrix then the diagonal elements of X are assigned to
- the diagonal elements of Y.
-
- The dimensions of X,Y need not match. In this case, the maximum
- assignable elements are considered.
-
- Returns X.
-"))
-
-(defmethod diag ((mat number))
- mat)
-
-(defmethod diag ((mat real-matrix))
- (if (row-or-col-vector-p mat)
- (let* ((nxm (number-of-elements mat))
- (result (make-real-matrix-dim nxm nxm)))
- (declare (type fixnum nxm))
- (dcopy nxm (store mat) 1 (store result) (1+ nxm))
- result)
- (let* ((n (nrows mat))
- (m (ncols mat))
- (p (min m n))
- (result (make-real-matrix-dim p 1)))
- (declare (type fixnum n m p))
- (dcopy p (store mat) (1+ n) (store result) 1)
- result)))
-
-(let ((diag-element (make-array 1 :element-type 'real-matrix-element-type)))
- (defmethod (setf diag) ((new-diag #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (mat real-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (p (min m n)))
- (declare (type fixnum n m p))
-
- (setf (aref diag-element 0) new-diag)
- (dcopy p diag-element 0 (store mat) (1+ n))
- mat)))
-
-(defmethod (setf diag) ((new-diag cl:real) (mat real-matrix))
- (setf (diag mat) (coerce new-diag 'real-matrix-element-type)))
-
-(defmethod (setf diag) ((new-diag complex) (mat real-matrix))
- (error "cannot set the diagonal of ~a to ~a, don't know how to
-coerce COMPLEX to REAL"
- mat
- new-diag))
-
-(defmethod (setf diag) ((new-diag real-matrix) (mat real-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (n-new (nrows new-diag))
- (m-new (ncols new-diag))
- (nxm-new (number-of-elements new-diag)))
- (declare (type fixnum n m n-new m-new nxm-new))
-
- (if (row-or-col-vector-p new-diag)
- (dcopy (min n m nxm-new) (store new-diag) 1 (store mat) (1+ n))
- (dcopy (min n m n-new m-new) (store new-diag) (1+ n-new) (store mat) (1+ n)))
-
- mat))
-
-(defmethod (setf diag) ((new-diag complex-matrix) (mat real-matrix))
- (error "cannot assign the COMPLEX matrix ~a to the diagonal of the REAL matrix ~a,
-don't know how to coerce COMPLEX to REAL"
- new-diag
- mat))
-
-(defmethod diag ((mat complex-matrix))
- (if (row-or-col-vector-p mat)
- (let* ((nxm (number-of-elements mat))
- (result (make-complex-matrix-dim nxm nxm)))
- (declare (type fixnum nxm))
- (zcopy nxm (store mat) 1 (store result) (1+ nxm))
- result)
- (let* ((n (nrows mat))
- (m (ncols mat))
- (p (min m n))
- (result (make-complex-matrix-dim p 1)))
-
- (declare (type fixnum n m p))
- (zcopy p (store mat) (1+ n) (store result) 1)
- result)))
-
-
-(defmethod (setf diag) ((new-diag complex-matrix) (mat complex-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (n-new (nrows new-diag))
- (m-new (ncols new-diag))
- (nxm-new (number-of-elements new-diag)))
- (declare (type fixnum n m n-new m-new nxm-new))
- (if (row-or-col-vector-p new-diag)
- (zcopy (min n m nxm-new) (store new-diag) 1 (store mat) (1+ n))
- (zcopy (min n m n-new m-new) (store new-diag) (1+ n-new) (store mat) (1+ n)))
- mat))
-
-
-(defmethod (setf diag) ((new-diag real-matrix) (mat complex-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (n-new (nrows new-diag))
- (m-new (ncols new-diag))
- (nxm-new (number-of-elements new-diag)))
- (declare (type fixnum n m n-new m-new nxm-new))
-
- (if (row-or-col-vector-p new-diag)
- (dotimes (i (min n m nxm-new))
- (declare (type fixnum i))
- (setf (matrix-ref mat (+ i (* n i)))
- (matrix-ref new-diag i)))
- (dotimes (i (min n m n-new m-new))
- (declare (type fixnum i))
- (setf (matrix-ref mat (+ i (* n i)))
- (matrix-ref new-diag (+ i (* n-new i))))))
-
- mat))
-
-(defmethod (setf diag) ((new-diag #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (mat complex-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (p (min n m)))
- (declare (type fixnum n m p))
-
- #-(or cmu sbcl) (setf new-diag (complex-coerce new-diag))
-
- (zcopy p new-diag 0 (store mat) (1+ n))
- mat))
-
-(defmethod (setf diag) ((new-diag number) (mat complex-matrix))
- (setf (diag mat) (complex-coerce new-diag)))
-
-
-
-
-
diff --git a/src/gemm.lisp b/src/gemm.lisp
index 419be54..e4d227c 100644
--- a/src/gemm.lisp
+++ b/src/gemm.lisp
@@ -248,7 +248,7 @@
(defmethod gemm! ((alpha number) (a real-matrix) (b real-matrix)
(beta complex) (c complex-matrix)
&optional (job :nn))
- (let ((r-c (realpart! c))
+ (let ((r-c (mrealpart c))
(c-be (complex-coerce beta)))
(declare (type real-matrix c)
(type complex-double-float c-al))
@@ -266,10 +266,10 @@
(defmethod gemm! ((alpha cl:real) (a real-matrix) (b complex-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-b (realpart! b))
- (i-b (imagpart! b))
- (r-c (realpart! c))
- (i-c (imagpart! c))
+ (let ((r-b (mrealpart b))
+ (i-b (mimagpart b))
+ (r-c (mrealpart c))
+ (i-c (mimagpart c))
(r-al (coerce alpha 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type real-matrix r-b i-b r-c i-c)
@@ -280,10 +280,10 @@
(defmethod gemm! ((alpha complex) (a real-matrix) (b complex-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-b (realpart! b))
- (i-b (imagpart! b))
- (r-c (realpart! c))
- (i-c (imagpart! c))
+ (let ((r-b (mrealpart b))
+ (i-b (mimagpart b))
+ (r-c (mrealpart c))
+ (i-c (mimagpart c))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
@@ -306,10 +306,10 @@
(defmethod gemm! ((alpha cl:real) (a complex-matrix) (b real-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-a (realpart! a))
- (i-a (imagpart! a))
- (r-c (realpart! c))
- (i-c (imagpart! c))
+ (let ((r-a (mrealpart a))
+ (i-a (mimagpart a))
+ (r-c (mrealpart c))
+ (i-c (mimagpart c))
(r-al (coerce alpha 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type real-matrix r-a i-a r-c i-c)
@@ -320,10 +320,10 @@
(defmethod gemm! ((alpha complex) (a complex-matrix) (b real-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-a (realpart! a))
- (i-a (imagpart! a))
- (r-c (realpart! c))
- (i-c (imagpart! c))
+ (let ((r-a (mrealpart a))
+ (i-a (mimagpart a))
+ (r-c (mrealpart c))
+ (i-c (mimagpart c))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
diff --git a/src/gemv.lisp b/src/gemv.lisp
index 88da9df..dbe8ce8 100644
--- a/src/gemv.lisp
+++ b/src/gemv.lisp
@@ -115,7 +115,7 @@
;
(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
(beta complex) (y complex-matrix) &optional (job :n))
- (let ((r-y (realpart! y)))
+ (let ((r-y (mrealpart y)))
(declare (type real-matrix r-y))
;; y <- \beta * y
(scal! (complex-coerce beta) y)
@@ -133,12 +133,12 @@
(beta cl:real) (y complex-matrix) &optional (job :n))
(let ((r-be (coerce beta 'double-float))
(r-al (coerce alpha 'double-float))
- (r-y (realpart! y)))
+ (r-y (mrealpart y)))
(declare (type double-float r-be r-al)
(type real-matrix r-y))
;; y <- \beta * y
(scal! r-be y)
- ;; (realpart! y) <- (realpart! y) + \alpha * A o x
+ ;; (mrealpart y) <- (mrealpart y) + \alpha * A o x
(real-double-gemv!-typed r-al A x 1d0 r-y job))
y)
@@ -147,13 +147,13 @@
(let ((r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float))
- (r-y (realpart! y))
- (i-y (imagpart! y)))
+ (r-y (mrealpart y))
+ (i-y (mimagpart y)))
(declare (type double-float r-al i-al r-be)
(type real-matrix r-y i-y))
- ;; (realpart! y) <- \beta * (realpart! y) + (realpart \alpha) . A o x
+ ;; (mrealpart y) <- \beta * (mrealpart y) + (realpart \alpha) . A o x
(real-double-gemv!-typed r-al A x r-be r-y job)
- ;; (imagpart! y) <- \beta * (imagpart! y) + (imagpart \alpha) . A o x
+ ;; (mimagpart y) <- \beta * (mimagpart y) + (imagpart \alpha) . A o x
(real-double-gemv!-typed i-al A x r-be i-y job))
y)
@@ -167,26 +167,26 @@
(defmethod gemv! ((alpha cl:real) (A real-matrix) (x complex-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-x (realpart! x))
- (i-x (imagpart! x))
- (r-y (realpart! y))
- (i-y (imagpart! y))
+ (let ((r-x (mrealpart x))
+ (i-x (mimagpart x))
+ (r-y (mrealpart y))
+ (i-y (mimagpart y))
(r-al (coerce (realpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type double-float r-al r-be)
(type real-matrix r-x i-x r-y i-y))
- ;; (realpart! y) <- \beta * (realpart! y) + \alpha . A o (realpart! x)
+ ;; (mrealpart y) <- \beta * (mrealpart y) + \alpha . A o (mrealpart x)
(real-double-gemv!-typed r-al A r-x r-be r-y job)
- ;; (imagpart! y) <- \beta * (imagpart! y) + \alpha . A o (realpart! x)
+ ;; (mimagpart y) <- \beta * (mimagpart y) + \alpha . A o (mrealpart x)
(real-double-gemv!-typed r-al A i-x r-be i-y job))
y)
(defmethod gemv! ((alpha complex) (A real-matrix) (x complex-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-x (realpart! x))
- (i-x (imagpart! x))
- (r-y (realpart! y))
- (i-y (imagpart! y))
+ (let ((r-x (mrealpart x))
+ (i-x (mimagpart x))
+ (r-y (mrealpart y))
+ (i-y (mimagpart y))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
@@ -209,26 +209,26 @@
(defmethod gemv! ((alpha cl:real) (A complex-matrix) (x real-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-A (realpart! A))
- (i-A (imagpart! A))
- (r-y (realpart! y))
- (i-y (imagpart! y))
+ (let ((r-A (mrealpart A))
+ (i-A (mimagpart A))
+ (r-y (mrealpart y))
+ (i-y (mimagpart y))
(r-al (coerce (realpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type double-float r-al r-be)
(type real-matrix r-A i-A r-y i-y))
- ;; (realpart! y) <- \beta * (realpart! y) + \alpha . A o (realpart! x)
+ ;; (mrealpart y) <- \beta * (mrealpart y) + \alpha . A o (mrealpart x)
(real-double-gemv!-typed r-al r-A x r-be r-y job)
- ;; (imagpart! y) <- \beta * (imagpart! y) + \alpha . A o (realpart! x)
+ ;; (mimagpart y) <- \beta * (mimagpart y) + \alpha . A o (mrealpart x)
(real-double-gemv!-typed r-al i-A x r-be i-y job))
y)
(defmethod gemv! ((alpha complex) (A complex-matrix) (x real-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-A (realpart! A))
- (i-A (imagpart! A))
- (r-y (realpart! y))
- (i-y (imagpart! y))
+ (let ((r-A (mrealpart A))
+ (i-A (mimagpart A))
+ (r-y (mrealpart y))
+ (i-y (mimagpart y))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
diff --git a/src/real-matrix.lisp b/src/real-matrix.lisp
index 7cb08a7..b6b90a2 100644
--- a/src/real-matrix.lisp
+++ b/src/real-matrix.lisp
@@ -18,18 +18,18 @@
:type (real-matrix-store-type *)))
(:documentation "A class of matrices with real elements."))
+(defclass sub-real-matrix (real-matrix)
+ ((parent-matrix
+ :initarg :parent
+ :accessor parent
+ :type real-matrix))
+ (:documentation "A class of matrices with real elements."))
+
;;
(defmethod initialize-instance ((matrix real-matrix) &rest initargs)
(setf (store-size matrix) (length (get-arg :store initargs)))
(call-next-method))
-(defmethod initialize-instance :after ((matrix real-matrix) &rest initargs)
- (declare (ignore initargs))
- (let ((ss (store-size matrix)))
- (declare (type fixnum ss))
- (unless (>= ss (number-of-elements matrix))
- (error "Store is not large enough to hold the matrix."))))
-
;;
(defmethod matrix-ref-1d ((matrix real-matrix) (idx fixnum))
(let ((store (store matrix)))
@@ -56,6 +56,75 @@
don't know how to coerce COMPLEX to REAL"))
;;
+(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 sub-matrix ((matrix real-matrix) (origin list) (dim list))
+ (destructuring-bind (o-i o-j) origin
+ (destructuring-bind (nr-s nc-s) dim
+ (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 *))))
+ (unless (and (< -1 o-i (+ o-j nr-s) nr) (< -1 o-j (+ o-j nc-s) nc))
+ (error "Bad index and/or size.
+Cannot create a sub-matrix of size (~a ~a) starting at (~a ~a)" nr-s nc-s o-i o-j))
+ (make-instance 'sub-real-matrix
+ :nrows nr-s :ncols nc-s
+ :store st
+ :head (store-indexing o-i o-j hd rs cs)
+ :row-stride rs :col-stride cs)))))
+
+;;
+(defmethod row ((matrix real-matrix) (i fixnum))
+ (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 *))))
+ (unless (< -1 i nr)
+ (error "Index ~a is outside the valid range for the given matrix." i))
+ (make-instance 'sub-real-matrix
+ :nrows 1 :ncols nc
+ :store st
+ :head (store-indexing i 0 hd rs cs)
+ :row-stride rs :col-stride cs)))
+
+;;
+(defmethod col ((matrix real-matrix) (j fixnum))
+ (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 *))))
+ (unless (< -1 j nc)
+ (error "Index ~a is outside the valid range for the given matrix." j))
+ (make-instance 'sub-real-matrix
+ :nrows nr :ncols 1
+ :store st
+ :head (store-indexing 0 j hd rs cs)
+ :row-stride rs :col-stride cs)))
+
+;;
+(defmethod diag ((matrix real-matrix) &optional (d 0))
+ (declare (type fixnum d))
+ (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 *)))
+ ((f-i f-j) (if (< d 0)
+ (values (- d) 0)
+ (values 0 d))
+ :type (fixnum fixnum)))
+ (unless (and (< -1 f-i nr) (< -1 f-j nc))
+ (error "Index ~a is outside the valid range for the given matrix." d))
+ (let ((d-s (min (- nr f-i) (- nc f-j))))
+ (declare (type fixnum d-s))
+ (make-instance 'sub-real-matrix
+ :nrows 1 :ncols d-s
+ :store st
+ :head (store-indexing f-i f-j hd rs cs)
+ :row-stride 1 :col-stride (+ rs cs)))))
+
+;;
(defun make-real-matrix-dim (n m &key (fill 0.0d0) (order :row-major))
"
Syntax
@@ -252,4 +321,4 @@ don't know how to coerce COMPLEX to REAL"))
"cannot make a ~A x ~A matrix" n m)
(make-real-matrix-dim n m)))
(t
- (error "require 1 or 2 arguments to make a matrix")))))
\ No newline at end of file
+ (error "require 1 or 2 arguments to make a matrix")))))
diff --git a/src/special.lisp b/src/special.lisp
index e00193f..f53ca0c 100644
--- a/src/special.lisp
+++ b/src/special.lisp
@@ -88,8 +88,7 @@
(error "the number of rows (~d) and columns (~d) must be positive integers" n m))
(let ((result (make-real-matrix-dim n m)))
- (setf (aref *1x1-real-array* 0) 1.0d0)
- (dcopy (min n m) *1x1-real-array* 0 (store result) (1+ n))
+ (setf (diag result) 1d0)
result))
(defun zeros (n &optional (m n))
@@ -160,5 +159,4 @@
(make-instance 'real-matrix
:nrows n :ncols m
- :row-stride m :col-stride 1 :store store))))
-
+ :store store))))
\ No newline at end of file
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index 4c918b8..c1ef52e 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -22,6 +22,44 @@ integer storage. Default INITIAL-ELEMENT = 0."
(declare (type (and fixnum (integer 0)) row col head row-stride col-stride))
(the fixnum (+ head (the fixnum (* row row-stride)) (the fixnum (* col col-stride)))))
+(defun blas-copyable-p (matrix)
+ (declare (optimize (safety 0) (speed 3))
+ (type (or real-matrix complex-matrix) matrix))
+ (mlet* ((nr (nrows matrix) :type fixnum)
+ (nc (ncols matrix) :type fixnum)
+ (rs (row-stride matrix) :type fixnum)
+ (cs (col-stride matrix) :type fixnum)
+ (ne (number-of-elements matrix) :type fixnum))
+ (cond
+ ((or (= nc 1) (= cs (* nr rs))) (values t rs ne))
+ ((or (= nr 1) (= rs (* nc cs))) (values t cs ne))
+ (t (values nil -1 -1)))))
+
+(defun blas-matrix-compatible-p (matrix &optional (op :n))
+ (declare (optimize (safety 0) (speed 3))
+ (type (or real-matrix complex-matrix) matrix))
+ (mlet* (((rs cs) (slot-values matrix '(row-stride col-stride))
+ :type (fixnum fixnum)))
+ (cond
+ ((= cs 1) (values :row-major rs (fortran-nop op)))
+ ((= rs 1) (values :col-major cs (fortran-op op)))
+ ;;Lets not confound lisp's type declaration.
+ (t (values nil -1 "?")))))
+
+(declaim (inline fortran-op))
+(defun fortran-op (op)
+ (ecase op (:n "N") (:t "T")))
+
+(declaim (inline fortran-nop))
+(defun fortran-nop (op)
+ (ecase op (:t "N") (:n "T")))
+
+(defun fortran-snop (sop)
+ (cond
+ ((string= sop "N") "T")
+ ((string= sop "T") "N")
+ (t (error "Unrecognised fortran-op."))))
+
;;
(defclass standard-matrix ()
((number-of-rows
@@ -76,30 +114,29 @@ that way."))
;;
(defmethod initialize-instance :after ((matrix standard-matrix) &rest initargs)
(declare (ignore initargs))
- (let* ((n (nrows matrix))
- (m (ncols matrix))
- (h (head matrix))
- (ss (store-size matrix))
- (nxm (* n m)))
- (declare (type fixnum n m h nxm))
- ;;Row-ordered by default.
- (unless (and (slot-boundp matrix 'row-stride) (slot-boundp matrix 'col-stride))
- (setf (row-stride matrix) m)
- (setf (col-stride matrix) 1))
- (let ((rs (row-stride matrix))
- (cs (row-stride matrix)))
- (declare (type fixnum rs cs))
- ;;Error checking is good if we use foreign-pointers as store types.
- (cond
- ((<= n 0) (error "Number of rows must be > 0. Initialized with ~A." n))
- ((<= m 0) (error "Number of columns must be > 0. Initialized with ~A." m))
- ;;
- ((< h 0) (error "Head of the store must be >= 0. Initialized with ~A." h))
- ((< rs 0) (error "Row-stride of the store must be >= 0. Initialized with ~A." rs))
- ((< cs 0) (error "Column-stride of the store must be >= 0. Initialized with ~A." cs))
- ((<= ss 0) (error "Store-size must be > 0. Initialized with ~A." ss))))
+ (mlet*
+ (((nr nc hd ss) (slot-values matrix '(number-of-rows number-of-cols head store-size))
+ :type (fixnum fixnum fixnum fixnum)))
+ ;;Row-ordered by default.
+ (unless (and (slot-boundp matrix 'row-stride) (slot-boundp matrix 'col-stride))
+ (setf (row-stride matrix) nc)
+ (setf (col-stride matrix) 1))
+ (let* ((rs (row-stride matrix))
+ (cs (col-stride matrix))
+ (l-idx (store-indexing (- nr 1) (- nc 1) hd rs cs)))
+ (declare (type fixnum rs cs))
+ ;;Error checking is good if we use foreign-pointers as store types.
+ (cond
+ ((<= nr 0) (error "Number of rows must be > 0. Initialized with ~A." nr))
+ ((<= nc 0) (error "Number of columns must be > 0. Initialized with ~A." nc))
+ ;;
+ ((< hd 0) (error "Head of the store must be >= 0. Initialized with ~A." hd))
+ ((< rs 0) (error "Row-stride of the store must be >= 0. Initialized with ~A." rs))
+ ((< cs 0) (error "Column-stride of the store must be >= 0. Initialized with ~A." cs))
+ ((<= ss l-idx) (error "Store is not large enough to hold the matrix.
+Initialized with ~A, but the largest possible index is ~A." ss l-idx))))
;;
- (setf (number-of-elements matrix) nxm)))
+ (setf (number-of-elements matrix) (* nr nc))))
;;
(defmacro matrix-ref (matrix row &optional col)
@@ -265,16 +302,17 @@ matrix and a number"))
(format stream "~%")))
;;
-(defun transpose-i! (matrix)
+(defun transpose! (matrix)
"
Syntax
======
- (transpose-i! matrix)
+ (transpose! matrix)
Purpose
=======
Exchange row and column strides so that effectively
- the matrix is transposed in place (without much effort).
+ the matrix is destructively transposed in place
+ (without much effort).
"
(cond
((typep matrix 'standard-matrix)
@@ -285,126 +323,144 @@ matrix and a number"))
((typep 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)))
+
;;
-(defun transpose! (matrix)
+(defgeneric transpose (matrix)
+ (:documentation
"
Syntax
======
- (transpose! matrix)
+ (transpose matrix)
Purpose
=======
Create a new matrix object which represents the transpose of the
- the given matrix, but shares the store with matrix.
-"
- (cond
- ((typep matrix 'standard-matrix)
- (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
- (make-instance (class-of matrix)
- :nrows nc :ncols nr
- :store (store matrix)
- :head hd
- :row-stride cs :col-stride rs)))
- ((typep matrix 'number) matrix)
- (t (error "Don't know how to take the transpose of ~A." matrix))))
+ the given matrix.
-;;
-(defmacro with-transpose! (matlst &rest body)
- `(progn
- ,@(mapcar #'(lambda (mat) `(transpose-i! ,mat)) matlst)
- ,@body
- ,@(mapcar #'(lambda (mat) `(transpose-i! ,mat)) matlst)))
+ Store is shared with \"matrix\".
-;;
-(defun sub! (matrix i j nrows ncols)
- (declare (type standard-matrix matrix)
- (type fixnum i j nrows ncols))
- (let ((hd (head matrix))
- (nr (nrows matrix))
- (nc (ncols matrix))
- (rs (row-stride matrix))
- (cs (col-stride matrix)))
- (declare (type fixnum hd nr nc rs cs))
- (unless (and (< -1 i (+ i nrows) nr) (< -1 j (+ j ncols) nc))
- (error "Bad index and/or size.
-Cannot create a sub-matrix of size (~a ~a) starting at (~a ~a)" nrows ncols i j))
- (make-instance (class-of matrix)
- :nrows nrows :ncols ncols
- :store (store matrix)
- :head (store-indexing i j hd rs cs)
- :row-stride rs :col-stride cs)))
-
-(defun (setf sub!) (mat-b mat-a i j nrows ncols)
- (copy! mat-b (sub! mat-a i j nrows ncols)))
+ Settable
+ ========
+ (setf (transpose matrix) value)
-;;
-(defun row! (matrix i)
- (declare (type standard-matrix matrix)
- (type fixnum i))
- (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
- (unless (< -1 i nr)
- (error "Index ~a is outside the valid range for the given matrix." i))
- (make-instance (class-of matrix)
- :nrows 1 :ncols nc
- :store (store matrix)
- :head (store-indexing i 0 hd rs cs)
- :row-stride rs :col-stride cs)))
-
-(defun (setf row!) (mat-b mat-a i)
- (copy! mat-b (row! mat-a i)))
+ is basically the same as
+
+ (copy! value (transpose matrix))
+"))
+
+(defun (setf transpose) (value matrix)
+ (copy! value (transpose matrix)))
+
+(defmethod transpose ((matrix number))
+ matrix)
;;
-(defun col! (matrix j)
- (declare (type standard-matrix matrix)
- (type fixnum j))
- (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
- (unless (< -1 j nc)
- (error "Index ~a is outside the valid range for the given matrix." j))
- (make-instance (class-of matrix)
- :nrows nr :ncols 1
- :store (store matrix)
- :head (store-indexing 0 j hd rs cs)
- :row-stride rs :col-stride cs)))
-
-(defun (setf col!) (mat-b mat-a j)
- (copy! mat-b (col! mat-a j)))
+(defgeneric sub-matrix (matrix origin dim)
+ (:documentation
+"
+ Syntax
+ ======
+ (sub-matrix matrix origin dimensions)
+
+ Purpose
+ =======
+ Create a block sub-matrix of \"matrix\" starting at \"origin\"
+ of dimension \"dim\", sharing the store.
+
+ origin, dim are lists with two elements.
+
+ Store is shared with \"matrix\"
+
+ Settable
+ ========
+ (setf (sub-matrix matrix origin dim) value)
+
+ is basically the same as
+
+ (copy! value (sub-matrix matrix origin dim))
+"))
+
+(defun (setf sub-matrix) (value matrix origin dim)
+ (copy! value (sub-matrix matrix origin dim)))
;;
-(defun blas-copyable-p (matrix)
- (declare (optimize (safety 0) (speed 3))
- (type (or real-matrix complex-matrix) matrix))
- (mlet* ((nr (nrows matrix) :type fixnum)
- (nc (ncols matrix) :type fixnum)
- (rs (row-stride matrix) :type fixnum)
- (cs (col-stride matrix) :type fixnum)
- (ne (number-of-elements matrix) :type fixnum))
- (cond
- ((or (= nc 1) (= cs (* nr rs))) (values t rs ne))
- ((or (= nr 1) (= rs (* nc cs))) (values t cs ne))
- (t (values nil -1 -1)))))
+(defgeneric row (matrix i)
+ (:documentation
+"
+ Syntax
+ ======
+ (row matrix i)
+
+ Purpose
+ =======
+ Returns the i'th row of the matrix.
+ Store is shared with \"matrix\".
+
+ Settable
+ ========
+ (setf (row matrix i) value)
+
+ is basically the same as
+
+ (copy! value (row matrix i))
+"))
+
+(defun (setf row) (value matrix i)
+ (copy! value (row matrix i)))
;;
-(defun blas-matrix-compatible-p (matrix &optional (op :n))
- (declare (optimize (safety 0) (speed 3))
- (type (or real-matrix complex-matrix) matrix))
- (mlet* (((rs cs) (slot-values matrix '(row-stride col-stride))
- :type (fixnum fixnum)))
- (cond
- ((= cs 1) (values :row-major rs (fortran-nop op)))
- ((= rs 1) (values :col-major cs (fortran-op op)))
- ;;Lets not confound lisp's type declaration.
- (t (values nil -1 "?")))))
+(defgeneric col (matrix j)
+ (:documentation
+"
+ Syntax
+ ======
+ (col matrix j)
+
+ Purpose
+ =======
+ Returns the j'th column of the matrix.
+ Store is shared with \"matrix\".
+
+ Settable
+ ========
+ (setf (col matrix j) value)
+
+ is basically the same as
+
+ (copy! value (col matrix j))
+"))
+
+(defun (setf col) (value matrix j)
+ (copy! value (col matrix j)))
+
;;
-(declaim (inline fortran-op))
-(defun fortran-op (op)
- (ecase op (:n "N") (:t "T")))
+(defgeneric diag (matrix &optional d)
+ (:documentation
+"
+ Syntax
+ ======
+ (diag matrix &optional (d 0))
-(declaim (inline fortran-nop))
-(defun fortran-nop (op)
- (ecase op (:t "N") (:n "T")))
+ Purpose
+ =======
+ Returns a row-vector representing the d'th diagonal of the matrix.
+ [a_{ij} : j - i = d]
-(defun fortran-snop (sop)
- (cond
- ((string= sop "N") "T")
- ((string= sop "T") "N")
- (t (error "Unrecognised fortran-op."))))
\ No newline at end of file
+ Store is shared with \"matrix\".
+
+ Settable
+ ========
+ (setf (diag matrix d) value)
+
+ is basically the same as
+
+ (copy! value (diag matrix d))
+"))
+
+(defun (setf diag) (value matrix &optional (d 0))
+ (copy! value (diag matrix d)))
\ No newline at end of file
diff --git a/src/trace.lisp b/src/trace.lisp
index 7bdeb3f..8490613 100644
--- a/src/trace.lisp
+++ b/src/trace.lisp
@@ -48,14 +48,12 @@
(in-package "MATLISP")
-#+nil (export '(tr))
-
-(defgeneric tr (a)
+(defgeneric mtrace (a)
(:documentation
"
Syntax
======
- (TR a)
+ (MTRACE a)
Purpose
=======
@@ -65,11 +63,11 @@
Notes
=====
- (TR a) is the same as (SUM (DIAG a))
+ (MTRACE a) is the same as (SUM (DIAG a))
"))
-(defmethod tr ((x number))
+(defmethod mtrace ((x number))
x)
-(defmethod tr ((a standard-matrix))
- (sum (diag a)))
+(defmethod mtrace ((a standard-matrix))
+ (sum (diag a 0)))
\ No newline at end of file
-----------------------------------------------------------------------
Summary of changes:
matlisp.asd | 1 -
packages.lisp | 77 +++++-------
src/axpy.lisp | 6 +-
src/blas.lisp | 16 +---
src/complex-matrix.lisp | 139 +++++++++++++++------
src/copy.lisp | 4 +-
src/diag.lisp | 235 ----------------------------------
src/gemm.lisp | 34 +++---
src/gemv.lisp | 54 ++++----
src/real-matrix.lisp | 85 +++++++++++--
src/special.lisp | 6 +-
src/standard-matrix.lisp | 314 +++++++++++++++++++++++++++-------------------
src/trace.lisp | 14 +--
13 files changed, 453 insertions(+), 532 deletions(-)
delete mode 100644 src/diag.lisp
hooks/post-receive
--
matlisp
|
