Author: mantoniotti
Date: Wed Aug 16 17:07:41 2006
New Revision: 1
Added:
trunk/
trunk/common-math/
trunk/common-math/common-math-pkg.lisp
trunk/common-math/common-math-user-pkg.lisp
trunk/common-math/common-math.lisp
trunk/common-math/common-math.system
trunk/common-math/numerics/
trunk/common-math/numerics/linear-algebra/
trunk/common-math/numerics/linear-algebra/common-math-extra.lisp
trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp
trunk/common-math/numerics/linear-algebra/linear-algebra.system
trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp
trunk/common-math/numerics/linear-algebra/matrix-pkg.lisp
trunk/common-math/numerics/linear-algebra/matrix.lisp
trunk/common-math/numerics/linear-algebra/vector-matrix-conditions.lisp
trunk/common-math/numerics/linear-algebra/vector-matrix.lisp
trunk/common-math/numerics/linear-algebra/vector.lisp
trunk/common-math/numerics/numerics.system
Log:
Initial import.
Added: trunk/common-math/common-math-pkg.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/common-math-pkg.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,47 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; common-math-pkg.lisp --
+
+(defpackage "COMMON-MATH" (:use "COMMON-LISP")
+ (:documentation "The COMMON MATH Package.
+This package contains `top level' definitions of a `generic math'
+package which can be used as a drop in replacement for the standard
+CL mathematics routines. Note however, that this will of course
+come at a price.
+
+The COMMON-MATH package addresses the problem of having a generic and
+as complete as possible mathematics interface layer over a number of
+different algorithms and representations. In this sense it aspires
+to bring to CL the flexibility of generic math operations, while
+preserving the low level hooks that could be used for (relative) efficiency.")
+
+ (:nicknames "CL.MATH" "CL.MATHEMATICS")
+
+ (:shadow "=" "+" "-" "*" "/" ">" "<" ">=" "<=")
+ (:shadow "ZEROP" "EXPT" "GCD")
+
+ (:export "=" "+" "-" "*" "/" ">" "<" ">=" "<=")
+ (:export "ZEROP" "EXPT" "GCD")
+
+ (:export "+POSITIVE-INFINITY+" "+NEGATIVE-INFINITY+")
+
+ (:export "INCREASING" "DECREASING")
+
+ (:export
+ "+%1" "+%2"
+ "*%1" "*%2"
+ "-%1" "-%2"
+ "/%1" "/%2"
+ "=%1" "=%2"
+ "<%2"
+ ">%2"
+ "<=%2"
+ ">=%2"
+ )
+
+ ;; Conditions
+ (:export "UNDEFINED-OPERATION")
+
+ )
+
+;;;; end of file -- common-math-pkg.lisp --
Added: trunk/common-math/common-math-user-pkg.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/common-math-user-pkg.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,12 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; common-math-user-pkg.lisp --
+;;;; The COMMON-MATH-USER package.
+
+(defpackage "COMMON-MATH-USER" (:use "COMMON-MATH" "CL" #| "LINALG" |# )
+ (:shadowing-import-from "COMMON-MATH" "+" "-" "*" "/")
+ (:shadowing-import-from "COMMON-MATH" "<" ">" "<=" ">=" "=")
+ (:shadowing-import-from "COMMON-MATH" "ZEROP" "EXPT" "GCD")
+ )
+
+;;;; end of file -- common-math-user-pkg.lisp --
\ No newline at end of file
Added: trunk/common-math/common-math.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/common-math.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,603 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; common-math.lisp --
+
+(in-package "CL.MATHEMATICS")
+
+;;; undefined-operation --
+
+(define-condition undefined-operation (error)
+ ((operator :reader undefined-operation-operator
+ :initarg :operator)
+ (arguments :reader undefined-operation-arguments
+ :initarg :arguments)
+ )
+ (:documentation "The Undefined Operation Condition.")
+ (:report (lambda (uoc stream)
+ (format stream "Undefined operation ~S called with ~S."
+ (undefined-operation-operator uoc)
+ (undefined-operation-arguments uoc))))
+ (:default-initargs :arguments () :operator nil))
+
+
+;;;---------------------------------------------------------------------------
+;;; Constants and parameters.
+
+(defconstant +negative-infinity+ '+negative-infinity+
+ "The symbolic +NEGATIVE-INFINITY+ constant.
+A representation of negative infinity.")
+(defconstant +positive-infinity+ '+positive-infinity+
+ "The symbolic +POSITIVE-INFINITY+ constant.
+A representation of positive infinity.")
+
+(defparameter *ignore-comparison-errors-p* nil)
+
+
+;;; (defconstant +negative-infinity+ most-negative-fixnum)
+;;; (defconstant +positive-infinity+ most-positive-fixnum)
+
+
+;;;---------------------------------------------------------------------------
+;;; Generic functions.
+;;; Note that the generic functions corresponding to traditional math
+;;; operators are named using a special notation to indicate their
+;;; arity. Such 'operator' generic function are named as
+;;;
+;;; <op>%<arity>
+;;;
+;;; The #\% character in the name is used to suggest "low level" as
+;;; usual. The Prolog-like #\/ could have been used instead but the
+;;; #\% has more tradition in CL circles. The #\@ would have made
+;;; more sense, but some CL implementations already use it for other
+;;; purposes and setting up a separate readtable seems a little bit
+;;; too much for the time being.
+
+(defgeneric <%2 (x y)
+ (:documentation "The <%2 generic function.
+The binary LESS generic function which is specialized for various
+combinations of argument types.
+
+At a minimum there is a guarantee that the method on two NUMBER
+arguments dispatching on CL:< is defined.")
+ (:method ((x number) (y number))
+ (cl:< x y))
+
+ (:method ((x number) (y (eql +negative-infinity+)))
+ nil)
+
+ (:method ((y (eql +negative-infinity+)) (x number))
+ T)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)))
+ T)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)))
+ nil)
+
+ (:method ((x number) (y (eql +positive-infinity+)))
+ T)
+
+ (:method ((y (eql +positive-infinity+)) (x number))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+ )
+
+
+(defgeneric >%2 (x y)
+ (:method ((x number) (y number))
+ (cl:> x y))
+
+ (:method ((x number) (y (eql +negative-infinity+)))
+ T)
+
+ (:method ((y (eql +negative-infinity+)) (x number))
+ nil)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)))
+ T)
+
+ (:method ((x number) (y (eql +positive-infinity+)))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x number))
+ T)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+ )
+
+
+(defmethod <%2 ((x symbol) (y symbol))
+ (string< x y))
+
+
+(defmethod >%2 ((x symbol) (y symbol))
+ (string> y x))
+
+
+(defgeneric <=%2 (x y)
+ (:method ((x number) (y number))
+ (cl:<= x y))
+
+ (:method ((x number) (y (eql +negative-infinity+)))
+ nil)
+
+ (:method ((y (eql +negative-infinity+)) (x number))
+ T)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)))
+ T)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)))
+ nil)
+
+ (:method ((x number) (y (eql +positive-infinity+)))
+ T)
+
+ (:method ((y (eql +positive-infinity+)) (x number))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+ )
+
+
+(defgeneric >=%2 (x y)
+ (:method ((x number) (y number))
+ (cl:>= x y))
+
+ (:method ((x number) (y (eql +negative-infinity+)))
+ T)
+
+ (:method ((y (eql +negative-infinity+)) (x number))
+ nil)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)))
+ T)
+
+ (:method ((x number) (y (eql +positive-infinity+)))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x number))
+ T)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator 'lessp%2)))
+ )
+
+
+(defgeneric +%2 (x y &optional r)
+ (:method ((x number) (y number) &optional r)
+ (declare (ignore r))
+ (cl:+ x y))
+
+ (:method ((x number) (y (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ (:method ((y (eql +negative-infinity+)) (x number) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))
+ &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))
+ &optional r)
+ (declare (ignore r))
+ (error 'undefined-operation :operator '+%2))
+
+ (:method ((x number) (y (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((y (eql +positive-infinity+)) (x number) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))
+ &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))
+ &optional r)
+ (declare (ignore r))
+ (error 'undefined-operation :operator '+%2))
+ )
+
+
+(defgeneric *%2 (x y &optional r)
+ (:method ((x number) (y number) &optional r)
+ (declare (ignore r))
+ (cl:* x y))
+
+ (:method ((x number) (y (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ (cond ((zerop x) (error 'undefined-operation :operator '*%2))
+ ((plusp x) +negative-infinity+)
+ (t +positive-infinity+)))
+
+ (:method ((y (eql +negative-infinity+)) (x number) &optional r)
+ (declare (ignore r))
+ (cond ((zerop x) (error 'undefined-operation :operator '*%2))
+ ((plusp x) +negative-infinity+)
+ (t +positive-infinity+)))
+
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))
+ &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))
+ &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ (:method ((x number) (y (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ (cond ((zerop x) (error 'undefined-operation :operator '*%2))
+ ((plusp x) +positive-infinity+)
+ (t +negative-infinity+)))
+
+ (:method ((y (eql +positive-infinity+)) (x number) &optional r)
+ (declare (ignore r))
+ (cond ((zerop x) (error 'undefined-operation :operator '*%2))
+ ((plusp x) +positive-infinity+)
+ (t +negative-infinity+)))
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))
+ &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))
+ &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+ )
+
+
+(defgeneric -%2 (x y &optional r)
+ (:method ((x t) (y t) &optional r)
+ (declare (ignore r))
+ (+%2 x (-%1 y))))
+
+
+(defgeneric /%2 (x y &optional r)
+ (:method ((x number) (y number) &optional r)
+ (declare (ignore r))
+ (cl:/ x y))
+
+ (:method ((x t) (y t) &optional r)
+ (declare (ignore r))
+ (*%2 x (/%1 y))))
+
+
+(defgeneric +%1 (x &optional r)
+ (:method ((x number) &optional r)
+ (declare (ignore r))
+ (cl:+ x))
+
+ (:method ((x (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ )
+
+(defgeneric *%1 (x &optional r)
+ (:method ((x number) &optional r)
+ (declare (ignore r))
+ (cl:* x))
+
+ (:method ((x (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ )
+
+(defgeneric -%1 (x &optional r)
+ (:method ((x number) &optional r)
+ (declare (ignore r))
+ (cl:- x))
+
+ (:method ((x (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ (:method ((x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ )
+
+(defgeneric /%1 (x &optional r)
+ (:method ((x number) &optional r)
+ (declare (ignore r))
+ (cl:/ x))
+
+ (:method ((x (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ 0)
+
+ (:method ((x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ 0)
+ )
+
+
+(defgeneric =%2 (x y)
+ (:method ((x number) (y number)) (cl:= x y))
+
+ (:method ((x number) (y (eql +negative-infinity+)))
+ nil)
+
+ (:method ((y (eql +negative-infinity+)) (x number))
+ nil)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator '=%2)))
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)))
+ nil)
+
+ (:method ((x number) (y (eql +positive-infinity+)))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x number))
+ nil)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)))
+ (unless *ignore-comparison-errors-p*
+ (error 'undefined-operation :operator '=%2)))
+ )
+
+
+(defgeneric =%1 (x)
+ (:method ((x t)) T)
+ )
+
+
+(defgeneric gcd%2 (x y)
+ (:method ((x integer) (y integer))
+ (cl:gcd x y)))
+
+
+;;;---------------------------------------------------------------------------
+;;; Redefined Common Lisp operators.
+
+(defun + (&rest args)
+ (if (null args)
+ (cl:+)
+ (let ((n-args (list-length args)))
+ (cond ((cl:= n-args 1)
+ (+%1 (first args)))
+ ((cl:= n-args 2)
+ (+%2 (first args) (second args)))
+ (t (+%2 (first args) (apply #'+ (rest args))))
+ ))))
+
+
+(defun - (arg1 &rest args)
+ (if (null args)
+ (-%1 arg1)
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (-%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of summation.
+ (loop for arg in args
+ for minuend =(-%2 arg1 arg) then (-%2 minuend arg)
+ finally (return minuend))))
+ ))
+
+
+(defun * (&rest args)
+ (if (null args)
+ (cl:*)
+ (let ((n-args (list-length args)))
+ (cond ((cl:= n-args 1)
+ (*%1 (first args)))
+ ((cl:= n-args 2)
+ (*%2 (first args) (second args)))
+ (t (*%2 (first args) (apply #'* (rest args))))
+ ))))
+
+
+(defun / (arg1 &rest args)
+ (if (null args)
+ (/%1 arg1)
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (/%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of divisions.
+ (loop for arg in args
+ for dividend = (/%2 arg1 arg) then (/%2 dividend arg)
+ finally (return dividend))))
+ ))
+
+
+(defun = (arg1 &rest args)
+ (if (null args)
+ (=%1 arg1)
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (=%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of the equality test.
+ (loop for arg in args
+ always (=%2 arg1 arg))))
+ ))
+
+
+(defun < (arg1 &rest args)
+ (if (null args)
+ ;; (<%1 arg1)
+ t
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (<%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of the equality test.
+ (loop for a1 = arg1 then a2
+ for a2 in args
+ always (<%2 a1 a2))))
+ ))
+
+
+(defun > (arg1 &rest args)
+ (if (null args)
+ t
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (>%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of the equality test.
+ (loop for a1 = arg1 then a2
+ for a2 in args
+ always (>%2 a1 a2))))
+ ))
+
+
+(defun <= (arg1 &rest args)
+ (if (null args)
+ t
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (<=%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of the equality test.
+ (loop for a1 = arg1 then a2
+ for a2 in args
+ always (<=%2 a1 a2))))
+ ))
+
+
+(defun >= (arg1 &rest args)
+ (if (null args)
+ t
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (>=%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of the equality test.
+ (loop for a1 = arg1 then a2
+ for a2 in args
+ always (>=%2 a1 a2))))
+ ))
+
+
+(defun gcd (&rest args)
+ (if (null args)
+ (cl:gcd)
+ (let ((n-args (list-length args)))
+ (cond ((cl:= n-args 1)
+ (gcd%2 (first args) (first args)))
+ ((cl:= n-args 1)
+ (gcd%2 (first args) (second args)))
+ (t
+ (apply #'gcd
+ (gcd%2 (first args) (second args))
+ (cddr args)))
+ ))))
+
+
+(defgeneric expt (base power)
+ (:method ((base number) (power number))
+ (cl:expt base power)))
+
+
+(defgeneric zerop (n)
+ (:method ((n number))
+ (cl:zerop n))
+
+ (:method ((x (eql +positive-infinity+)))
+ nil)
+
+ (:method ((x (eql +negative-infinity+)))
+ nil)
+ )
+
+
+
+;;;---------------------------------------------------------------------------
+;;; Other functions.
+
+
+(defun increasing (x)
+ (declare (type sequence x))
+ (loop for i from 0 below (1- (length x))
+ for a1 = (elt x 0) then a2
+ for a2 = (elt x (1+ i))
+ always (<=%2 a1 a2)))
+
+
+(defun decreasing (x)
+ (declare (type sequence x))
+ (loop for i from 0 below (1- (length x))
+ for a1 = (elt x 0) then a2
+ for a2 = (elt x (1+ i))
+ always (>=%2 a1 a2)))
+
+
+
+(defun strictly-increasing (x)
+ (declare (type sequence x))
+ (loop for i from 0 below (1- (length x))
+ for a1 = (elt x 0) then a2
+ for a2 = (elt x (1+ i))
+ always (<%2 a1 a2)))
+
+
+(defun strictly-decreasing (x)
+ (declare (type sequence x))
+ (loop for i from 0 below (1- (length x))
+ for a1 = (elt x 0) then a2
+ for a2 = (elt x (1+ i))
+ always (>%2 a1 a2)))
+
+;;; end of file -- common-math.lisp --
Added: trunk/common-math/common-math.system
==============================================================================
--- (empty file)
+++ trunk/common-math/common-math.system Wed Aug 16 17:07:41 2006
@@ -0,0 +1,24 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; common-math.system --
+
+(eval-when (:load-toplevel :execute)
+ (mk:add-registry-location (make-pathname :name nil
+ :type nil
+ :defaults *load-truename*))
+ )
+
+
+(mk:defsystem "common-math"
+ :components ((:file "common-math-pkg")
+ (:file "common-math"
+ :depends-on ("common-math-pkg"))
+ (:system "numerics")
+ ;; (:system "polynomials")
+ ;; (:system "formulas")
+ (:file "common-math-user-pkg"
+ :depends-on ("numerics"))
+ )
+ )
+
+;;;; end of file -- common-math.system --
Added: trunk/common-math/numerics/linear-algebra/common-math-extra.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/common-math-extra.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,98 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; common-math-extra.lisp --
+
+(in-package "CL.MATH.LINEAR-ALGEBRA")
+
+;;;---------------------------------------------------------------------------
+;;; More functions. Essentially the top level element-wise array operations.
+
+(defgeneric .*%2 (x y &optional r)
+ (:method ((x number) (y number) &optional r)
+ (declare (ignore r))
+ (*%2 x y))
+
+ (:method ((x (eql +negative-infinity+)) (y number) &optional r)
+ (declare (ignore r))
+ (*%2 x y))
+
+ (:method ((x (eql +positive-infinity+)) (y number) &optional r)
+ (declare (ignore r))
+ (*%2 x y))
+
+
+ (:method ((y number) (x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ (*%2 x y))
+
+
+ (:method ((y number) (x (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ (*%2 x y))
+
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+
+ (:method ((x (eql +negative-infinity+)) (y (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+ )
+
+
+(defgeneric .*%1 (x &optional r)
+ (:method ((x number) &optional r)
+ (declare (ignore r))
+ (*%1 x))
+
+ (:method ((x (eql +positive-infinity+)) &optional r)
+ (declare (ignore r))
+ +positive-infinity+)
+
+ (:method ((x (eql +negative-infinity+)) &optional r)
+ (declare (ignore r))
+ +negative-infinity+)
+ )
+
+
+(defgeneric ./%2 (x y &optional r))
+
+(defgeneric ./%1 (x &optional r))
+
+
+;;;---------------------------------------------------------------------------
+;;; The generic top-level functions.
+
+(defun .* (&rest args)
+ (if (null args)
+ (cl:*)
+ (let ((n-args (list-length args)))
+ (cond ((cl:= n-args 1) (.*%1 (first args)))
+ ((cl:= n-args 2) (.*%2 (first args) (second args)))
+ (t (.*%2 (first args) (apply #'.* (rest args))))
+ ))))
+
+
+(defun ./ (arg1 &rest args)
+ (if (null args)
+ (./%1 arg1)
+ (let ((n-args (list-length args)))
+ (if (cl:= n-args 1)
+ (./%2 arg1 (first args))
+ ;; The loop is done in order not to assume anything about
+ ;; the nature of divisions.
+ (loop for arg in args
+ for dividend = (./%2 arg1 arg) then (./%2 dividend arg)
+ finally (return dividend))))
+ ))
+
+;;;; end of file -- common-math-extra.lisp --
Added: trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,19 @@
+;;; -*- Mode: Lisp -*-
+
+(defpackage "CL.MATH.LINEAR-ALGEBRA" (:use "CL.MATH" "COMMON-LISP")
+ (:nicknames "CL-MATH-LINALG" "LINALG")
+
+ (:shadow cl:zerop cl:expt cl:gcd)
+
+ (:shadowing-import-from "CL.MATH" "+" "-" "*" "/" "=" ">" "<" ">=" "<=")
+
+ (:export "=" "+" "-" "*" "/" ">" "<" ">=" "<=")
+ (:export "ZEROP" "EXPT" "GCD")
+
+ (:export ".*" "./" "EXPT*")
+
+ (:export ".*%2" ".*%1" "./%2" "./%1")
+
+ )
+
+;;; end of file -- linear-algebra-pkg.lisp --
Added: trunk/common-math/numerics/linear-algebra/linear-algebra.system
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/linear-algebra.system Wed Aug 16 17:07:41 2006
@@ -0,0 +1,13 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; linear-algebra.system --
+
+(mk:defsystem "linear-algebra"
+ :components ("linear-algebra-pkg"
+ "vector-matrix-conditions"
+ "common-math-extra"
+ "vector"
+ "matrix"
+ "lu-decomposition"))
+
+;;;; end of file -- linear-algebra.system --
Added: trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,326 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; lu-decomposition.lisp --
+
+(in-package "CL.MATH.LINEAR-ALGEBRA")
+
+(define-condition singular-matrix-in-operation (arithmetic-error)
+ ()
+ (:report (lambda (cnd stream)
+ (format stream "Singular matrix passed to operation ~S."
+ (arithmetic-error-operation cnd)))))
+
+
+(defparameter *epsilon* 1.0e-20) ; Change this to CL constant.
+
+
+(defgeneric square-matrix-p (m)
+ (:method ((m array))
+ (and (matrix-array-p m)
+ (cl:= (array-dimension m 0) (array-dimension m 1))))
+ (:method ((m matrix))
+ (cl:= (matrix-row-number m) (matrix-column-number m))))
+
+
+
+;;; Different base type.
+;;; Probably a macrolet is in order.
+
+(defun lu-decompose (m &optional
+ (lu (make-array (array-dimensions m)
+ :element-type (array-element-type m)
+ :initial-element (coerce 0 (array-element-type m))))
+ (permutations (make-array (array-dimension m 0)
+ :initial-element 0
+ :element-type 'fixnum)))
+ (declare (type (array number (cl:* cl:*)) m lu)
+ (type (vector fixnum) permutations))
+
+ (assert (square-matrix-p m))
+ (assert (square-matrix-p lu))
+ (assert (cl:= (length permutations) (array-dimension m 0)))
+
+ (let* ((n (array-dimension m 0))
+ (row-scaling-vector (make-array n
+ :initial-element (coerce 0 (array-element-type m))
+ :element-type (array-element-type m)))
+ (parity 1)
+ (imax 0)
+ )
+ (declare (dynamic-extent row-scaling-vector)
+ (type fixnum imax n)
+ (type (member -1 1) parity)
+ (type (simple-array number (cl:*)) row-scaling-vector))
+
+ ;; Copy M into LU and then just operate the in-place algorithm on LU.
+ ;; If M is EQ to LU, that means that the argumen was supplied and
+ ;; that the user actually wants a destructive operation. Copying
+ ;; can be avoided in this case.
+ (unless (eq m lu)
+ (dotimes (i n)
+ (dotimes (j n)
+ (setf (aref lu i j) (aref m i j)))))
+
+
+ ;; Loop over rows to get the implicit scaling information.
+ ;; Meanwhile check that the matrix is non-singular.
+ (dotimes (i n)
+ (loop for k of-type fixnum from 0 below n
+ maximize (abs (aref lu i k)) into big ;; of-type (number) ; not CL.
+ finally
+ (when (cl:zerop big)
+ (error 'singular-matrix-in-operation
+ :operation 'lu-decompose
+ :operands (list m)))
+ (setf (aref row-scaling-vector i) (/ big))))
+
+ ;; Loop over the columns of Crout's method.
+ (dotimes (j n)
+ (block operate-on-column
+ ;; (format t "Column ~D IMAX = ~D~%" j imax)
+ (dotimes (i j)
+ (let ((sum (aref lu i j)))
+ (declare (type number sum))
+ (dotimes (k i)
+ (decf sum (* (aref lu i k) (aref lu k j))))
+ (setf (aref lu i j) sum)))
+
+ ;; Initialize the search for the largest pivot.
+ (loop with big of-type number = 0.0
+ for i from j below n
+ do (let ((sum (aref lu i j)))
+ (declare (type number sum))
+ (dotimes (k j)
+ (decf sum (* (aref lu i k) (aref lu k j))))
+ (setf (aref lu i j) sum)
+ (let ((temp (* (aref row-scaling-vector i) (abs sum))))
+ (declare (type number temp))
+ (when (>= temp big)
+ (setf big temp
+ imax i)))
+ ))
+
+ ;; (format t "Now IMAX = ~D~%" imax)
+ (when (/= imax j)
+ ;; We need to interchange rows.
+ (dotimes (k n)
+ (rotatef (aref lu imax k) (aref lu j k)))
+ (setf parity (- parity))
+ (setf (aref row-scaling-vector imax)
+ (aref row-scaling-vector j)))
+ (setf (aref permutations j) imax)
+ ;; (format t "Now J = ~D IMAX = ~D~2%" j imax)
+ (when (cl:zerop (aref lu j j))
+ (setf (aref lu j j) *epsilon*))
+
+ (when (/= j (1- n))
+ (loop with pivot of-type number = (/ (aref lu j j))
+ for i from (1+ j) below n
+ do (setf (aref lu i j) (* pivot (aref lu i j)))))
+
+ ))
+ (values lu permutations parity)
+ ))
+
+
+(defun lu-decompose/double-float (m &optional
+ (lu (make-array (array-dimensions m)
+ :element-type 'double-float
+ :initial-element 0.0d0))
+ (permutations (make-array (array-dimension m 0)
+ :initial-element 0
+ :element-type 'fixnum)))
+ (declare (type (array double-float (cl:* cl:*)) m lu)
+ (type (vector fixnum) permutations))
+
+ (assert (square-matrix-p m))
+ (assert (square-matrix-p lu))
+ (assert (cl:= (length permutations) (array-dimension m 0)))
+
+ (let* ((n (array-dimension m 0))
+ (row-scaling-vector (make-array n
+ :initial-element (coerce 0 (array-element-type m))
+ :element-type (array-element-type m)))
+ (parity 1.0d0)
+ (imax 0)
+ )
+ (declare (dynamic-extent row-scaling-vector)
+ (type fixnum n imax)
+ (type (double-float -1.01d0 1.01d0) parity)
+ (type (simple-array double-float (cl:*)) row-scaling-vector))
+
+ ;; Copy M into LU and then just operate the in-place algorithm on LU.
+ ;; If M is EQ to LU, that means that the argumen was supplied and
+ ;; that the user actually wants a destructive operation. Copying
+ ;; can be avoided in this case.
+ (unless (eq m lu)
+ (dotimes (i n)
+ (dotimes (j n)
+ (setf (aref lu i j) (aref m i j)))))
+
+
+ ;; Loop over rows to get the implicit scaling information.
+ ;; Meanwhile check that the matrix is non-singular.
+ (dotimes (i n)
+ (loop for k of-type fixnum from 0 below n
+ maximize (abs (aref lu i k)) into big of-type double-float
+ finally
+ (when (cl:zerop big)
+ (error 'singular-matrix-in-operation
+ :operation 'lu-decompose
+ :operands (list m)))
+ (setf (aref row-scaling-vector i) (/ big))))
+
+ ;; Loop over the columns of Crout's method.
+ (dotimes (j n)
+ (block operate-on-column
+ ;; (format t "Column ~D IMAX = ~D~%" j imax)
+ (dotimes (i j)
+ (let ((sum (aref lu i j)))
+ (declare (type double-float sum))
+ (dotimes (k i)
+ (decf sum (* (aref lu i k) (aref lu k j))))
+ (setf (aref lu i j) sum)))
+
+ ;; Initialize the search for the largest pivot.
+ (loop with big of-type double-float = 0.0d0
+ for i of-type fixnum from j below n
+ do (let ((sum (aref lu i j)))
+ (declare (type double-float sum))
+ (dotimes (k j)
+ (decf sum (* (aref lu i k) (aref lu k j))))
+ (setf (aref lu i j) sum)
+ (let ((temp (* (aref row-scaling-vector i) (abs sum))))
+ (declare (type double-float temp))
+ (when (>= temp big)
+ (setf big temp
+ imax i)))
+ ))
+
+ ;; (format t "Now IMAX = ~D~%" imax)
+ (when (/= imax j)
+ ;; We need to interchange rows.
+ (dotimes (k n)
+ (rotatef (aref lu imax k) (aref lu j k)))
+ (setf parity (- parity))
+ (setf (aref row-scaling-vector imax)
+ (aref row-scaling-vector j)))
+ (setf (aref permutations j) imax)
+ ;; (format t "Now J = ~D IMAX = ~D~2%" j imax)
+ (when (cl:zerop (aref lu j j))
+ (setf (aref lu j j) *epsilon*))
+
+ (when (/= j (1- n))
+ (loop with pivot of-type double-float = (/ (aref lu j j))
+ for i from (1+ j) below n
+ do (setf (aref lu i j) (* pivot (aref lu i j)))))
+
+ ))
+ (values lu permutations parity)
+ ))
+
+
+(defun split-lu-matrix (lu)
+ (let ((l (make-array (array-dimensions lu) :initial-element 0))
+ (u (make-array (array-dimensions lu) :initial-element 0))
+ (n (array-dimension lu 0))
+ )
+ (dotimes (i n)
+ (dotimes (j i)
+ (setf (aref l i j) (aref lu i j)))
+ (setf (aref l i i) 1))
+
+ (dotimes (i n)
+ (loop for j from i below n
+ do (setf (aref u i j) (aref lu i j))))
+
+ (values l u)))
+
+
+;;;---------------------------------------------------------------------------
+;;; lu-back-substitution
+
+(defun lu-back-substitution (lu permutations b)
+ (let ((n (array-dimension lu 0)))
+ (loop with ii of-type fixnum = 0
+ for i of-type fixnum from 0 below n
+ for ip of-type number = (aref permutations i)
+ for sum of-type number = (aref b ip)
+ do (setf (aref b ip) (aref b i))
+ do (cond ((not (cl:zerop ii))
+ (loop for j of-type fixnum from ii below (1- i)
+ do (decf sum (* (aref lu i j) (aref b i)))))
+ ((not (cl:zerop sum))
+ (setf (aref b i) sum)))
+ )
+ (loop for i of-type fixnum from (1- n) downto 0
+ for sum of-type number = (aref b i)
+ do (loop for j of-type fixnum from (1+ i) below n
+ do (decf sum (* (aref lu i j) (aref b i))))
+ do (setf (aref b i) (/ sum (aref lu i i))))
+ )
+ b)
+
+
+;;;---------------------------------------------------------------------------
+;;; solve
+
+(defmethod solve ((a array) (b vector) &optional (result (copy-seq b)))
+ (multiple-value-bind (lu permutations)
+ (lu-decompose a)
+ (lu-back-substitution lu permutations result)))
+
+(defmethod solve ((a matrix) (b vector) &optional (result (copy-seq b)))
+ (solve (matrix-data a) b result))
+
+
+(defmethod solve/double-float ((a array) (b vector) &optional (result (copy-seq b)))
+ (multiple-value-bind (lu permutations)
+ (lu-decompose/double-float a)
+ (lu-back-substitution lu permutations result)))
+
+(defmethod solve/double-float ((a matrix) (b vector) &optional (result (copy-seq b)))
+ (solve/double-float (matrix-data a) b result))
+
+
+;;;---------------------------------------------------------------------------
+;;; det
+
+(defmethod det ((a array))
+ (assert (square-matrix-p a))
+ (multiple-value-bind (lu permutations parity)
+ (lu-decompose a)
+ (declare (ignore permutations))
+ (dotimes (i (array-dimension lu 0) parity)
+ (setf parity (* parity (aref lu i i))))))
+
+(defmethod det ((a matrix))
+ (det (matrix-data a)))
+
+
+(defmethod det/double-float ((a array))
+ (assert (square-matrix-p a))
+ (multiple-value-bind (lu permutations parity)
+ (lu-decompose a)
+ (declare (type (array double-float (cl:* cl:*)) lu)
+ (ignore permutations)
+ (type double-float parity))
+ (dotimes (i (array-dimension lu 0) parity)
+ (setf parity (* parity (aref lu i i))))))
+
+(defmethod det/double-float ((a matrix))
+ (det/double-float (matrix-data a)))
+
+
+;;;===========================================================================
+;;; Test
+
+(defvar A #2A((1 2 3)
+ (2 -1 1)
+ (3 4 -1)))
+
+(defvar LU-A #2A((1 2 3)
+ (2 -5 -5)
+ (3 0.4 -8)))
+
+;;; end of file -- lu-decomposition.lisp --
Added: trunk/common-math/numerics/linear-algebra/matrix-pkg.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/matrix-pkg.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,19 @@
+;;; -*- Mode: Lisp -*-
+
+(defpackage "CL.MATH.POLYNOMIALS" (:use "CL.MATH" "COMMON-LISP")
+ (:nicknames "POLY")
+
+ (:shadow cl:zerop cl:expt cl:gcd)
+
+ (:shadowing-import-from "CL.MATH" "+" "-" "*" "/" "=" ">" "<" ">=" "<=")
+
+ (:export "=" "+" "-" "*" "/" ">" "<" ">=" "<=")
+ (:export "ZEROP" "EXPT" "GCD")
+
+ (:export "MAKE-POLYNOMIAL"
+ "POLYNOMIAL"
+ "UNIVARIATE-POLYNOMIAL"
+ )
+ )
+
+;;; end of file -- polymonials-pkg.lisp --
Added: trunk/common-math/numerics/linear-algebra/matrix.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/matrix.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,950 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; matrix.lisp --
+
+(in-package "CL.MATH.LINEAR-ALGEBRA")
+
+(deftype matrix-array (&optional (type 'number) (m 'cl:*) (n 'cl:*))
+ `(array ,type (,m ,n)))
+
+
+(defparameter *default-matrix-element-type* 'double-float)
+
+(defclass matrix ()
+ ((data :reader matrix-data :initarg :data))
+ (:default-initargs :data (make-array '(0 0) :element-type 'number)))
+
+(defgeneric matrix-p (x
+ #|&optional
+ (element-type *default-matrix-element-type*)
+ |#
+ )
+ (:method ((x matrix)) t)
+ (:method ((x t)) nil)
+ (:method ((x vector)) t) ; Note that this is technically correct, but
+ ; it may cause problems.
+ (:method ((x array)) (cl:<= (array-rank x) 2)) ; Same here.
+ )
+
+
+(defvar *check-element-type* nil)
+
+
+;;;---------------------------------------------------------------------------
+;;; Protocol.
+
+(defgeneric matrix-array-p (x &key check-element-type)
+ (:method ((x array) &key (check-element-type *check-element-type*))
+ (and (cl:<= (array-rank x) 2)
+ (or (not check-element-type)
+ (nth-value 0 (subtypep (array-element-type x) 'number)))))
+
+ (:method ((x t) &key check-element-type)
+ (declare (ignore check-element-type))
+ nil))
+
+
+
+(defgeneric matrix-element-type (m))
+
+(defgeneric matrix-row-number (m)
+ (:method ((x matrix)) (array-dimension (matrix-data x) 0))
+ (:method ((x vector)) 1)
+ (:method ((x array))
+ (assert (matrix-array-p x))
+ (array-dimension x 0)))
+
+
+(defgeneric matrix-column-number (m)
+ (:method ((x matrix)) (array-dimension (matrix-data x) 1))
+ (:method ((x vector)) (length x))
+ (:method ((x array))
+ (assert (matrix-array-p x))
+ (array-dimension x 1)))
+
+
+(defgeneric matrix-dimensions (m)
+ (:method ((x matrix)) (array-dimensions (matrix-data x)))
+ (:method ((x array))
+ (assert (matrix-array-p x))
+ (array-dimensions x)))
+
+
+(defgeneric matrix-dimension (m axis)
+ (:method ((m matrix) axis)
+ (declare (type fixnum axis))
+ (array-dimension (matrix-data m) axis))
+ (:method ((m array) axis)
+ (declare (type fixnum axis))
+ (array-dimension m axis)))
+
+
+(defmethod matrix-dimension :before ((m array) axis)
+ (assert (matrix-array-p m)))
+
+
+
+(defmethod matrix-element-type ((m matrix)) (array-element-type (matrix-data m)))
+
+(defmethod matrix-element-type ((m array)) (array-element-type m))
+
+(defmethod matrix-element-type :before ((m array))
+ (assert (matrix-array-p m)))
+
+
+(defgeneric shape-equal-p (m1 m2)
+ (:documentation "Returns non-NIL when M1 and M2 have the same `shape'."))
+
+
+(defmethod shape-equal-p ((m1 matrix) (m2 matrix))
+ (with-slots ((m1d data))
+ m1
+ (with-slots ((m2d data))
+ m2
+ (and (cl:= (array-dimension m1d 0)
+ (array-dimension m2d 0))
+ (cl:= (array-dimension m1d 1)
+ (array-dimension m2d 1))))))
+
+
+(defmethod shape-equal-p ((m1 array) (m2 matrix))
+ (assert (matrix-array-p m1 :check-element-type nil))
+ (let ((m1d m1))
+ (with-slots ((m2d data))
+ m2
+ (and (cl:= (array-dimension m1d 0)
+ (array-dimension m2d 0))
+ (cl:= (array-dimension m1d 1)
+ (array-dimension m2d 1))))))
+
+
+(defmethod shape-equal-p ((m2 matrix) (m1 array))
+ (shape-equal-p m1 m2))
+
+
+(defmethod shape-equal-p ((m1 array) (m2 array))
+ (assert (matrix-array-p m1 :check-element-type nil))
+ (assert (matrix-array-p m2 :check-element-type nil))
+ (let ((m1d m1)
+ (m2d m2)
+ )
+ (and (cl:= (array-dimension m1d 0)
+ (array-dimension m2d 0))
+ (cl:= (array-dimension m1d 1)
+ (array-dimension m2d 1)))))
+
+
+(defmethod shape-equal-p ((m1 vector) (m2 vector))
+ (cl:= (cl:length m1) (cl:length m2)))
+
+
+(defmethod shape-equal-p ((m1 column) (m2 vector))
+ (cl:= 1 (vector-length m1) (cl:length m2)))
+
+
+(defmethod shape-equal-p ((m2 vector) (m1 column))
+ (shape-equal-p m1 m2))
+
+
+;;;---------------------------------------------------------------------------
+;;; Creation.
+
+(defun make-matrix (n m
+ &rest array-keys
+ &key
+ (element-type *default-matrix-element-type*)
+ (initial-element (coerce 0 element-type) ie-supplied-p)
+ (initial-contents () ic-supplied-p)
+ &allow-other-keys)
+ (when (and ie-supplied-p ic-supplied-p)
+ (error "The MAKE-MATRIX :INITIAL-ELEMENT and :INITIAL-CONTENTS keywords may not be simultaneously supplied."))
+
+ (let ((data (cond ((and ic-supplied-p
+ (matrix-array-p initial-contents))
+ initial-contents)
+ (ic-supplied-p
+ (apply #'make-array (list n m)
+ (append array-keys
+ (list :element-type element-type))))
+ (t
+ (apply #'make-array (list n m)
+ (append array-keys
+ (list :element-type element-type
+ :initial-element initial-element)))
+ )))
+ )
+ (make-instance 'matrix :data data)))
+
+
+(defun make-identity-matrix (n &rest array-keys &key &allow-other-keys)
+ (when (getf array-keys :initial-contents)
+ (warn "MAKE-IDENTITY-MATRIX passed an :INITIAL-CONTENTS argument.")
+ (remf array-keys :initial-contents))
+ (when (getf array-keys :initial-element)
+ (warn "MAKE-IDENTITY-MATRIX passed an :INITIAL-ELEMENT argument.")
+ (remf array-keys :initial-element))
+ (let* ((id (apply #'make-matrix n n array-keys))
+ (data (matrix-data id))
+ (one (coerce 1 (array-element-type data)))
+ )
+ (dotimes (i n id)
+ (setf (aref data i i) one))))
+
+
+(defun eye (n &rest array-keys &key &allow-other-keys)
+ (apply #'make-identity-matrix n array-keys))
+
+
+(defun make-zero-matrix (n m
+ &rest array-keys
+ &key (element-type *default-matrix-element-type*)
+ &allow-other-keys)
+ (when (getf array-keys :initial-element)
+ (warn "MAKE-ZERO-MATRIX passed and :INITIAL-ELEMENT argument.")
+ (remf array-keys :initial-element))
+ (apply #'make-matrix n m (append array-keys
+ (list :initial-element (coerce 0 element-type)
+ :element-type element-type))))
+
+
+(defun zeros (n m &rest array-keys &key &allow-other-keys)
+ (apply #'make-zero-matrix n m array-keys))
+
+
+;;;---------------------------------------------------------------------------
+;;; copy-matrix --
+;;; Uses the code from K Pitman and B. Margolin appeared on CLL a long
+;;; long time ago.
+;;; Note that at least on LW this code has some of the usual problems
+;;; with DOUBLE-FLOATS.
+#|
+(defun copy-array (array)
+ (adjust-array (make-array (array-dimensions array)
+ :displaced-to array
+ :element-type (array-element-type array))
+ (array-dimensions array)
+ :displaced-to nil))
+|#
+
+
+(defmethod copy-matrix ((m matrix))
+ (let* ((array (matrix-data m))
+ (result-data
+ (adjust-array (make-array (array-dimensions array)
+ :displaced-to array
+ :element-type (array-element-type array))
+ (array-dimensions array)
+ :displaced-to nil))
+ )
+ (make-instance 'matrix :data result-data)))
+
+
+(defmethod copy-matrix ((m array))
+ (check-type m matrix-array)
+ (adjust-array (make-array (array-dimensions m)
+ :displaced-to m
+ :element-type (array-element-type m))
+ (array-dimensions m)
+ :displaced-to nil))
+
+
+;;;---------------------------------------------------------------------------
+;;; Equality.
+
+(defmethod =%2 ((x array) (y array))
+ (and (matrix-array-p x)
+ (matrix-array-p y)
+ ;; (cl:= (matrix-row-number x) (matrix-row-number y))
+ ;; (cl:= (matrix-column-number x) (matrix-column-number y))
+ (cl:= (array-dimension x 0) (array-dimension y 0))
+ (cl:= (array-dimension x 1) (array-dimension y 1))
+ (loop for i from 0 below (array-total-size x)
+ always (=%2 (row-major-aref x i) (row-major-aref y i)))))
+
+
+(defmethod =%2 ((x matrix) (y matrix))
+ (or (eq x y)
+ (=%2 (matrix-data x) (matrix-data y))))
+
+
+(defmethod =%2 ((x matrix) (y array))
+ (=%2 (matrix-data x) y))
+
+
+(defmethod =%2 ((y array) (x matrix))
+ (=%2 (matrix-data x) y))
+
+
+;;;---------------------------------------------------------------------------
+;;; zerop
+
+(defmethod zerop ((x array))
+ (and (matrix-array-p x)
+ (loop for i from 0 below (array-total-size x)
+ always (zerop (row-major-aref x i)))))
+
+
+(defmethod zerop ((x matrix))
+ (zerop (matrix-data x)))
+
+
+;;;---------------------------------------------------------------------------
+;;; Addition.
+
+(defmethod +%2 :before ((x matrix) (y matrix)
+ &optional
+ (r nil r-supplied-p))
+ (assert (shape-equal-p x y))
+ (when r-supplied-p
+ (assert (matrix-p x))
+ (assert (shape-equal-p x r))))
+
+
+(defmethod +%2 :before ((x matrix) (y array)
+ &optional
+ (r nil r-supplied-p))
+ (assert (shape-equal-p x y))
+ (when r-supplied-p
+ (assert (matrix-p x))
+ (assert (shape-equal-p x r))))
+
+
+(defmethod +%2 :before ((y array) (x matrix)
+ &optional
+ (r nil r-supplied-p))
+ (assert (shape-equal-p x y))
+ (when r-supplied-p
+ (assert (matrix-p x))
+ (assert (shape-equal-p x r))))
+
+
+(defmethod +%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (when r-supplied-p (assert (shape-equal-p y r)))
+ (with-slots (data) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size data) r)
+ (setf (row-major-aref result i) (+%2 x (row-major-aref data i))))
+ )))
+
+
+(defmethod +%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p))
+ (when r-supplied-p (assert (shape-equal-p y r)))
+ (with-slots (data) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size data) r)
+ (setf (row-major-aref result i) (+%2 x (row-major-aref data i))))
+ )))
+
+
+
+#+with-slots-slow
+(defmethod +%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (with-slots ((m1 data)) x
+ (with-slots ((m2 data)) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size m1))
+ (setf (row-major-aref result i)
+ (cl:+ (row-major-aref m1 i) (row-major-aref m2 i))))
+ )))
+ r
+ )
+
+
+(defmethod +%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (let ((m1 (matrix-data x))
+ (m2 (matrix-data y))
+ (result (matrix-data r))
+ )
+ (dotimes (i (array-total-size m1) r)
+ (setf (row-major-aref result i)
+ (cl:+ (row-major-aref m1 i) (row-major-aref m2 i))))
+ ))
+
+
+#+with-slots-slow
+(defmethod +%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (let ((m1 x))
+ (with-slots ((m2 data)) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size m1))
+ (setf (row-major-aref result i)
+ (cl:+ (row-major-aref m1 i) (row-major-aref m2 i))))
+ )))
+ r)
+
+
+
+(defmethod +%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (let ((m1 x)
+ (m2 (matrix-data y))
+ (result (matrix-data r))
+ )
+ (dotimes (i (array-total-size m1) r)
+ (setf (row-major-aref result i)
+ (cl:+ (row-major-aref m1 i) (row-major-aref m2 i))))
+ ))
+
+
+#+with-slots-slow
+(defmethod +%2 ((y matrix) (x array) &optional (r (copy-matrix y) r-supplied-p))
+ (let ((m1 x))
+ (with-slots ((m2 data)) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size m1))
+ (setf (row-major-aref result i)
+ (cl:+ (row-major-aref m1 i) (row-major-aref m2 i))))
+ )))
+ r)
+
+
+
+(defmethod +%2 ((y matrix) (x array) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (let ((m1 x)
+ (m2 (matrix-data y))
+ (result (matrix-data r))
+ )
+ (dotimes (i (array-total-size m1) r)
+ (setf (row-major-aref result i)
+ (cl:+ (row-major-aref m1 i) (row-major-aref m2 i))))
+ ))
+
+
+;;;---------------------------------------------------------------------------
+;;; Subtraction.
+
+
+(defmethod -%2 :before ((x matrix) (y matrix)
+ &optional
+ (r nil r-supplied-p))
+ (assert (shape-equal-p x y))
+ (when r-supplied-p
+ (assert (matrix-p x))
+ (assert (shape-equal-p x r))))
+
+
+(defmethod -%2 :before ((x matrix) (y array)
+ &optional
+ (r nil r-supplied-p))
+ (assert (shape-equal-p x y))
+ (when r-supplied-p
+ (assert (matrix-p x))
+ (assert (shape-equal-p x r))))
+
+
+(defmethod -%2 :before ((y array) (x matrix)
+ &optional
+ (r nil r-supplied-p))
+ (assert (shape-equal-p x y))
+ (when r-supplied-p
+ (assert (matrix-p x))
+ (assert (shape-equal-p x r))))
+
+
+(defmethod -%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (with-slots ((m1 data)) x
+ (with-slots ((m2 data)) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size m1))
+ (setf (row-major-aref result i)
+ (cl:- (row-major-aref m1 i) (row-major-aref m2 i))))
+ )))
+ r
+ )
+
+(defmethod -%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (let ((m1 x))
+ (with-slots ((m2 data)) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size m1))
+ (setf (row-major-aref result i)
+ (cl:- (row-major-aref m1 i) (row-major-aref m2 i))))
+ )))
+ r)
+
+(defmethod -%2 ((y matrix) (x array) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (with-slots ((m1 data)) y
+ (let ((m2 x))
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size m1))
+ (setf (row-major-aref result i)
+ (cl:- (row-major-aref m1 i) (row-major-aref m2 i))))
+ )))
+ r)
+
+
+;;;---------------------------------------------------------------------------
+;;; Multiplication.
+
+(defgeneric make-*-result-matrix (x y &optional element-type)
+ (:documentation
+ "Creates a matrix serve as result of a multiplication.")
+
+ (:method ((x matrix) (y matrix)
+ &optional (element-type *default-matrix-element-type*))
+ (make-matrix (matrix-row-number x)
+ (matrix-column-number y)
+ :element-type element-type))
+
+ (:method ((x array) (y matrix)
+ &optional (element-type *default-matrix-element-type*))
+ (make-matrix (matrix-row-number x)
+ (matrix-column-number y)
+ :element-type element-type))
+
+ (:method ((x matrix) (y array)
+ &optional (element-type *default-matrix-element-type*))
+ (make-matrix (matrix-row-number x)
+ (matrix-column-number y)
+ :element-type element-type))
+
+ (:method ((x vector) (y matrix)
+ &optional (element-type *default-matrix-element-type*))
+ (make-matrix 1
+ (matrix-column-number y)
+ :element-type element-type))
+
+ (:method ((x matrix) (y vector)
+ &optional (element-type *default-matrix-element-type*))
+ (make-matrix (matrix-row-number x)
+ 1
+ :element-type element-type)))
+
+
+
+(defgeneric conforming-*-dimensions-p (x y r)
+ (:documentation
+ "Checks whether the dimensions of R are correct to hold XxY.")
+
+ (:method ((x t) (y t) (r t))
+ (error "Result object is not a matrix: ~S." r))
+
+ ;; 1. Matrix results.
+
+ ;; 1.1. Number methods.
+
+ (:method ((x number) (y matrix) (r matrix))
+ (shape-equal-p y r))
+
+ (:method ((x number) (y array) (r matrix))
+ (shape-equal-p y r))
+
+ (:method ((y matrix) (x number) (r matrix))
+ (conforming-*-dimensions-p x y r))
+
+ (:method ((y array) (x number) (r matrix))
+ (conforming-*-dimensions-p x y r))
+
+ ;; 1.2. Matrix and Array arguments.
+
+ (:method ((x array) (y array) (r matrix))
+ (and (matrix-array-p x)
+ (matrix-array-p y)
+ (cl:= (matrix-row-number x) (matrix-column-number y))
+ (cl:= (matrix-row-number x) (matrix-row-number r))
+ (cl:= (matrix-column-number y) (matrix-column-number r))))
+
+ (:method ((x vector) (y array) (r matrix))
+ (and (matrix-array-p y)
+ (cl:= (matrix-row-number x) (matrix-column-number y))
+ (cl:= (matrix-row-number x) (matrix-row-number r))
+ (cl:= (matrix-column-number y) (matrix-column-number r))))
+
+
+ (:method ((x array) (y vector) (r matrix))
+ (and (matrix-array-p x)
+ (cl:= (matrix-row-number x) (matrix-column-number y))
+ (cl:= (matrix-row-number x) (matrix-row-number r))
+ (cl:= (matrix-column-number y) (matrix-column-number r))))
+
+ (:method ((x matrix) (y matrix) (r matrix))
+ (conforming-*-dimensions-p (matrix-data x) (matrix-data y) r))
+
+ (:method ((x array) (y matrix) (r matrix))
+ (conforming-*-dimensions-p x (matrix-data y) r))
+
+ (:method ((x matrix) (y array) (r matrix))
+ (conforming-*-dimensions-p (matrix-data x) y r))
+
+
+ ;; 2. Array results.
+
+ ;; 2.1. Number methods.
+
+ (:method ((x number) (y matrix) (r array))
+ (shape-equal-p y r))
+
+ (:method ((x number) (y array) (r array))
+ (shape-equal-p y r))
+
+ (:method ((y matrix) (x number) (r array))
+ (conforming-*-dimensions-p x y r))
+
+ (:method ((y array) (x number) (r array))
+ (conforming-*-dimensions-p x y r))
+
+ ;; 2.2. Matrix and Array arguments.
+
+ (:method ((x array) (y array) (r array))
+ (and (matrix-array-p x)
+ (matrix-array-p y)
+ (matrix-array-p r)
+ (cl:= (matrix-row-number x) (matrix-column-number y))
+ (cl:= (matrix-row-number x) (matrix-row-number r))
+ (cl:= (matrix-column-number y) (matrix-column-number r))))
+
+ (:method ((x vector) (y array) (r array))
+ (and (matrix-array-p y)
+ (cl:= (matrix-row-number x) (matrix-column-number y))
+ (cl:= (matrix-row-number x) (matrix-row-number r))
+ (cl:= (matrix-column-number y) (matrix-column-number r))))
+
+ (:method ((x array) (y vector) (r array))
+ (and (matrix-array-p x)
+ (cl:= (matrix-row-number x) (matrix-column-number y))
+ (cl:= (matrix-row-number x) (matrix-row-number r))
+ (cl:= (matrix-column-number y) (matrix-column-number r))))
+
+ (:method ((x matrix) (y matrix) (r array))
+ (conforming-*-dimensions-p (matrix-data x) (matrix-data y) r))
+
+ (:method ((x array) (y matrix) (r array))
+ (conforming-*-dimensions-p x (matrix-data y) r))
+
+ (:method ((x matrix) (y array) (r array))
+ (conforming-*-dimensions-p (matrix-data x) y r))
+
+ ;; 3. Other methods.
+
+ (:method ((x vector) (y vector) (r t))
+ ;; Maybe a warning here is needed.
+ (cl:= (length x) (length y)))
+ )
+
+
+;;; *%2 Methods.
+
+(defmethod *%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (when r-supplied-p (assert (conforming-*-dimensions-p x y r)))
+ (with-slots ((m data)) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size m) r)
+ (setf (row-major-aref result i) (*%2 x (row-major-aref m i))))
+ )))
+
+
+(defmethod *%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (*%2 x y r))
+
+
+(defmethod *%2 ((x matrix) (y matrix)
+ &optional (r (make-*-result-matrix x y) r-supplied-p))
+ (when r-supplied-p (assert (conforming-*-dimensions-p x y r)))
+ ;; X is a NxK matrix
+ ;; Y is a KxM matrix
+ ;; R is a NxM matrix
+ (let* ((m1 (matrix-data x))
+ (m2 (matrix-data y))
+ (result (matrix-data r))
+ (n (array-dimension m1 0))
+ (k (array-dimension m1 1))
+ (m (array-dimension m2 0))
+ (z (coerce 0 (array-element-type result)))
+ )
+ (declare (type fixnum n k m))
+ (dotimes (i n r)
+ (dotimes (j m)
+ (let ((temp z))
+ (dotimes (l k)
+ (setf (aref result i j)
+ (setf temp (+%2 temp (*%2 (aref m1 i l) (aref m2 l j)))))
+ ))
+ ))))
+
+
+(defmethod *%2 ((x array) (y matrix)
+ &optional (r (make-*-result-matrix x y) r-supplied-p))
+ (when r-supplied-p (assert (conforming-*-dimensions-p x y r)))
+ ;; X is a NxK matrix
+ ;; Y is a KxM matrix
+ ;; R is a NxM matrix
+ (let* ((m1 x)
+ (m2 (matrix-data y))
+ (result (matrix-data r))
+ (n (array-dimension m1 0))
+ (k (array-dimension m1 1))
+ (m (array-dimension m2 0))
+ (z (coerce 0 (array-element-type result)))
+ )
+ (declare (type fixnum n k m))
+ (dotimes (i n r)
+ (dotimes (j m)
+ (let ((temp z))
+ (dotimes (l k)
+ (setf (aref result i j)
+ (incf temp (* (aref m1 i l) (aref m2 l j))))
+ ))
+ ))))
+
+
+(defmethod *%2 ((x matrix) (y array)
+ &optional (r (make-*-result-matrix x y) r-supplied-p))
+ (when r-supplied-p (assert (conforming-*-dimensions-p x y r)))
+ ;; X is a NxK matrix
+ ;; Y is a KxM matrix
+ ;; R is a NxM matrix
+ (let* ((m1 (matrix-data x))
+ (m2 y)
+ (result (matrix-data r))
+ (n (array-dimension m1 0))
+ (k (array-dimension m1 1))
+ (m (array-dimension m2 1))
+ (z (coerce 0 (array-element-type result)))
+ )
+ (declare (type fixnum n k m))
+ (dotimes (i n r)
+ (dotimes (j m)
+ (let ((temp z))
+ (dotimes (l k)
+ (setf (aref result i j)
+ (incf temp (* (aref m1 i l) (aref m2 l j))))
+ ))
+ ))))
+
+
+(defmethod *%2 ((x vector) (y matrix)
+ &optional (r (make-*-result-matrix x y) r-supplied-p))
+ (when r-supplied-p (assert (conforming-*-dimensions-p x y r)))
+ (let* ((m1 x)
+ (m2 (matrix-data y))
+ (result (matrix-data r))
+ (n 1)
+ (k (length m1))
+ (m (array-dimension m2 1))
+ (z (coerce 0 (array-element-type result)))
+ )
+ (declare (type fixnum k m) (type (integer 1 1) n))
+ (dotimes (i n r)
+ (dotimes (j m)
+ (let ((temp z))
+ (dotimes (l k)
+ (setf (aref result i j)
+ (incf temp (* (aref m1 l) (aref m2 l j))))
+ ))
+ ))))
+
+
+(defmethod *%2 ((x matrix) (y vector)
+ &optional (r (make-*-result-matrix x y) r-supplied-p))
+ (when r-supplied-p (assert (conforming-*-dimensions-p x y r)))
+ (let* ((m1 y )
+ (m2 (matrix-data x))
+ (result (matrix-data r))
+ (n (array-dimension m1 0))
+ (k (length y))
+ (m k)
+ (z (coerce 0 (array-element-type result)))
+ )
+ (declare (type fixnum n k m))
+ (dotimes (i n r)
+ (dotimes (j m)
+ (let ((temp z))
+ (dotimes (l k)
+ (setf (aref result i j)
+ (incf temp (* (aref m1 l) (aref m2 l j))))
+ ))
+ ))))
+
+
+(defmethod *%2 ((x number) (y array) &optional (r (copy-matrix y) r-supplied-p))
+ (when r-supplied-p (assert (conforming-*-dimensions-p x y r)))
+ (let ((m y)
+ (result r)
+ )
+ (dotimes (i (array-total-size m) r)
+ (setf (row-major-aref result i) (*%2 x (row-major-aref m i))))
+ ))
+
+
+(defmethod *%2 ((y array) (x number) &optional (r (copy-matrix y) r-supplied-p))
+ (*%2 x y r))
+
+
+;;; The next one breaks the return type convention.
+#| Defined in 'vector.lisp'.
+(defmethod *%2 ((x vector) (y vector) &optional r)
+ (declare (ignore r))
+ (assert (conforming-*-dimensions-p x y nil))
+ (let ((result 0))
+ (dotimes (i (length x) result)
+ (setf result (+%2 result (* (aref x i) (aref y i)))))))
+|#
+
+
+;;;---------------------------------------------------------------------------
+;;; Division.
+;;; Only the simple form of division is implemented here.
+;;; The general form requires the inverse, which requires extra machinery
+
+(defmethod /%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (*%2 (/ x) y r))
+
+
+(defmethod /%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (*%2 (/ x) y r))
+
+
+;;;---------------------------------------------------------------------------
+;;; Exponentiation.
+;;; ... missing ...
+
+
+;;;---------------------------------------------------------------------------
+;;; Element-wise operations.
+
+(defmethod .*%2 ((x number) (y matrix)
+ &optional (r (copy-matrix y) r-supplied-p))
+ (when r-supplied-p (assert (shape-equal-p y r)))
+ (with-slots (data) y
+ (with-slots ((result data)) r
+ (dotimes (i (array-total-size data) r)
+ (setf (row-major-aref result i) (*%2 x (row-major-aref data i))))
+ )))
+
+
+(defmethod .*%2 ((y matrix) (x number) &optional (r (copy-matrix y)))
+ (.*%2 x y r))
+
+
+(defmethod .*%2 ((x number) (y array)
+ &optional (r (copy-matrix y) r-supplied-p))
+ (assert (matrix-array-p y))
+ (when r-supplied-p (assert (shape-equal-p y r)))
+ (let* ((data y)
+ (result r)
+ )
+ (dotimes (i (array-total-size data) r)
+ (setf (row-major-aref result i) (*%2 x (row-major-aref data i))))
+ ))
+
+
+(defmethod .*%2 ((y array) (x number) &optional (r (copy-matrix y)))
+ (.*%2 x y r))
+
+
+
+(defmethod .*%2 ((x array) (y array) &optional (r (copy-matrix y) r-supplied-p))
+ (when r-supplied-p (assert (shape-equal-p y r)))
+ (assert (shape-equal-p x y))
+
+ (dotimes (i (array-total-size y) r)
+ (setf (row-major-aref r i)
+ (*%2 (row-major-aref x i) (row-major-aref y i)))))
+
+
+;;;---------------------------------------------------------------------------
+;;; Transpose
+
+(defgeneric conforming-transpose-dimensions-p (x r)
+ (:documentation
+ "Checks whether the dimensions of R are correct to hold X'.")
+
+ (:method ((x t) (r t))
+ (error "Result object is not a matrix: ~S." r))
+
+ ;; 1. Matrix results.
+
+ (:method ((x array) (r matrix))
+ (and (matrix-array-p x)
+ (cl:= (matrix-row-number x) (matrix-column-number r))
+ (cl:= (matrix-column-number x) (matrix-row-number r))))
+
+ (:method ((x matrix) (r matrix))
+ (and (cl:= (matrix-row-number x) (matrix-column-number r))
+ (cl:= (matrix-column-number x) (matrix-row-number r))))
+
+
+ ;; 2. Array results.
+
+ (:method ((x array) (r array))
+ (and (matrix-array-p x)
+ (matrix-array-p r)
+ (cl:= (matrix-row-number x) (matrix-column-number r))
+ (cl:= (matrix-column-number x) (matrix-row-number r))))
+
+ (:method ((x matrix) (r array))
+ (and (matrix-array-p r)
+ (cl:= (matrix-row-number x) (matrix-column-number r))
+ (cl:= (matrix-column-number x) (matrix-row-number r))))
+ )
+
+
+
+(defmethod transpose ((m matrix)
+ &optional
+ (r (make-matrix (matrix-column-number m)
+ (matrix-row-number m)
+ :element-type
+ (matrix-element-type m))
+ r-supplied-p))
+ (when r-supplied-p (assert (conforming-transpose-dimensions-p m r)))
+ (let ((rows (matrix-row-number m))
+ (cols (matrix-column-number m))
+ (m (matrix-data m))
+ (result (matrix-data r))
+ )
+ (dotimes (i rows r)
+ (dotimes (j cols)
+ (setf (aref result j i) (aref m i j))))))
+
+
+(defmethod transpose ((m array)
+ &optional
+ (r (make-array (array-dimensions m)
+ :element-type
+ (array-element-type m)
+ )
+ r-supplied-p))
+ (when r-supplied-p (assert (conforming-transpose-dimensions-p m r)))
+ (let ((rows (matrix-row-number m))
+ (cols (matrix-column-number m))
+ )
+ (dotimes (i rows r)
+ (dotimes (j cols)
+ (setf (aref r j i) (aref m i j))))))
+
+
+(defmethod transpose :before ((m array) &optional r)
+ (declare (ignore r))
+ (assert (matrix-array-p m)))
+
+
+;;;---------------------------------------------------------------------------
+;;; Utilities.
+
+#+simple
+(defmethod print-object ((m matrix) stream)
+ (print-unreadable-object (m stream :identity t :type t)
+ (format stream "~S" (matrix-data m))))
+
+
+(defmethod print-object ((m matrix) stream)
+ (let ((d (matrix-data m)))
+ (format stream "#M(~S ~S ~S)"
+ (array-element-type d)
+ (array-dimensions d)
+ d)))
+
+;;; Missing:
+;;; Define #M read macro.
+
+
+;;;; end of file -- matrix.lisp --
Added: trunk/common-math/numerics/linear-algebra/vector-matrix-conditions.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/vector-matrix-conditions.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,27 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; vector-matrix-conditions.lisp --
+
+(in-package "CL.MATH.LINEAR-ALGEBRA")
+
+;;;---------------------------------------------------------------------------
+;;; Conditions relative to the vectors and matrix subsystem..
+
+(define-condition incompatible-shapes (error)
+ ((a :accessor incompatible-shape-a
+ :initarg :a)
+ (b :accessor incompatible-shape-b
+ :initarg :b)
+ )
+ (:documentation "The Incompatible Shapes Error.
+This error is signalled by functions that get arguments -- mainly
+matrices and/or vectors -- whose 'shape' is not compatible with
+respect to the semantics of the operation.")
+
+ (:report (lambda (isc stream)
+ (declare (ignore isc))
+ (format stream "The shapes passed have incompatible dimensions.")))
+ )
+
+
+;;;; end of file -- vector-matrix-conditions.lisp --
Added: trunk/common-math/numerics/linear-algebra/vector-matrix.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/vector-matrix.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,613 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; vector-matrix.lisp --
+;;;; Protocols and functions common to and dependent on vector and matrix
+;;;; definitions.
+
+(in-package "CL.MATH.LINEAR-ALGEBRA")
+
+;;;;===========================================================================
+;;;; Protocol.
+
+;;; join and stack
+;;; Note: no 'result' machinery in place yet.
+
+(defgeneric join%2 (d1 d2))
+
+(defgeneric stack%2 (d1 d2))
+
+
+;;; slice
+
+(defgeneric slice (e shape-designator))
+
+
+;;; range
+
+(defstruct (range (:constructor %range (start &optional end step)))
+ (start 0 :type (or fixnum (member * cl:*)))
+ (end nil :type (or null fixnum (member * cl:*)))
+ (step 1 :type fixnum))
+
+
+(defstruct (decreasing-range (:include range)))
+
+(defstruct (increasing-range (:include range)))
+
+
+(defun range (start &optional end (step 1))
+ (if end
+ (if (<= start end)
+ (make-increasing-range :start start :end end :step step)
+ (make-decreasing-range :start start :end end :step step))
+ (if (plusp step)
+ (make-increasing-range :start start :end end :step step)
+ (make-decreasing-range :start start :end end :step step))
+ ))
+
+
+
+;;; ref
+
+(defgeneric ref (e index &rest indexes))
+
+(defgeneric ref%2 (e i j))
+
+
+;;;---------------------------------------------------------------------------
+;;; Implementation.
+
+;;; join%2 --
+
+(defmethod join%2 ((d1 vector) (d2 vector))
+ (concatenate 'vector d1 d2))
+
+
+(defmethod join%2 ((d1 string) (d2 string))
+ (concatenate 'string d1 d2))
+
+
+(defmethod join%2 ((d1 list) (d2 list))
+ (concatenate 'list d1 d2))
+
+
+(defmethod join%2 ((d1 vector) (d2 sequence))
+ (concatenate 'vector d1 d2))
+
+
+(defmethod join%2 ((d1 sequence) (d2 vector))
+ (concatenate (type-of d1) d1 d2))
+
+
+(defmethod join%2 ((d1 column) (d2 column))
+ (let ((initial-contents
+ (map 'list #'vector (column-vector d1) (column-vector d2))))
+ (declare (dynamic-extent initial-contents))
+ (make-array (list (length (column-vector d1)) 2)
+ :initial-contents initial-contents)))
+
+
+(defmethod join%2 ((d1 vector) (d2 column))
+ (unless (cl:= (length (column-vector d2)) 1)
+ (error "Trying to JOIN a vector (1x~D) to a column (~Dx1)."
+ (length d1)
+ (length (column-vector d2))))
+ (concatenate 'vector d1 (column-vector d2)))
+
+
+(defmethod join%2 ((d1 column) (d2 vector))
+ (unless (cl:= (length (column-vector d1)) 1)
+ (error "Trying to JOIN a column (~Dx1) to a vector (1x~D)."
+ (length (column-vector d1))
+ (length d2)))
+ (concatenate 'vector (column-vector d1) d2))
+
+
+(defmethod join%2 ((d1 column) (d2 array))
+ (let ((nr (array-dimension d2 0)))
+ (assert (cl:= (length (column-vector d1)) nr))
+ (let* ((nc (1+ (array-dimension d2 1)))
+ (result (make-array (list nr nc)
+ :initial-element (aref d2 0 0)
+ :element-type (array-element-type d2)))
+ (col (column-vector d1))
+ )
+ (dotimes (i nr result)
+ (setf (aref result i 0) (aref col i))
+ (loop for j from 1 below nc
+ do (setf (aref result i j) (aref d2 i (1- j)))))
+ )))
+
+
+(defmethod join%2 ((d1 array) (d2 column))
+ (let ((nr (array-dimension d1 0)))
+ (assert (cl:= (length (column-vector d2)) nr))
+ (let* ((nc (1+ (array-dimension d1 1)))
+ (nc-1 (1- nc))
+ (result (make-array (list nr nc)
+ :initial-element (aref d1 0 0)
+ :element-type (array-element-type d1)))
+ (col (column-vector d2))
+ )
+
+ (dotimes (i nr result)
+ (setf (aref result i nc-1) (aref col i))
+ (dotimes (j nc-1)
+ (setf (aref result i j) (aref d1 i j))))
+ )))
+
+
+(defmethod join%2 ((d1 array) (d2 array))
+ (let ((nr1 (array-dimension d1 0))
+ (nr2 (array-dimension d2 0))
+ )
+ (assert (cl:= nr1 nr2))
+ (let* ((nc1 (array-dimension d1 1))
+ (nc2 (array-dimension d2 1))
+ (result (make-array (list nr1 (+ nc1 nc2))
+ :initial-element (aref d1 0 0)
+ :element-type (array-element-type d1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref d1 i j))))
+
+ (dotimes (i nr2 result)
+ (dotimes (j nc2)
+ (setf (aref result i (+ j nc1)) (aref d2 i j))))
+ )))
+
+
+(defmethod join%2 ((d1 matrix) (d2 matrix))
+ (let ((nr1 (matrix-row-number d1))
+ (nr2 (matrix-row-number d2))
+ )
+ (assert (cl:= nr1 nr2))
+ (let* ((nc1 (matrix-column-number d1))
+ (nc2 (matrix-column-number d2))
+ (data1 (matrix-data d1))
+ (data2 (matrix-data d2))
+ (result (make-array (list nr1 (+ nc1 nc2))
+ :initial-element (aref data1 0 0)
+ :element-type (array-element-type data1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref data1 i j))))
+
+ (dotimes (i nr2 (make-instance 'matrix :data result))
+ (dotimes (j nc2)
+ (setf (aref result i (+ j nc1)) (aref data2 i j))))
+ )))
+
+
+(defmethod join%2 ((d1 array) (d2 matrix))
+ (let ((nr1 (matrix-row-number d1))
+ (nr2 (matrix-row-number d2))
+ )
+ (assert (cl:= nr1 nr2))
+ (let* ((nc1 (matrix-column-number d1))
+ (nc2 (matrix-column-number d2))
+ (data1 d1)
+ (data2 (matrix-data d2))
+ (result (make-array (list nr1 (+ nc1 nc2))
+ :initial-element (aref data1 0 0)
+ :element-type (array-element-type data1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref data1 i j))))
+
+ (dotimes (i nr2 (make-instance 'matrix :data result))
+ (dotimes (j nc2)
+ (setf (aref result i (+ j nc1)) (aref data2 i j))))
+ )))
+
+
+(defmethod join%2 ((d1 matrix) (d2 array))
+ (let ((nr1 (matrix-row-number d1))
+ (nr2 (matrix-row-number d2))
+ )
+ (assert (cl:= nr1 nr2))
+ (let* ((nc1 (matrix-column-number d1))
+ (nc2 (matrix-column-number d2))
+ (data1 (matrix-data d1))
+ (data2 d2)
+ (result (make-array (list nr1 (+ nc1 nc2))
+ :initial-element (aref data1 0 0)
+ :element-type (array-element-type data1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref data1 i j))))
+
+ (dotimes (i nr2 (make-instance 'matrix :data result))
+ (dotimes (j nc2)
+ (setf (aref result i (+ j nc1)) (aref data2 i j))))
+ )))
+
+
+(defmethod join%2 ((d1 vector) (d2 matrix))
+ (let ((result (join%2 d1 (matrix-data d2))))
+ (make-instance 'matrix :data result)))
+
+
+(defmethod join%2 ((d2 matrix) (d1 vector))
+ (let ((result (join%2 (matrix-data d2) d1)))
+ (make-instance 'matrix :data result)))
+
+
+(defmethod join%2 ((d1 column) (d2 matrix))
+ (let ((result (join%2 d1 (matrix-data d2))))
+ (make-instance 'matrix :data result)))
+
+
+(defmethod join%2 ((d2 matrix) (d1 column))
+ (let ((result (join%2 (matrix-data d2) d1)))
+ (make-instance 'matrix :data result)))
+
+
+;;; stack%2 --
+
+(defmethod stack%2 ((d1 column) (d2 column))
+ (column (concatenate 'vector (column-vector d1) (column-vector d2))))
+
+
+(defmethod stack%2 ((d1 vector) (d2 array))
+ (let ((nc (array-dimension d2 1)))
+ (assert (cl:= (length d1) nc))
+ (let* ((nr (1+ (array-dimension d2 0)))
+ (result (make-array (list nr nc)
+ :initial-element (aref d2 0 0)
+ :element-type (array-element-type d2)))
+ )
+ (dotimes (i nc)
+ (setf (aref result 0 i) (aref d1 i)))
+ (loop for i from 1 below nr
+ do (dotimes (j nc)
+ (setf (aref result i j) (aref d2 (1- i) j))))
+ result)))
+
+
+(defmethod stack%2 ((d1 sequence) (d2 array))
+ (let ((nc (array-dimension d2 1)))
+ (assert (cl:= (length d1) nc))
+ (let* ((nr (1+ (array-dimension d2 0)))
+ (result (make-array (list nr nc)
+ :initial-element (aref d2 0 0)
+ :element-type (array-element-type d2)))
+ )
+ (dotimes (i nc)
+ (setf (aref result 0 i) (elt d1 i)))
+ (loop for i from 1 below nr
+ do (dotimes (j nc)
+ (setf (aref result i j) (aref d2 (1- i) j))))
+ result)))
+
+
+(defmethod stack%2 ((d1 array) (d2 vector))
+ (let ((nc (array-dimension d1 1)))
+ (assert (cl:= (length d2) nc))
+ (let* ((nr (1+ (array-dimension d1 0)))
+ (nr-1 (1- nr))
+ (result (make-array (list nr nc)
+ :initial-element (aref d1 0 0)
+ :element-type (array-element-type d1)))
+ )
+
+ (dotimes (i nr-1)
+ (dotimes (j nc)
+ (setf (aref result i j) (aref d1 i j))))
+ (dotimes (i nc)
+ (setf (aref result nr-1 i) (aref d2 i)))
+ result)))
+
+
+(defmethod stack%2 ((d1 array) (d2 sequence))
+ (let ((nc (array-dimension d1 1)))
+ (assert (cl:= (length d2) nc))
+ (let* ((nr (1+ (array-dimension d1 0)))
+ (nr-1 (1- nr))
+ (result (make-array (list nr nc)
+ :initial-element (aref d1 0 0)
+ :element-type (array-element-type d1)))
+ )
+
+ (dotimes (i nr-1)
+ (dotimes (j nc)
+ (setf (aref result i j) (aref d1 i j))))
+ (dotimes (i nc)
+ (setf (aref result nr-1 i) (elt d2 i)))
+ result)))
+
+
+(defmethod stack%2 ((d1 array) (d2 array))
+ (let ((nc1 (array-dimension d1 1))
+ (nc2 (array-dimension d2 1))
+ )
+ (assert (cl:= nc1 nc2))
+ (let* ((nr1 (array-dimension d1 0))
+ (nr2 (array-dimension d2 0))
+ (result (make-array (list (+ nr1 nr2) nc1)
+ :initial-element (aref d1 0 0)
+ :element-type (array-element-type d1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref d1 i j))))
+
+ (dotimes (i nr2 result)
+ (dotimes (j nc2)
+ (setf (aref result (+ i nr1) j) (aref d2 i j))))
+ )))
+
+
+(defmethod stack%2 ((d1 matrix) (d2 matrix))
+ (let ((nc1 (matrix-column-number d1))
+ (nc2 (matrix-column-number d2))
+ )
+ (assert (cl:= nc1 nc2))
+ (let* ((nr1 (matrix-row-number d1))
+ (nr2 (matrix-row-number d2))
+ (data1 (matrix-data d1))
+ (data2 (matrix-data d2))
+ (result (make-array (list (+ nr1 nr2) nc1)
+ :initial-element (aref data1 0 0)
+ :element-type (array-element-type data1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref data1 i j))))
+
+ (dotimes (i nr2 (make-instance 'matrix :data result))
+ (dotimes (j nc2)
+ (setf (aref result (+ i nr1) j) (aref data2 i j))))
+ )))
+
+
+(defmethod stack%2 ((d1 array) (d2 matrix))
+ (let ((nc1 (matrix-column-number d1))
+ (nc2 (matrix-column-number d2))
+ )
+ (assert (cl:= nc1 nc2))
+ (let* ((nr1 (matrix-row-number d1))
+ (nr2 (matrix-row-number d2))
+ (data1 d1)
+ (data2 (matrix-data d2))
+ (result (make-array (list (+ nr1 nr2) nc1)
+ :initial-element (aref data1 0 0)
+ :element-type (array-element-type data1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref data1 i j))))
+
+ (dotimes (i nr2 (make-instance 'matrix :data result))
+ (dotimes (j nc2)
+ (setf (aref result (+ i nr1) j) (aref data2 i j))))
+ )))
+
+
+
+(defmethod stack%2 ((d1 matrix) (d2 array))
+ (let ((nc1 (matrix-column-number d1))
+ (nc2 (matrix-column-number d2))
+ )
+ (assert (cl:= nc1 nc2))
+ (let* ((nr1 (matrix-row-number d1))
+ (nr2 (matrix-row-number d2))
+ (data1 (matrix-data d1))
+ (data2 d2)
+ (result (make-array (list (+ nr1 nr2) nc1)
+ :initial-element (aref data1 0 0)
+ :element-type (array-element-type data1)))
+ )
+
+ (dotimes (i nr1)
+ (dotimes (j nc1)
+ (setf (aref result i j) (aref data1 i j))))
+
+ (dotimes (i nr2 (make-instance 'matrix :data result))
+ (dotimes (j nc2)
+ (setf (aref result (+ i nr1) j) (aref data2 i j))))
+ )))
+
+
+(defmethod stack%2 ((d1 vector) (d2 matrix))
+ (let ((result (stack%2 d1 (matrix-data d2))))
+ (make-instance 'matrix :data result)))
+
+
+(defmethod stack%2 ((d2 matrix) (d1 vector))
+ (let ((result (stack%2 (matrix-data d2) d1)))
+ (make-instance 'matrix :data result)))
+
+
+(defmethod stack%2 ((d1 column) (d2 matrix))
+ (let ((result (stack%2 d1 (matrix-data d2))))
+ (make-instance 'matrix :data result)))
+
+
+(defmethod stack%2 ((d2 matrix) (d1 column))
+ (let ((result (stack%2 (matrix-data d2) d1)))
+ (make-instance 'matrix :data result)))
+
+
+;;; slice --
+#|
+(defmethod slice ((v vector) (x fixnum))
+ (aref v x))
+
+
+(defmethod slice ((v vector) (sd list))
+ (assert (<= (list-length sd) 2))
+ (destructuring-bind (start &optional end)
+ sd
+ (let ((start (cond ((or (eq start '*) (eq start 'cl:*)) 0)
+ ((typep start `(mod ,array-total-size-limit)) start)
+ (t (error "Unrecognized start ~S in slice designator."
+ start))))
+ (end (cond ((or (null end) (eq end '*) (eq end 'cl:*)) nil)
+ ((typep end `(mod ,array-total-size-limit)) end)
+ (t (error "Unrecognized end ~S in slice designator."
+ end))))
+ )
+ (subseq v start end))))
+
+
+(defmethod slice ((v column) (sd list))
+ (column (slice (column-vector v) sd)))
+
+
+(defmethod slice ((v vector) (sd (eql '*)))
+ (copy-vector v))
+
+
+(defmethod slice ((v vector) (sd (eql 'cl:*)))
+ (copy-vector v))
+
+
+(defmethod slice ((v column) (sd (eql '*)))
+ (copy-vector v))
+
+
+(defmethod slice ((v column) (sd (eql 'cl:*)))
+ (copy-vector v))
+
+
+(defmethod slice (()))
+ |#
+
+
+
+;;; ref --
+
+(defmethod ref ((v vector) (i fixnum) &rest indexes)
+ (when indexes
+ (error "REF of a vector requires only one index but ~D were supplied."
+ (1+ (list-length indexes))))
+
+
+ (aref v i))
+
+
+(defmethod ref ((v vector) (i increasing-range) &rest indexes)
+ (when indexes
+ (error "REF of a vector requires only one index or range but a range and ~D indexes were supplied."
+ (list-length indexes)))
+
+ (loop for vi from (range-start i)
+ below (or (range-end i) (length v))
+ by (range-step i)
+ collect (aref v vi) into vs
+ finally (return (coerce vs 'vector))))
+
+
+(defmethod ref ((v vector) (i decreasing-range) &rest indexes)
+ (when indexes
+ (error "REF of a vector requires only one index or range but a range and ~D indexes were supplied."
+ (list-length indexes)))
+
+ (loop for vi from (or (range-end i) (1- (length v)))
+ downto (range-start i)
+ by (abs (range-step i))
+ collect (aref v vi) into vs
+ finally (return (coerce vs 'vector))))
+
+
+(defmethod ref ((v vector) (i vector) &rest indexes)
+ (when indexes
+ (error "REF of a vector requires only one index or subset but a subset and ~D indexes were supplied."
+ (list-length indexes)))
+
+ (loop for x across i
+ collect (aref v x) into vs
+ finally (return (coerce vs 'vector))))
+
+
+(defmethod ref ((v vector) (i list) &rest indexes)
+ (when indexes
+ (error "REF of a vector requires only one index or subset but a subset and ~D indexes were supplied."
+ (list-length indexes)))
+
+ (loop for x in i
+ collect (aref v x) into vs
+ finally (return (coerce vs 'vector))))
+
+
+(defmethod ref ((v column) (i fixnum) &rest indexes)
+ (apply #'ref (column-vector v) i indexes))
+
+
+(defmethod ref ((v column) (i range) &rest indexes)
+ (column (apply #'ref (column-vector v) i indexes)))
+
+
+(defmethod ref ((v column) (i vector) &rest indexes)
+ (column (apply #'ref (column-vector v) i indexes)))
+
+
+(defmethod ref ((v column) (i list) &rest indexes)
+ (column (apply #'ref (column-vector v) i indexes)))
+
+
+(defmethod ref ((m matrix) (i t) &rest indexes)
+ (when (cl:/= 1 (list-length indexes))
+ (error "........"))
+
+ (ref%2 m i (first indexes)))
+
+
+
+(defmethod ref%2 ((m matrix) (i fixnum) (j fixnum))
+ (aref (matrix-data m) i j))
+
+
+(defmethod ref%2 ((m matrix) (i fixnum) (j vector))
+ (loop with md = (matrix-data m)
+ for x across j
+ collecting (aref md i x) into result
+ finally return (coerce result 'vector)))
+
+
+(defmethod ref%2 ((m matrix) (i vector) (j fixnum))
+ (loop with md = (matrix-data m)
+ for x across i
+ collecting (aref md x j) into result
+ finally return (coerce result 'vector)))
+
+
+(defmethod ref%2 ((m matrix) (i vector) (j vector))
+ (loop with md = (matrix-data m)
+ for x across i
+ collecting (loop for y across j
+ collecting (aref md x y))
+ into result
+ finally
+ (return (make-instance 'matrix
+ :data (make-array (list (length i)
+ (length j))
+ :initial-contents result)))))
+
+
+
+
+
+;;;;===========================================================================
+;;;; Generic Interface Functions.
+
+(defun join (&rest elements)
+ (cl:reduce #'join%2 elements))
+
+
+(defun stack (&rest elements)
+ (cl:reduce #'stack%2 elements))
+
+;;;; end of file -- vector-matrix.lisp --
Added: trunk/common-math/numerics/linear-algebra/vector.lisp
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/linear-algebra/vector.lisp Wed Aug 16 17:07:41 2006
@@ -0,0 +1,461 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; vector.lisp --
+;;;; There is no wrapper for vectors, but we do have a class (struct)
+;;;; for 'column' vectors.
+;;;; Apart from that there are some specialized generic functions on
+;;;; vectors.
+
+(in-package "CL.MATH.LINEAR-ALGEBRA")
+
+(defparameter *default-vector-element-type* 'double-float)
+
+
+(deftype column-array (&optional (type 'cl:number) (length 'cl:*))
+ `(array ,type (,length 1)))
+
+
+(defstruct (column (:constructor column (vector)))
+ "The Column Structure.
+A representation for `column' vectors."
+ (vector #() :type (or list vector)))
+
+
+(defmethod print-object ((c column) (s stream))
+ (format s "(~S ~S)" 'column (column-vector c)))
+
+
+;;;---------------------------------------------------------------------------
+;;; Protocol.
+
+(defgeneric vector-p (v)
+ (:method ((v vector)) t)
+ (:method ((v column)) t)
+ (:method ((v t)) nil)
+ )
+
+
+(defgeneric vector-element-type (v)
+ (:method ((v vector)) (array-element-type v))
+ (:method ((v column)) (array-element-type (column-vector v)))
+ )
+
+
+(defgeneric vector-length (m)
+ (:method ((x vector)) (length x))
+ (:method ((x column)) (length (column-vector x)))
+ )
+
+
+;;;---------------------------------------------------------------------------
+;;; Creation.
+
+(defun make-vector (n
+ &rest array-keys
+ &key
+ (element-type *default-vector-element-type*)
+ (column-p nil)
+ &allow-other-keys
+ )
+ (remf array-keys :column-p)
+ (let ((new-v (apply #'make-array n
+ (append array-keys
+ (list :initial-element (coerce 0 element-type)
+ :element-type
+ element-type)))))
+ (if column-p
+ (column new-v)
+ new-v)))
+
+
+(defun make-zero-vector (n
+ &rest array-keys
+ &key
+ (element-type *default-matrix-element-type*)
+ &allow-other-keys)
+ (when (getf array-keys :initial-element)
+ (warn "MAKE-ZERO-MATRIX passed and :INITIAL-ELEMENT argument.")
+ (remf array-keys :initial-element))
+ (apply #'make-vector
+ n
+ (append array-keys
+ (list :initial-element (coerce 0 element-type)
+ :element-type element-type))))
+
+
+;;;---------------------------------------------------------------------------
+;;; copy-vector
+;;; Uses the code from K Pitman and B. Margolin appeared on CLL a long
+;;; long time ago.
+;;; Note that at least on LW this code has some of the usual problems
+;;; with DOUBLE-FLOATS.
+#|
+(defun copy-array (array)
+ (adjust-array (make-array (array-dimensions array)
+ :displaced-to array
+ :element-type (array-element-type array))
+ (array-dimensions array)
+ :displaced-to nil))
+|#
+
+(defgeneric copy-vector (v)
+ (:documentation "Returns a (shallow) copy of vector V.
+The effect are the same as CL:COPY-SEQ for CL vectors while columns
+are treated accordingly."))
+
+
+(defmethod copy-vector ((v vector))
+ "COPY-VECTOR method for VECTORs."
+ (cl:copy-seq v))
+
+
+(defmethod copy-vector ((v column))
+ "COPY-VECTOR method for COLUMN vectors."
+ (column (cl:copy-seq (column-vector v))))
+
+
+(defmethod copy-vector ((v array))
+ (check-type v column-array)
+ (adjust-array (make-array (array-dimensions v)
+ :displaced-to v
+ :element-type (array-element-type v))
+ (array-dimensions v)
+ :displaced-to nil))
+
+
+
+#|
+(defmethod copy-vector ((v vector)
+ &optional
+ (result
+ (adjust-array
+ (make-array (length v)
+ :displaced-to v
+ :element-type (array-element-type v))
+ (array-dimensions v)
+ :displaced-to nil)
+ result-supplied-p))
+ "COPY-VECTOR method for VECTORs."
+ (when result-supplied-p
+ (assert (vector-p result))
+ (assert (cl:= (length result) (vector-length v)))
+ (map-into result #'identity v))
+ result)
+
+
+(defmethod copy-vector ((c column)
+ &optional
+ (result
+ (let ((v (column-vector c)))
+ (column
+ (adjust-array
+ (make-array (length v)
+ :displaced-to v
+ :element-type (array-element-type v))
+ (array-dimensions v)
+ :displaced-to nil)))
+ result-supplied-p))
+ "COPY-VECTOR method for COLUMN vectors."
+ (when result-supplied-p
+ (assert (column-p result))
+ (assert (cl:= (vector-length result) (vector-length c)))
+ (map-into (column-vector result) #'identity (column-vector c)))
+ result)
+|#
+
+
+
+;;;---------------------------------------------------------------------------
+;;; Equality.
+
+(defmethod =%2 ((x vector) (y vector))
+ (and (cl:= (length x) (length y))
+ (every #'=%2 x y)))
+
+
+(defmethod =%2 ((x vector) (y column))
+ (and (or (= 1 (length x) (vector-length y)) ; Only cases when this may happen.
+ (= 0 (length x) (vector-length y)))
+ (=%2 x (column-vector y))))
+
+
+(defmethod =%2 ((y column) (x vector))
+ (=%2 x y))
+
+
+;;;---------------------------------------------------------------------------
+;;; zerop
+
+(defmethod zerop ((x vector))
+ (cl:every #'zerop x))
+
+
+(defmethod zerop ((x column))
+ (cl:every #'zerop (column-vector x)))
+
+
+;;;---------------------------------------------------------------------------
+;;; Addition.
+
+(defmethod +%2 :before ((x vector) (y vector)
+ &optional
+ (r nil r-supplied-p))
+ (unless (eq y r)
+ (let ((xl (length x)))
+ (assert (cl:= xl (length y)))
+ (when r-supplied-p
+ (assert (vector-p r))
+ (assert (cl:= xl (vector-length r)))
+ ))))
+
+
+(defmethod +%2 :before ((x vector) (y column)
+ &optional
+ (r nil r-supplied-p))
+ (declare (ignore r r-supplied-p))
+ (assert (or (cl:= 1 (length x) (vector-length y))
+ (cl:= 0 (length x) (vector-length y))))
+ )
+
+
+(defmethod +%2 :before ((y column) (x vector)
+ &optional
+ (r nil r-supplied-p))
+ (declare (ignore r r-supplied-p))
+ (assert (or (cl:= 1 (length x) (vector-length y))
+ (cl:= 0 (length x) (vector-length y))))
+ )
+
+
+(defmethod +%2 :before ((x number) (y vector)
+ &optional
+ (r nil r-supplied-p))
+ (unless (eq y r)
+ (when r-supplied-p
+ (assert (vector-p r))
+ (assert (cl:= (length y) (vector-length r)))
+ )))
+
+
+(defmethod +%2 ((x number) (y vector) &optional (r (copy-vector y)))
+ (dotimes (i (length y) r)
+ (setf (aref r i) (+%2 x (aref y i)))))
+
+
+(defmethod +%2 ((y vector) (x number) &optional (r (copy-vector y)))
+ (+%2 x y r))
+
+
+(defmethod +%2 ((x vector) (y vector) &optional (r (copy-vector y)))
+ (dotimes (i (length x) r)
+ (setf (aref r i) (+%2 (aref x i) (aref y i)))))
+
+
+(defmethod +%2 ((x vector) (y column)
+ &optional (r (copy-vector (column-vector y))))
+ ;; This is valid only for vectors and columns of length 1 (or 0).
+ (+%2 x (column-vector y) r))
+
+
+(defmethod +%2 ((x column) (y vector) &optional (r (copy-vector y)))
+ ;; This is valid only for vectors and columns of length 1 (or 0).
+ (+%2 (column-vector x) y r))
+
+
+;;;---------------------------------------------------------------------------
+;;; Subtraction.
+
+(defmethod -%2 :before ((x vector) (y vector)
+ &optional
+ (r nil r-supplied-p))
+ (let ((xl (length x)))
+ (assert (cl:= xl (length y)))
+ (when r-supplied-p
+ (assert (vector-p r))
+ (assert (cl:= xl (vector-length r)))
+ )))
+
+
+(defmethod -%2 :before ((x vector) (y column)
+ &optional
+ (r nil r-supplied-p))
+ (declare (ignore r r-supplied-p))
+ (assert (or (cl:= 1 (length x) (vector-length y))
+ (cl:= 0 (length x) (vector-length y))))
+ )
+
+
+(defmethod -%2 :before ((y column) (x vector)
+ &optional
+ (r nil r-supplied-p))
+ (declare (ignore r r-supplied-p))
+ (assert (or (cl:= 1 (length x) (vector-length y))
+ (cl:= 0 (length x) (vector-length y))))
+ )
+
+
+(defmethod -%2 :before ((x number) (y vector)
+ &optional
+ (r nil r-supplied-p))
+ (unless (eq y r)
+ (when r-supplied-p
+ (assert (vector-p r))
+ (assert (cl:= (length y) (vector-length r)))
+ )))
+
+
+(defmethod -%2 ((x vector) (y vector) &optional (r (copy-vector y)))
+ (dotimes (i (length x) r)
+ (setf (aref r i) (-%2 (aref x i) (aref y i)))))
+
+
+
+(defmethod -%2 ((x vector) (y column)
+ &optional (r (copy-vector (column-vector y))))
+ ;; This is valid only for vectors and columns of length 1 (or 0).
+ (-%2 x (column-vector y) r))
+
+
+(defmethod -%2 ((x column) (y vector) &optional (r (copy-vector y)))
+ ;; This is valid only for vectors and columns of length 1 (or 0).
+ (-%2 (column-vector x) r r))
+
+
+;;;---------------------------------------------------------------------------
+;;; Multiplication.
+;;; Note that left and right multiplications are different!
+;;; They follow the matrix multiplication rules.
+;;; Here we treat only the vector x column multiplication as the
+;;; column x vector has a matrix as a result.
+
+(defmethod *%2 :before ((x vector) (y column) &optional r)
+ (declare (ignore r))
+ (assert (cl:= (cl:length x) (vector-length y))))
+
+
+(defmethod *%2 ((x vector) (y column) &optional r)
+ (declare (ignore r))
+ (let ((result (coerce 0 (array-element-type x)))
+ (y (column-vector y))
+ )
+ (dotimes (i (cl:length x) result)
+ (setf result (+%2 result (*%2 (aref x i) (aref y i)))))
+ ))
+
+;;; The next one breaks the correct semantics passing of arguments,
+;;; but it is provided as a convenience.
+
+(defmethod *%2 ((x vector) (y vector) &optional r)
+ (declare (ignore r))
+ (assert (cl:= (cl:length x) (cl:length y)))
+ (let ((result (coerce 0 (array-element-type x))))
+ (dotimes (i (cl:length x) result)
+ (setf result (+%2 result (*%2 (aref x i) (aref y i)))))))
+
+
+;;; Scalar multiplication.
+
+(defmethod *%2 :before ((x number) (v vector) &optional (r nil r-supplied-p))
+ (when r-supplied-p
+ (assert (cl:vectorp r))
+ (assert (cl:= (cl:length v) (vector-length r)))))
+
+
+(defmethod *%2 :before ((x number) (v column) &optional (r nil r-supplied-p))
+ (when r-supplied-p
+ (assert (column-p r))
+ (assert (cl:= (vector-length v) (vector-length r)))))
+
+
+(defmethod *%2 ((x number) (v vector) &optional (r (copy-vector v)))
+ (map-into r (lambda (e) (*%2 x e)) v))
+
+
+(defmethod *%2 ((v vector) (x number) &optional (r (copy-vector v)))
+ (map-into r (lambda (e) (*%2 x e)) v))
+
+
+(defmethod *%2 ((x number) (v column) &optional (r (copy-vector v)))
+ (map-into (column-vector r) (lambda (e) (*%2 x e)) v)
+ r)
+
+
+(defmethod *%2 ((v column) (x number) &optional (r (copy-vector v)))
+ (*%2 x v r))
+
+
+;;;---------------------------------------------------------------------------
+;;; Division.
+;;; Only the simple form of division is implemented here.
+;;; The general form requires the inverse, which requires extra machinery
+
+#|
+(defmethod /%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (*%2 (/ x) y r))
+
+
+(defmethod /%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p))
+ (declare (ignore r-supplied-p))
+ (*%2 (/ x) y r))
+|#
+
+;;;---------------------------------------------------------------------------
+;;; Exponentiation.
+;;; All errors here ftb.
+
+
+;;;---------------------------------------------------------------------------
+;;; Transpose
+
+(defmethod transpose ((v vector) &optional (r (column (copy-vector v)) r-supplied-p))
+ (if r-supplied-p
+ (map-into (column-vector r) #'identity v))
+ r)
+
+
+(defmethod transpose ((v column)
+ &optional
+ (r (copy-vector (column-vector v)) r-supplied-p))
+ (if r-supplied-p
+ (map-into r #'identity (column-vector v))
+ r))
+
+
+(defmethod transpose :before ((v vector) &optional (r nil r-supplied-p))
+ (when r-supplied-p
+ (assert (column-p r))
+ (assert (cl:= (cl:length v) (vector-length r)))))
+
+
+(defmethod transpose :before ((v column) &optional (r nil r-supplied-p))
+ (when r-supplied-p
+ (assert (cl:vectorp r))
+ (assert (cl:= (cl:length r) (vector-length v)))))
+
+
+;;;---------------------------------------------------------------------------
+;;; Other functions.
+
+(defun iota (n &optional (start 0) (step 1 step-supplied-p)
+ &aux (initial-contents ()))
+ (cond ((< n start)
+ (when step-supplied-p (assert (minusp step)))
+ (loop for x from start above n by step
+ collect x into xs
+ finally (setf initial-contents xs)))
+ ((> n start)
+ (when step-supplied-p (assert (plusp step)))
+ (loop for x from start below n by step
+ collect x into xs
+ finally (setf initial-contents xs)))
+ (t ()))
+ (make-array n :initial-contents initial-contents))
+
+
+
+;;;---------------------------------------------------------------------------
+;;; Utilities.
+
+
+;;;; end of file -- matrix.lisp --
Added: trunk/common-math/numerics/numerics.system
==============================================================================
--- (empty file)
+++ trunk/common-math/numerics/numerics.system Wed Aug 16 17:07:41 2006
@@ -0,0 +1,18 @@
+;;;; -*- Mode: Lisp -*-
+
+;;;; numerics.system --
+
+(mk:defsystem "numerics"
+ :components ((:system "linear-algebra")
+ )
+ )
+
+
+(eval-when (:load-toplevel :execute)
+ (mk:add-registry-location (make-pathname :name nil
+ :type nil
+ :defaults *load-truename*))
+ )
+
+;;;; end of file -- numerics.system --
+
