[common-math-cvs] r1 - in trunk: . common-math common-math/numerics common-math/numerics/linear-algebra

16 Aug 2006

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 --
+

    

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

[common-math-cvs] r1 - in trunk: . common-math common-math/numerics common-math/numerics/linear-algebra