Menu
▾
▴
matlisp-commit
|
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))
+ ...
[truncated message content] |
Thanks for helping keep SourceForge clean.
X