From: Akshay S. <ak...@us...> - 2012-07-14 05:16:32
|
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 d18665bf3b836e17d2ff75065b384b5ff07059e3 (commit) from 2b87e86f1392efee853a1807d7c9299fee1f7958 (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 d18665bf3b836e17d2ff75065b384b5ff07059e3 Author: Akshay Srinivasan <aks...@gm...> Date: Sat Jul 14 10:41:46 2012 +0530 Started work on LAPACK wrappers. diff --git a/matlisp.asd b/matlisp.asd index f5cb98e..e0e8283 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -118,7 +118,7 @@ (:file "dot" :depends-on ("realimag")) (:file "axpy" - :depends-on ("copy")))) + :depends-on ("copy" "scal")))) (:module "matlisp-level-2" :pathname "level-2" :depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1") @@ -126,7 +126,15 @@ (:module "matlisp-level-3" :pathname "level-3" :depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1") - :components ((:file "gemm"))))) + :components ((:file "gemm"))) + (:module "matlisp-lapack" + :pathname "lapack" + :depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3") + :components ((:file "gesv"))) + (:module "matlisp-sugar" + :pathname "sugar" + :depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3") + :components ((:file "mplusminus"))))) ;; (defclass f2cl-cl-source-file (asdf:cl-source-file) diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index d5f18df..017072b 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -300,6 +300,7 @@ (apply func-a (permute! (multiple-value-list (funcall func-b args)) perm)))) ;; +;;Optimize: pick different pivot. (defun idx-sort-permute (seq predicate) " (sort-permute seq predicate) diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 37e7cd8..8dcf4d7 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -31,6 +31,24 @@ (definline idxv (&rest contents) (make-index-store contents)) + +;; +(deftype integer4 () + '(signed-byte 32)) + +(deftype integer4-array (size) + `(simple-array integer4-array (,size))) + +(make-array-allocator allocate-integer4-store 'integer4 0 + " + Syntax + ====== + (ALLOCATE-INTEGER4-STORE SIZE [INITIAL-ELEMENT 0]) + + Purpose + ======= + Allocates integer4 (32-bits) storage.") + ;; (defclass standard-tensor () diff --git a/src/foreign-core/lapack.lisp b/src/foreign-core/lapack.lisp index 5d19674..9ce7a66 100644 --- a/src/foreign-core/lapack.lisp +++ b/src/foreign-core/lapack.lisp @@ -307,7 +307,7 @@ The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. - A (input) DOUBLE PRECISION array, dimension (LDA,N) + A (input/output) DOUBLE PRECISION array, dimension (LDA,N) The factors L and U from the factorization A = P*L*U as computed by DGETRF. diff --git a/src/old/gesv.lisp b/src/lapack/gesv.lisp similarity index 70% rename from src/old/gesv.lisp rename to src/lapack/gesv.lisp index 732bfda..e31c7b1 100644 --- a/src/old/gesv.lisp +++ b/src/lapack/gesv.lisp @@ -25,51 +25,10 @@ ;;; ENHANCEMENTS, OR MODIFICATIONS. ;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;; -;;; Originally written by Raymond Toy -;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;; -;;; $Id: gesv.lisp,v 1.4 2000/07/11 18:02:03 simsek Exp $ -;;; -;;; $Log: gesv.lisp,v $ -;;; 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") -#+:nil (use-package "BLAS") -#+:nil (use-package "LAPACK") -#+:nil (use-package "FORTRAN-FFI-ACCESSORS") +(in-package #:matlisp) -#+:nil (export '(gesv! - gesv)) - - -#+:pre-allocate-workspaces -(defvar *ipiv* (make-array *ipiv-size* :element-type '(unsigned-byte 32))) - -(defgeneric gesv! (a b &key ipiv) +(defgeneric gesv! (a b) (:documentation " Syntax @@ -108,50 +67,19 @@ used in the computation has been completed, but the factor U is exactly singular. Solution could not be computed. -")) - -(defgeneric gesv (a b) - (:documentation - " - Sytnax - ====== - (GESV a b) - - Purpose - ======= - Same as GESV! except that A,B are not overwritten. -")) - -(defmethod gesv! :before ((a standard-matrix) (b standard-matrix) &key ipiv) - (let ((n-a (nrows a)) - (m-a (ncols a)) - (n-b (nrows b))) - (if (not (= n-a m-a n-b)) - (error "dimensions of A,B given to GESV do not match")) - (if ipiv - (progn - (check-type ipiv (simple-array (unsigned-byte 32) (*))) - (if (< (length ipiv) n-a) - (error "argument IPIV given to GESV! must dimension >= N, -where NxN is the dimension of argument A given to GESV!")))))) - - -(defmethod gesv! ((a real-matrix) (b real-matrix) &key ipiv) - - (let* ((n (nrows a)) - (m (ncols b)) - (ipiv #+:pre-allocate-workspaces - (or ipiv *ipiv*) - #-:pre-allocate-workspaces - (or ipiv (make-array n :element-type '(unsigned-byte 32))))) - - (declare (type fixnum n m) - (type (simple-array (unsigned-byte 32) (*)) ipiv)) - - (multiple-value-bind (factors - ipiv - x - info) +") + (:before ((a standard-matrix) (b standard-matrix) &key ipiv) + (assert (= (nrows a) (ncols a) (nrows b)) nil 'tensor-dimension-mismatch))) + +(defmethod gesv! ((a real-matrix) (b real-matrix)) + + (let* ((nrc-a (nrows a)) + (nc-b (ncols b)) + (ipiv (make-integer4-store nrc-a))) + (declare (type fixnum nrc-a nc-b) + (type (integer4-array *) ipiv)) + + (multiple-value-bind (lu ipiv x info) (dgesv n m (store a) @@ -214,6 +142,20 @@ where NxN is the dimension of argument A given to GESV!")))))) (gesv! a b :ipiv ipiv))) + +(defgeneric gesv (a b) + (:documentation + " + Sytnax + ====== + (GESV a b) + + Purpose + ======= + Same as GESV! except that A,B are not overwritten. +")) + + (defmethod gesv :before ((a standard-matrix) (b standard-matrix)) (let ((n-a (nrows a)) (m-a (ncols a)) diff --git a/src/old/getrf.lisp b/src/lapack/getrf.lisp similarity index 100% rename from src/old/getrf.lisp rename to src/lapack/getrf.lisp diff --git a/src/old/getrs.lisp b/src/lapack/getrs.lisp similarity index 100% rename from src/old/getrs.lisp rename to src/lapack/getrs.lisp diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp index 0a3b9cc..e789d7a 100644 --- a/src/level-1/axpy.lisp +++ b/src/level-1/axpy.lisp @@ -212,3 +212,16 @@ (defmethod axpy ((alpha number) (x complex-tensor) (y complex-tensor)) (let ((ret (copy y))) (axpy! alpha x ret))) + +(defmethod axpy ((alpha number) (x (eql nil)) (y complex-tensor)) + (let ((ret (copy y))) + (axpy! alpha nil ret))) + +(defmethod axpy ((alpha number) (x (eql nil)) (y real-tensor)) + (let ((ret (if (complexp alpha) + (copy! y (apply #'make-complex-tensor (idx->list (dimensions y)))) + (copy y)))) + (axpy! alpha nil ret))) + +(defmethod axpy ((alpha number) (x standard-tensor) (y (eql nil))) + (scal alpha x)) diff --git a/src/old/mminus.lisp b/src/old/mminus.lisp deleted file mode 100644 index c769c32..0000000 --- a/src/old/mminus.lisp +++ /dev/null @@ -1,118 +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: mminus.lisp,v 1.4 2000/07/11 18:02:03 simsek Exp $ -;;; -;;; $Log: mminus.lisp,v $ -;;; 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") - -#+nil (use-package "BLAS") -#+nil (use-package "LAPACK") -#+nil (use-package "FORTRAN-FFI-ACCESSORS") - -#+nil (export '(m- - m.-)) - -(defgeneric m- (a b) - (:documentation - " - Syntax - ====== - (M- a b) - - Purpose - ======= - Create a new matrix which is the difference of A and B. - B may be a scalar, in which case the subtraction - is elementwise. -")) - -(defgeneric m.- (a b) - (:documentation - " - Syntax - ====== - (M- a b) - - Purpose - ======= - Same as M- -")) - -(defmethod m.- (a b) - (m- a b)) - -(defmethod m- :before ((a standard-matrix) (b standard-matrix)) - (let ((n-a (nrows a)) - (m-a (ncols a)) - (n-b (nrows b)) - (m-b (ncols b))) - (declare (type fixnum n-a m-a n-b m-b)) - - (unless (and (= n-a n-b) - (= m-a m-b)) - (error "Cannot subtract a ~d x ~d matrix from a ~d x ~d matrix" - n-b m-b - n-a m-a)))) - - -(defmethod m- ((a standard-matrix) (b standard-matrix)) - (axpy -1.0d0 b a)) - -(defmethod m- ((a number) (b standard-matrix)) - (error "cannot M- a matrix from a scalar")) - -(defmethod m- ((a standard-matrix) (b number)) - (m+ a (- b))) diff --git a/src/old/mplus.lisp b/src/old/mplus.lisp deleted file mode 100644 index bbe1229..0000000 --- a/src/old/mplus.lisp +++ /dev/null @@ -1,293 +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. -;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(in-package #:matlisp) - -(defgeneric m+ (a b) - (:documentation - " - Syntax - ====== - (M+ a b) - - Purpose - ======= - Create a new matrix which is the sum of A and B. - A or B (but not both) may be a scalar, in which - case the addition is element-wise. -") - (:method ((a number) (b number)) - (+ a b)) - (:method ((a standard-tensor) (b standard-tensor)) - (axpy 1 a b)) - (:method ( - -(definline m.+ (a b) - " - Syntax - ====== - (M.+ a b) - - Purpose - ======= - Same as M+ -" - (m+ a b)) - -(defgeneric m+! (a b) - (:documentation - " - Syntax - ====== - (M+! a b) - - Purpose - ======= - Desctructive version of M+: - - B <- A + B -") - (:method ((a number) (b number)) - (+ a b))) - -(definline m.+! (a b) - " - Syntax - ====== - (M.+! a b) - - Purpose - ======= - Same as M+! -" - (m+! a b)) - -(defmethod m+ :before ((a standard-matrix) (b standard-matrix)) - (let ((n-a (nrows a)) - (m-a (ncols a)) - (n-b (nrows b)) - (m-b (ncols b))) - (declare (type fixnum n-a m-a n-b m-b)) - - (unless (and (= n-a n-b) - (= m-a m-b)) - (error "Cannot add a ~d x ~d matrix and a ~d x ~d matrix" - n-a m-a - n-b m-b)))) - - -(defmethod m+ ((a standard-matrix) (b standard-matrix)) - (axpy 1.0d0 a b)) - -(let ((b-array (make-array 1 :element-type 'real-matrix-element-type))) - (defmethod m+ ((a real-matrix) (b #+(or cmu sbcl) double-float #-(or cmu sbcl) float)) - (let ((nxm (number-of-elements a)) - (result (copy a))) - (declare (type fixnum nxm)) - - (setf (aref b-array 0) b) - (daxpy nxm 1.0d0 b-array 0 (store result) 1) - result))) - -(defmethod m+ ((a real-matrix) (b cl:real)) - (m+ a (coerce b 'real-matrix-element-type))) - -(defmethod m+ ((a #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (b real-matrix)) - (m+ b a)) - -(defmethod m+ ((a cl:real) (b real-matrix)) - (m+ b (coerce a 'real-matrix-element-type))) - -(defmethod m+ ((a real-matrix) (b #+:cmu kernel::complex-double-float - #+:sbcl sb-kernel::complex-double-float - #-(or cmu sbcl) complex)) - (let* ((n (nrows a)) - (m (ncols a)) - #-(or cmu sbcl) (b (complex-coerce b)) - (result (make-complex-matrix-dim n m b))) - (declare (type fixnum n m)) - - (axpy! 1.0d0 a result))) - -#+(or :cmu :sbcl) -(defmethod m+ ((a real-matrix) (b complex)) - (m+ a (complex-coerce b))) - -(defmethod m+ ((a #+:cmu kernel::complex-double-float - #+:sbcl sb-kernel::complex-double-float - #-(or cmu sbcl) complex) (b real-matrix)) - (m+ b a)) - -#+(or :cmu :sbcl) -(defmethod m+ ((a complex) (b real-matrix)) - (m+ b (complex-coerce a))) - -;;; -(let ((b-array (make-array 1 :element-type 'real-matrix-element-type))) - (defmethod m+ ((a complex-matrix) (b #+(or cmu sbcl) double-float #-(or cmu sbcl) float)) - (let ((nxm (number-of-elements a)) - (result (copy a))) - (declare (type fixnum nxm)) - - (setf (aref b-array 0) b) - (daxpy nxm 1.0d0 b-array 0 (store result) 2) - result))) - -(defmethod m+ ((a complex-matrix) (b cl:real)) - (m+ a (coerce b 'complex-matrix-element-type))) - -(defmethod m+ ((a #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (b complex-matrix)) - (m+ b a)) - -(defmethod m+ ((a cl:real) (b complex-matrix)) - (m+ b (coerce a 'complex-matrix-element-type))) - -(defmethod m+ ((a complex-matrix) (b #+:cmu kernel::complex-double-float - #+:sbcl sb-kernel::complex-double-float - #-(or cmu sbcl) complex)) - (let* ((n (nrows a)) - (m (ncols a)) - #-(or cmu sbcl) (b (complex-coerce b)) - (result (make-complex-matrix-dim n m b))) - (declare (type fixnum n m)) - - (axpy! 1.0d0 a result))) - -#+(or :cmu :sbcl) -(defmethod m+ ((a complex-matrix) (b complex)) - (m+ a (complex-coerce b))) - -(defmethod m+ ((a #+:cmu kernel::complex-double-float - #+:sbcl sb-kernel::complex-double-float - #-(or cmu sbcl) complex) (b complex-matrix)) - (m+ b a)) - -#+(or :cmu :sbcl) -(defmethod m+ ((a complex) (b complex-matrix)) - (m+ b (complex-coerce a))) - - -(defmethod m+! :before ((a standard-matrix) (b standard-matrix)) - (let ((n-a (nrows a)) - (m-a (ncols a)) - (n-b (nrows b)) - (m-b (ncols b))) - (declare (type fixnum n-a m-a n-b m-b)) - (unless (and (= n-a n-b) - (= m-a m-b)) - (error "Cannot add a ~d x ~d matrix and a ~d x ~d matrix" - n-a m-a - n-b m-b)))) - -(defmethod m+! ((a standard-matrix) (b standard-matrix)) - (axpy! 1.0d0 a b)) - -(defmethod m+! ((a complex-matrix) (b real-matrix)) - (error "cannot M+! a COMPLEX-MATRIX A and a REAL-MATRIX B, -don't know how to coerce COMPLEX to REAL.")) - -;;; -(let ((b-array (make-array 1 :element-type 'real-matrix-element-type))) - (defmethod m+! ((a real-matrix) (b #+(or cmu sbcl) double-float #-(or cmu sbcl) float)) - (let ((nxm (number-of-elements a))) - (declare (type fixnum nxm)) - - (setf (aref b-array 0) b) - (daxpy nxm 1.0d0 b-array 0 (store a) 1) - a))) - -(defmethod m+! ((a real-matrix) (b cl:real)) - (m+! a (coerce b 'real-matrix-element-type))) - -(defmethod m+! ((a #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (b real-matrix)) - (m+! b a)) - -(defmethod m+! ((a cl:real) (b real-matrix)) - (m+! b (coerce a 'real-matrix-element-type))) - -(defmethod m+! ((a real-matrix) (b complex)) - (error "cannon M+! a REAL-MATRIX and a COMPLEX, -don't know how to coerce COMPLEX to REAL")) - -(defmethod m+! ((a complex) (b real-matrix)) - (error "cannon M+! a REAL-MATRIX and a COMPLEX, -don't know how to coerce COMPLEX to REAL")) - -(let ((b-array (make-array 1 :element-type 'real-matrix-element-type))) - (defmethod m+! ((a complex-matrix) (b #+(or cmu sbcl) double-float #-(or cmu sbcl) float)) - (let ((nxm (number-of-elements a))) - (declare (type fixnum nxm)) - - (setf (aref b-array 0) b) - (daxpy nxm 1.0d0 b-array 0 (store a) 2) - a))) - -(defmethod m+! ((a complex-matrix) (b cl:real)) - (m+! a (coerce b 'complex-matrix-element-type))) - -(defmethod m+! ((a #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (b complex-matrix)) - (m+! b a)) - -(defmethod m+! ((a cl:real) (b complex-matrix)) - (m+! b (coerce a 'complex-matrix-element-type))) - -#-:sbcl ;; sbcl doesn't like constant arrays -(defconstant *complex-unity-as-array* - (make-array 2 :element-type 'complex-matrix-element-type - :initial-contents '(1.0d0 0.0d0))) - -#+:sbcl -(defvar *complex-unity-as-array* - (make-array 2 :element-type 'complex-matrix-element-type - :initial-contents '(1.0d0 0.0d0))) - -(defmethod m+! ((a complex-matrix) (b #+:cmu kernel::complex-double-float - #+:sbcl sb-kernel::complex-double-float - #-(or cmu sbcl) complex)) - (let* ((nxm (number-of-elements a))) - (declare (type fixnum nxm)) - - #-(or cmu sbcl) (setq b (complex-coerce b)) - - (setf (aref *1x1-complex-array* 0) (realpart b)) - (setf (aref *1x1-complex-array* 1) (imagpart b)) - (zaxpy nxm #c(1d0 0) b 0 (store a) 1) - a)) - -#+(or :cmu :sbcl) -(defmethod m+! ((a complex-matrix) (b complex)) - (m+! a (complex-coerce b))) - -(defmethod m+! ((a #+:cmu kernel::complex-double-float - #+:sbcl sb-kernel::complex-double-float - #-(or cmu sbcl) complex) (b complex-matrix)) - (m+! b a)) - -#+(or :cmu :sbcl) -(defmethod m+! ((a complex) (b complex-matrix)) - (m+! b (complex-coerce a))) - diff --git a/src/sugar/mplusminus.lisp b/src/sugar/mplusminus.lisp new file mode 100644 index 0000000..1c65dff --- /dev/null +++ b/src/sugar/mplusminus.lisp @@ -0,0 +1,164 @@ +;;; -*- 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) + +(defgeneric m+ (a b) + (:documentation + " + Syntax + ====== + (M+ a b) + + Purpose + ======= + Create a new matrix which is the sum of A and B. + A or B (but not both) may be a scalar, in which + case the addition is element-wise. +") + (:method ((a number) (b number)) + (+ a b)) + (:method ((a standard-tensor) (b standard-tensor)) + (axpy 1 a b)) + (:method ((a number) (b standard-tensor)) + (axpy a nil b)) + (:method ((a standard-tensor) (b number)) + (axpy b nil a))) + +(definline m.+ (a b) + " + Syntax + ====== + (M.+ a b) + + Purpose + ======= + Same as M+ +" + (m+ a b)) +;;---------------------------------------------------------------;; + +(defgeneric m+! (a b) + (:documentation + " + Syntax + ====== + (M+! a b) + + Purpose + ======= + Desctructive version of M+: + + B <- A + B +") + (:method ((a number) (b number)) + (+ a b)) + (:method ((a standard-tensor) (b standard-tensor)) + (axpy! 1 a b)) + (:method ((a number) (b standard-tensor)) + (axpy! a nil b))) + +(definline m.+! (a b) + " + Syntax + ====== + (M.+! a b) + + Purpose + ======= + Same as M+! +" + (m+! a b)) +;;---------------------------------------------------------------;; + +(defgeneric m- (a b) + (:documentation + " + Syntax + ====== + (M- a b) + + Purpose + ======= + Create a new matrix which is the difference of A and B. + B may be a scalar, in which case the subtraction + is elementwise. +") + (:method ((a number) (b number)) + (- a b)) + (:method ((a standard-tensor) (b standard-tensor)) + (axpy -1 b a)) + (:method ((a number) (b standard-tensor)) + (scal! -1 (axpy (- a) nil b))) + (:method ((a standard-tensor) (b number)) + (axpy (- b) nil a))) + +(definline m.- (a b) + " + Syntax + ====== + (M.- a b) + + Purpose + ======= + Same as M- +" + (m- a b)) +;;---------------------------------------------------------------;; + +(defgeneric m-! (a b) + (:documentation + " + Syntax + ====== + (M-! a b) + + Purpose + ======= + Desctructive version of M-: + + B <- A - B +") + (:method ((a number) (b number)) + (- a b)) + (:method ((a standard-tensor) (b standard-tensor)) + (scal! -1 (axpy! -1 a b))) + (:method ((a number) (b standard-tensor)) + (scal! -1 (axpy! (- a) nil b)))) + +(definline m.-! (a b) + " + Syntax + ====== + (M.-! a b) + + Purpose + ======= + Same as M-! +" + (m+! a b)) +;;---------------------------------------------------------------;; ----------------------------------------------------------------------- Summary of changes: matlisp.asd | 12 ++- src/base/permutation.lisp | 1 + src/base/standard-tensor.lisp | 18 +++ src/foreign-core/lapack.lisp | 2 +- src/{old => lapack}/gesv.lisp | 116 ++++------------ src/{old => lapack}/getrf.lisp | 0 src/{old => lapack}/getrs.lisp | 0 src/level-1/axpy.lisp | 13 ++ src/old/mminus.lisp | 118 ---------------- src/old/mplus.lisp | 293 ---------------------------------------- src/sugar/mplusminus.lisp | 164 ++++++++++++++++++++++ 11 files changed, 236 insertions(+), 501 deletions(-) rename src/{old => lapack}/gesv.lisp (70%) rename src/{old => lapack}/getrf.lisp (100%) rename src/{old => lapack}/getrs.lisp (100%) delete mode 100644 src/old/mminus.lisp delete mode 100644 src/old/mplus.lisp create mode 100644 src/sugar/mplusminus.lisp hooks/post-receive -- matlisp |