matlisp-commit

[Akshay S. <ak...@us...>]

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".

The branch, tensor has been updated
       via  71aca48b041b5be2cd4c6ab8d514b260bdc02b19 (commit)
       via  174d27300595c21a466a330fa34ab66fa7131bdf (commit)
       via  f25a68740987eeac4539d48b7a58d189da5b28e7 (commit)
       via  536e528120dcf8631dbb9a8d4efd9af5541e55e0 (commit)
      from  d5f7ad309ca59d41c6e405c512f9a3544be01ea2 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 71aca48b041b5be2cd4c6ab8d514b260bdc02b19
Author: Akshay Srinivasan <aks...@gm...>
Date:   Thu Jul 5 17:16:14 2012 +0530

    o Loopy has a new macro to recurse through lists of lists of lists ...
    o Added a tensor-maker macro.
    o dot now works.

diff --git a/README.org b/README.org
index 071d4fc..05bc347 100644
--- a/README.org
+++ b/README.org
@@ -26,6 +26,8 @@ This is the development branch of Matlisp.
    * ODEPACK: Add abstraction for DLSODE, and DLSODAR.
    * Tensor contraction: Hard to do very quickly.
      Might have to copy stuff into a contiguous array; like Femlisp.
+*** Syntactic sugar
+   * Add array slicing macros
 
 *** Python-bridge
     (C)Python has far too many things, that we cannot even begin to hope to replicate.
diff --git a/matlisp.asd b/matlisp.asd
index 5db4d00..ec882aa 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -80,12 +80,10 @@
 	    :components ((:file "blas")
 			 (:file "lapack")
 			 (:file "dfftpack")))
-   (:module "matlisp-essentials"
+   (:module "matlisp-base"
 	    :pathname "src/"
-	    :depends-on ("foreign-interface"
-			 "foreign-functions")
+	    :depends-on ("foreign-functions")
 	    :components ((:file "conditions")
-			 ;;
 			 (:file "standard-tensor"
 				:depends-on ("conditions"))
 			 ;;
@@ -94,28 +92,26 @@
 			 (:file "permutation"
 				:depends-on ("standard-tensor"))
 			 (:file "blas-helpers"
-				:depends-on ("standard-tensor" "permutation"))			 
-			 ;;
-			 (:file "real-tensor"
-				:depends-on ("standard-tensor"))
-			 (:file "complex-tensor"
-				:depends-on ("standard-tensor"))
-			 (:file "standard-matrix"
-				:depends-on ("standard-tensor" "real-tensor" "complex-tensor"))
-			 ;; (:file "real-matrix"
-			 ;;	    :depends-on ("standard-matrix"))
-			 ;; (:file "complex-matrix"
-			 ;;	    :depends-on ("standard-matrix"))
+				:depends-on ("standard-tensor" "permutation"))
 			 (:file "print"
-				:depends-on ("standard-tensor" "standard-matrix"))
-			 ;;Copy, Scal
-			 (:file "copy"
-				:depends-on ("real-tensor" "complex-tensor" "loopy"))
+				:depends-on ("standard-tensor"))))
+   (:module "matlisp-classes"
+	    :pathname "src/"
+	    :depends-on ("matlisp-base")
+	    :components ((:file "real-tensor")
+			 (:file "complex-tensor")
+			 (:file "matrix"
+				:depends-on ("real-tensor" "complex-tensor"))))
+   (:module "matlisp-level-1"
+	    :pathname "src/"
+	    :depends-on ("matlisp-base" "matlisp-classes" "foreign-functions")
+	    :components ((:file "tensor-maker")
+			 (:file "copy")
+			 (:file "dot")
 			 (:file "scal"
-				:depends-on ("copy" "loopy"))
+				:depends-on ("copy"))
 			 (:file "realimag"
-				:depends-on ("real-tensor" "complex-tensor" "copy"))
-			 ))))
+				:depends-on ("copy"))))))
 
 
 ;; (defclass f2cl-cl-source-file (asdf:cl-source-file)
diff --git a/packages.lisp b/packages.lisp
index 5226d68..d202ce8 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -162,7 +162,7 @@
 	   #:recursive-append #:unquote-args #:flatten
 	   #:format-to-string #:string+
 	   #:linear-array-type
-	   #:seq-max #:seq-min
+	   #:list-dimensions
 	   ;;Macros
 	   #:when-let #:if-let #:if-ret #:with-gensyms #:let-rec
 	   #:mlet* #:make-array-allocator
diff --git a/src/complex-tensor.lisp b/src/complex-tensor.lisp
index 4a636cf..aecac89 100644
--- a/src/complex-tensor.lisp
+++ b/src/complex-tensor.lisp
@@ -1,24 +1,24 @@
 (in-package :matlisp)
 
-(eval-when (load eval compile)
-  (deftype complex-base-type ()
-    "The type of the elements stored in a COMPLEX-MATRIX"
-    'double-float)
-  
-  (deftype complex-base-array (size)
-    "The type of the storage structure for a COMPLEX-MATRIX"
-    `(simple-array complex-base-type (,size)))
- 
-  (deftype complex-type ()
-    "Complex number with Re, Im parts in complex-base-type."
-    '(cl:complex complex-base-type))
-  )
+(deftype complex-base-type ()
+  "The type of the elements stored in a COMPLEX-MATRIX"
+  'double-float)
+
+(deftype complex-base-array (size)
+  "The type of the storage structure for a COMPLEX-MATRIX"
+  `(simple-array complex-base-type (,size)))
+
+(deftype complex-type ()
+  "Complex number with Re, Im parts in complex-base-type."
+  '(cl:complex complex-base-type))
 ;;
 
 (definline allocate-complex-store (size)
-"(allocate-complex-store size)
-Allocates real storage of size (* SIZE 2).
-Default initial-element = 0d0."
+  "
+  (allocate-complex-store size)
+  Allocates real storage of size (* SIZE 2).
+  Default initial-element = 0d0.
+"
   (make-array (* 2 size) :element-type 'complex-base-type
 	      :initial-element (coerce 0 'complex-base-type)))
 
@@ -87,11 +87,3 @@ Cannot hold complex numbers."))
 		       "~11,5,,,,,'Eg"
 		       "#C(~11,4,,,,,'Ee ~11,4,,,,,'Ee)")
 	    realpart imagpart)))
-
-;;
-(defun make-complex-tensor-dims (&rest subs)
-  (let* ((dims (make-index-store subs))
-	 (ss (reduce #'* dims))
-	 (store (allocate-complex-store ss)))
-    (make-instance 'complex-tensor :store store :dimensions dims)))
-
diff --git a/src/conditions.lisp b/src/conditions.lisp
index 27986e0..842624e 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -33,7 +33,7 @@
 ;;Generic conditions---------------------------------------------;;
 (defcondition generic-error (error)
   ((message :reader message :initarg :message :initform ""))
-  (:method print-method ((c generic-error) stream)
+  (:method print-object ((c generic-error) stream)
 	   (format stream (message c))))
 
 (defcondition invalid-type (generic-error)
@@ -64,9 +64,25 @@
    (to :reader to :initarg :to))
   (:documentation "Cannot coerce one type into another.")
   (:method print-object ((c coercion-error) stream)
-	   (format stream "Cannot coerce ~a into ~a." (from c) (to c))
+	   (format stream "Cannot coerce ~a into ~a.~%" (from c) (to c))
 	   (call-next-method)))
 
+(defcondition out-of-bounds-error (generic-error)
+  ((requested :reader requested :initarg :requested)
+   (bound :reader bound :initarg :bound))
+  (:documentation "General out-of-bounds error")
+  (:method print-object ((c out-of-bounds-error) stream)
+	   (format stream "Out-of-bounds error, requested index : ~a, bound : ~a.~%" (requested c) (bound c))
+	   (call-next-method)))
+
+(defcondition non-uniform-bounds-error (generic-error)
+  ((assumed :reader assumed :initarg :assumed)
+   (found :reader found :initarg :found))
+  (:documentation "Bounds are not uniform")
+  (:method print-object ((c out-of-bounds-error) stream)
+	   (format stream "The bounds are not uniform, assumed bound : ~a, now found to be : ~a.~%" (assumed c) (found c))
+	   (call-next-method)))
+	   
 ;;Tensor conditions----------------------------------------------;;
 (define-condition tensor-error (error)
   ;;Optional argument for error-handling.
@@ -92,6 +108,12 @@
   (:report (lambda (c stream)
 	     (format stream "Given tensor with rank ~A, is not a matrix." (rank c)))))
 
+(define-condition tensor-not-vector (tensor-error)
+  ((tensor-rank :reader rank :initarg :rank))
+  (:documentation "Given tensor is not a vector.")
+  (:report (lambda (c stream)
+	     (format stream "Given tensor with rank ~A, is not a vector." (rank c)))))
+
 (define-condition tensor-index-out-of-bounds (tensor-error)
   ((argument :reader argument :initarg :argument)
    (index :reader index :initarg :index)
diff --git a/src/dot.lisp b/src/dot.lisp
index 8983caa..1affebc 100644
--- a/src/dot.lisp
+++ b/src/dot.lisp
@@ -97,149 +97,67 @@
   otherwise.
 "))
 
-(defmethod dot ((x number) (y number) &optional (conjugate-p nil))
+(defmethod dot ((x number) (y number) &optional (conjugate-p t))
   (if conjugate-p
       (* (conjugate x) y)
       (* x y)))
 
-(defmethod dot :before ((x standard-matrix) (y standard-matrix) &optional conjugate-p)
+(defmethod dot :before ((x standard-tensor) (y standard-tensor) &optional (conjugate-p t))
   (declare (ignore conjugate-p))
-  (if (not (row-or-col-vector-p x))
-      (error "argument X to DOT is not a row or column vector")
-    (if (not (row-or-col-vector-p y))
-	(error "argument Y to DOT is not a row or column vector")
-      (let ((nxm-x (number-of-elements x))
-	    (nxm-y (number-of-elements y)))
-	(declare (type fixnum nxm-x nxm-y))
-	(if (not (= nxm-x nxm-y))
-	    (error "arguments X,Y to DOT are not of the same size"))))))
-
-(defmethod dot ((x real-matrix) (y real-matrix) &optional conjugate-p)
+  (unless (and (vector-p x) (vector-p y))
+    (error 'tensor-not-vector
+	   :rank (cond
+		   ((not (vector-p x))
+		    (rank x))
+		   ((not (vector-p y))
+		    (rank y)))))
+  (unless (idx= (dimensions x) (dimensions y))
+    (error 'tensor-dimension-mismatch)))
+
+(defmethod dot ((x real-tensor) (y real-tensor) &optional (conjugate-p t))
   (declare (ignore conjugate-p))
-  (let ((nxm (number-of-elements x)))
-    (declare (type fixnum nxm))
-    (ddot nxm (store x) 1 (store y) 1)))
+  (ddot (number-of-elements x)
+	(store x) (aref (strides x) 0)
+	(store y) (aref (strides y) 0)
+	(head x) (head y)))
 
-;;#+(or :cmu :sbcl)
-(defmethod dot ((x real-matrix) (y complex-matrix) &optional conjugate-p)
+(defmethod dot ((x real-tensor) (y complex-tensor) &optional (conjugate-p t))
   (declare (ignore conjugate-p))
-  (let ((nxm (number-of-elements x))
-	(store-x (store x))
-	(store-y (store y)))
-    (declare (type fixnum nxm)
-	     (type (real-matrix-store-type *) store-x)
-	     (type (complex-matrix-store-type *) store-y))
-
-    (let ((realpart (ddot nxm store-x 1 store-y 2))
-	  (imagpart (with-vector-data-addresses ((addr-x store-x)
-						 (addr-y store-y))
-			  (incf-sap addr-y :double-float)
-			  (ddot nxm addr-x 1 addr-y 2))))
-
-      (declare (type complex-matrix-element-type realpart imagpart))
-
-      #+:complex-arg-implies-complex-result
-      (complex realpart imagpart)
-      #-:complex-arg-implies-comples-result
-      (if (zerop imagpart)
-	  realpart
-	(complex realpart imagpart))
-      )))
-
-
-;;#+:allegro
-#+nil
-(defmethod dot ((x real-matrix) (y complex-matrix) &optional conjugate-p)
-  (declare (ignore conjugate-p))
-  (let ((nxm (number-of-elements x)))
-    (declare (type fixnum nxm))
-
-    (let ((realpart 0.0d0)
-	  (imagpart 0.0d0))
-      (declare (type complex-matrix-element-type realpart imagpart))
-
-      (dotimes (i nxm)
-        (declare (type fixnum i))
-        (let ((x-elt (matrix-ref x i))
-	      (y-elt (matrix-ref y i)))
-	  (incf realpart (+ x-elt (realpart y-elt)))
-	  (incf imagpart (+ x-elt (imagpart y-elt)))))
-
-
-      #+:complex-arg-implies-complex-result
-      (complex realpart imagpart)
-      #-:complex-arg-implies-comples-result
-      (if (zerop imagpart)
-	  realpart
-	(complex realpart imagpart))
-      )))
-
-;;#+(or :cmu :sbcl)
-(defmethod dot ((x complex-matrix) (y real-matrix) &optional (conjugate-p t))
-  (let ((nxm (number-of-elements x))
-	(store-x (store x))
-	(store-y (store y)))
-    (declare (type fixnum nxm)
-	     (type (real-matrix-store-type *) store-y)
-	     (type (complex-matrix-store-type *) store-x))
-
-    (let ((realpart (ddot nxm store-x 2 store-y 1))
-	  (imagpart (with-vector-data-addresses ((addr-x store-x)
-						  (addr-y store-y))
-			  (incf-sap addr-x :double-float)
-			  (ddot nxm addr-x 2 addr-y 1))))
-
-      (declare (type complex-matrix-element-type realpart imagpart))
-
-      (if conjugate-p
-	  (setq imagpart (- imagpart)))
-
-      #+:complex-arg-implies-complex-result
-      (complex realpart  imagpart)
-      #-:complex-arg-implies-comples-result
-      (if (zerop imagpart)
-	  realpart
-	(complex realpart imagpart))
-      )))
-
-#+nil
-(defmethod dot ((x complex-matrix) (y real-matrix) &optional (conjugate-p t))
-  (let ((nxm (number-of-elements x)))
-    (declare (type fixnum nxm))
-
-    (let ((realpart 0.0d0)
-	  (imagpart 0.0d0))
-      (declare (type complex-matrix-element-type realpart imagpart))
-
-      (dotimes (i nxm)
-        (declare (type fixnum i))
-        (let ((x-elt (matrix-ref x i))
-	      (y-elt (matrix-ref y i)))
-	  (incf realpart (+ y-elt (realpart x-elt)))
-	  (incf imagpart (+ y-elt (imagpart x-elt)))))
-
-      (if conjugate-p
-	  (setq imagpart (- imagpart)))
-
-      #+:complex-arg-implies-complex-result
-      (complex realpart  imagpart)
-      #-:complex-arg-implies-comples-result
-      (if (zerop imagpart)
-	  realpart
-	(complex realpart imagpart))
-      )))
-
-(defmethod dot ((x complex-matrix) (y complex-matrix) &optional (conjugate-p t))
-  (let* ((nxm (number-of-elements x))
-	 (store-x (store x))
-	 (store-y (store y))
-	 (dot (if conjugate-p
-		  (zdotc nxm store-x 1 store-y 1)
-		  (zdotu nxm store-x 1 store-y 1))))
-    
-    #-:complex-arg-implies-complex-result
-    (if (zerop (imagpart dot))
-	(realpart dot)
-	dot)
-    #+:complex-arg-implies-complex-result
-    dot))
+  (let ((nele (number-of-elements x))
+	(std-x (aref (strides x) 0))
+	(hd-x (head x))
+	(std-y (aref (strides y) 0))
+	(hd-y (head y)))
+    (declare (type index-type nele std-x std-y hd-x hd-y))
+    (let ((rpart (ddot nele (store x) std-x (store y) (* 2 std-y) hd-x (* 2 hd-y)))
+	  (ipart (ddot nele (store x) std-x (store y) (* 2 std-y) hd-x (1+ (* 2 hd-y)))))
+      (declare (type complex-base-type rpart ipart))
+      (if (zerop ipart)
+	  rpart
+	  (complex rpart ipart)))))
+
+(defmethod dot ((x complex-tensor) (y real-tensor) &optional (conjugate-p t))
+  (let ((cres (dot y x)))
+    (if conjugate-p
+	(conjugate cres)
+	cres)))
+
+(defmethod dot ((x complex-tensor) (y complex-tensor) &optional (conjugate-p t))
+  (let ((nele (number-of-elements x))
+	(std-x (aref (strides x) 0))
+	(hd-x (head x))
+	(std-y (aref (strides y) 0))
+	(hd-y (head y)))
+    (declare (type index-type nele std-x hd-x std-y hd-y))
+    (let ((ret (if conjugate-p
+		   (zdotc nele
+			  (store x) std-x
+			  (store y) std-y
+			  hd-x hd-y)
+		   (zdotu nele
+			  (store x) std-x
+			  (store y) std-y
+			  hd-x hd-y))))
+      (if (zerop (imagpart ret))
+	  (realpart ret)
+	  ret))))
diff --git a/src/loopy.lisp b/src/loopy.lisp
index 32037bc..1828f66 100644
--- a/src/loopy.lisp
+++ b/src/loopy.lisp
@@ -8,38 +8,28 @@
     {with (loop-order {:row-major :col-major})}
     {with (linear-sums
 	      {(offsets {stride-seq})}*)}
-    {with (variables
-	      {(vars init &key type)}*)}
     {do ({code}*)})
 
   Examples:
-  > (mod-dotimes (idx (vidx 2 2))
-      with (linear-sums (of (vidx 2 1)))
+  > (mod-dotimes (idx (idxv 2 2))
+      with (linear-sums (of (idxv 2 1)))
       do (format t \"~a ~a~%\" idx of))
   #(0 0) 0
   #(0 1) 1
   #(1 0) 2
   #(1 1) 3
 
-  > (mod-dotimes (idx (vidx 2 2))
+  > (mod-dotimes (idx (idxv 2 2))
       with (loop-order :col-major)
-      with (linear-sums (of (vidx 2 1)))
+      with (linear-sums (of (idxv 2 1)))
       do (format t \"~a ~a~%\" idx of))
   #(0 0) 0
   #(1 0) 2
   #(0 1) 1
   #(1 1) 3
 
-  > (mod-dotimes (idx (vidx 2 2))
-      with (variables (tmp 1d0 :type double-float))
-      with (linear-sums (of (vidx 2 1)))
-      do (progn
-	   (format t \"~a ~a ~a~%\" idx of tmp)
-	   (incf tmp)))
-  #(0 0) 0 0d0
-  #(0 1) 1 1d0
-  #(1 0) 2 2d0
-  #(1 1) 3 3d0
+  Make sure that \"do\" is specified at the end. Parser stops
+  at the first 'do it finds.
 "
   (check-type idx symbol)
   (labels ((parse-code (body ret)
@@ -50,6 +40,11 @@
 		(multiple-value-bind (indic decl) (parse-with (cadr body))
 		  (setf (getf ret indic) decl))
 		(parse-code (cddr body) ret))
+	       ;;Let's not do too much
+	       #+nil
+	       ((eq (car body) 'finally)
+		(setf (getf ret :finally) (second body))
+		(parse-code (cddr body) ret))
 	       ((eq (car body) 'do)
 		(values (cadr body) ret))
 	       (t (error 'unknown-token :token (car body) :message "Error in macro: mod-dotimes -> parse-code.~%"))))
@@ -66,6 +61,8 @@
 	       ((and (eq (car code) 'loop-order)
 		     (member (cadr code) '(:row-major :col-major)))
 		(values :loop-order (second code)))
+	       ;;Useless without a finally clause
+	       #+nil
 	       ((eq (car code) 'variables)
 		(values :variables
 			(loop for decl in (cdr code)
@@ -119,4 +116,111 @@
 						  (unless (= ,cstrd 0)
 						    (incf ,(getf decl :offset-sym) ,cstrd)))))
 				(return t)))
-			  finally (return nil))))))))))))
+			  finally (return nil))))
+	      ,@(unless (null (getf sdecl :finally))
+			`(finally (,@(getf sdecl :finally))))))))))))
+
+(defmacro list-loop ((idx ele lst) &rest body)
+  "
+  (list-loop (idx ele {list}) compound-form*)
+
+  Examples:
+  > (list-loop (idx ele '((1 2) (4 5)))
+      with (linear-sums (of (idxv 2 1)))
+      do (format t \"~a ~a ~a~%\" idx of ele))
+  #(0 0) 0 1
+  #(0 1) 1 2
+  #(1 0) 2 4
+  #(1 1) 3 5
+"
+  (check-type idx symbol)
+  (check-type ele symbol)
+  (labels ((parse-code (body ret)
+	     (cond
+	       ((null body)
+		(values nil ret))
+	       ((eq (car body) 'with)
+		(multiple-value-bind (indic decl) (parse-with (cadr body))
+		  (setf (getf ret indic) decl))
+		(parse-code (cddr body) ret))
+	       ;;Let's not do too much.
+	       #+nil
+	       ((eq (car body) 'finally)
+		(setf (getf ret :finally) (second body))
+		(parse-code (cddr body) ret))
+	       ((eq (car body) 'do)
+		(values (cadr body) ret))
+	       (t (error 'unknown-token :token (car body) :message "Error in macro: mod-dotimes -> parse-code.~%"))))
+	   (parse-with (code)
+	     (cond
+	       ((eq (car code) 'linear-sums)
+		(values :linear-sums
+			(loop for decl in (cdr code)
+			   collect (destructuring-bind (offst strds &optional (init 0)) decl
+				     (list :offset-sym offst
+					   :offset-init init
+					   :stride-sym (gensym (string+ (symbol-name offst) "-stride"))
+					   :stride-expr strds)))))
+	       ;;Traversing the list the other way is far too inefficient and/or too hard to do.
+	       #+nil
+	       ((and (eq (car code) 'loop-order)
+		     (member (cadr code) '(:row-major :col-major)))
+		(values :loop-order (second code)))
+	       ;;Useless without a finally clause.
+	       #+nil
+	       ((eq (car code) 'variables)
+		(values :variables
+			(loop for decl in (cdr code)
+			   collect (destructuring-bind (sym init &key type) decl
+				     (list :variable sym
+					   :init init
+					   :type type)))))
+	       (t (error 'unknown-token :token (car code) :message "Error in macro: mod-dotimes -> parse-with.~%")))))
+    (multiple-value-bind (code sdecl) (parse-code body nil)
+      (with-gensyms (lst-sym dims-sym rank-sym lst-rec-sym lst-rec-count-sym lst-rec-lst-sym)
+	`(let ((,lst-sym ,lst))
+	   (declare (type list ,lst-sym))
+	   (let ((,dims-sym (make-index-store (list-dimensions ,lst-sym))))
+	     (declare (type (index-array *) ,dims-sym))
+	     (let ((,rank-sym (array-dimension ,dims-sym 0)))
+	       (declare (type index-type ,rank-sym))
+	       (let ((,idx (allocate-index-store ,rank-sym))
+		     ,@(mapcar #'(lambda (x) `(,(getf x :offset-sym) ,(getf x :offset-init))) (getf sdecl :linear-sums))
+		     ,@(mapcar #'(lambda (x) `(,(getf x :stride-sym) ,(getf x :stride-expr))) (getf sdecl :linear-sums))
+		     ,@(mapcar #'(lambda (x) `(,(getf x :variable) ,(getf x :init))) (getf sdecl :variables)))
+		 (declare (type (index-array *) ,idx)
+			  ,@(when (getf sdecl :linear-sums)
+				  `((type (index-array *) ,@(mapcar #'(lambda (x) (getf x :stride-sym)) (getf sdecl :linear-sums)))
+				    (type index-type ,@(mapcar #'(lambda (x) (getf x :offset-sym)) (getf sdecl :linear-sums)))))
+			  ,@(loop for x in (getf sdecl :variables)
+			       unless (null (getf x :type))
+			       collect `(type ,(getf x :type) ,(getf x :variable))))
+		 (labels ((,lst-rec-sym (,lst-rec-count-sym ,lst-rec-lst-sym)
+			    (if (null ,lst-rec-lst-sym)
+				(progn
+				  (unless (= (aref ,idx ,lst-rec-count-sym) (aref ,dims-sym ,lst-rec-count-sym))
+				    (error 'non-uniform-bounds-error :assumed (aref ,dims-sym ,lst-rec-count-sym) :found ,lst-rec-count-sym
+					   :message "Error in list-loop, given list is not uniform in dimensions."))
+				  (setf (aref ,idx ,lst-rec-count-sym) 0)
+				  ,@(loop
+				       for decl in (getf sdecl :linear-sums)
+				       collect `(decf ,(getf decl :offset-sym) (* (aref ,(getf decl :stride-sym) ,lst-rec-count-sym) (aref ,dims-sym ,lst-rec-count-sym))))
+				  ,@(if (null (getf sdecl :finally))`(nil)
+					`((when (= ,lst-rec-count-sym 0)
+					    ,(getf sdecl :finally)))))
+				(progn
+				  ;;list-dimensions does not parse the entire list, just goes through caaa..r's to find out the
+				  ;;dimensions if it is uniform.
+				  (unless (< -1 (aref ,idx ,lst-rec-count-sym) (aref ,dims-sym ,lst-rec-count-sym))
+				    (error 'out-of-bounds-error :requested ,lst-rec-count-sym :bound (aref ,dims-sym ,lst-rec-count-sym)
+					   :message "Error in list-loop, given list is not uniform in dimensions."))
+				  (if (consp (car ,lst-rec-lst-sym))
+				      (,lst-rec-sym (1+ ,lst-rec-count-sym) (car ,lst-rec-lst-sym))
+				      (let ((,ele (car ,lst-rec-lst-sym)))
+					,code))
+				  (incf (aref ,idx ,lst-rec-count-sym))
+				  ,@(loop
+				       for decl in (getf sdecl :linear-sums)
+				       collect `(incf ,(getf decl :offset-sym) (aref ,(getf decl :stride-sym) ,lst-rec-count-sym)))
+				  (,lst-rec-sym ,lst-rec-count-sym (cdr ,lst-rec-lst-sym))))))
+		   (,lst-rec-sym 0 ,lst-sym))))))))))
diff --git a/src/matrix.lisp b/src/matrix.lisp
index 4454c6f..e8fa449 100644
--- a/src/matrix.lisp
+++ b/src/matrix.lisp
@@ -1,511 +1,87 @@
-;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Copyright (c) 2000 The Regents of the University of California.
-;;; All rights reserved. 
-;;; 
-;;; Permission is hereby granted, without written agreement and without
-;;; license or royalty fees, to use, copy, modify, and distribute this
-;;; software and its documentation for any purpose, provided that the
-;;; above copyright notice and the following two paragraphs appear in all
-;;; copies of this software.
-;;; 
-;;; IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
-;;; FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
-;;; ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
-;;; THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
-;;; SUCH DAMAGE.
-;;;
-;;; THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
-;;; INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-;;; MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
-;;; PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
-;;; CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
-;;; ENHANCEMENTS, OR MODIFICATIONS.
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; Originally written by Raymond Toy
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;;
-;;; $Id: matrix.lisp,v 1.16 2011/01/25 18:36:56 rtoy Exp $
-;;;
-;;; $Log: matrix.lisp,v $
-;;; Revision 1.16  2011/01/25 18:36:56  rtoy
-;;; Merge changes from automake-snapshot-2011-01-25-1327 to get the new
-;;; automake build infrastructure.
-;;;
-;;; Revision 1.15.2.1  2011/01/25 18:16:53  rtoy
-;;; Use cl:real instead of real.
-;;;
-;;; Revision 1.15  2010/12/12 02:07:31  rtoy
-;;; matrix.lisp:
-;;;
-;;; o Apply patch from Nicolas Neuss for matrices with 0 dimensions.  (See
-;;;   <http://sourceforge.net/mailarchive/message.php?msg_id=1124576>)
-;;;
-;;; print.lisp:
-;;; o Apparently the above patch to print was also applied previously.  We
-;;;   just fix a bug in printing 0xm matrices.  Just exit early if the
-;;;   matrix has no dimensions.
-;;;
-;;; Revision 1.14  2004/05/24 16:34:22  rtoy
-;;; More SBCL support from Robert Sedgewick.  The previous SBCL support
-;;; was incomplete.
-;;;
-;;; Revision 1.13  2003/05/31 22:20:26  rtoy
-;;; o Add some support for CMUCL with Gerd's PCL so we can inline
-;;;   accessors and such for the matrix classes.
-;;; o Only use one, system-independent, standard-matrix class.
-;;; o Try to declare the types of the slots of the matrix classes
-;;; o FORTRAN-MATRIX-INDEXING changed to use fixnum arithmetic.
-;;;
-;;; Revision 1.12  2003/03/09 14:26:30  rtoy
-;;; Forgot one more :type 'fixnum bug.  From Gerd Moellmann.
-;;;
-;;; Revision 1.11  2003/02/19 21:59:52  rtoy
-;;; Correct the slot type declarations.
-;;;
-;;; Revision 1.10  2001/10/29 18:00:28  rtoy
-;;; Updates from M. Koerber to support QR routines with column pivoting:
-;;;
-;;; o Add an integer4 type and allocate-integer4-store routine.
-;;; o Add the necessary Fortran routines
-;;; o Add Lisp interface to the Fortran routines
-;;; o Update geqr for the new routines.
-;;;
-;;; Revision 1.9  2001/06/22 12:51:49  rtoy
-;;; o Added ALLOCATE-REAL-STORE and ALLOCATE-COMPLEX-STORE functions to
-;;;   allocate appropriately sized arrays for holding real and complex
-;;;   matrix elements.
-;;; o Use it to allocate space.
-;;;
-;;; Revision 1.8  2000/10/04 15:56:50  simsek
-;;; o Fixed bug in (MAKE-COMPLEX-MATRIX n)
-;;;
-;;; Revision 1.7  2000/07/11 18:02:03  simsek
-;;; o Added credits
-;;;
-;;; Revision 1.6  2000/07/11 02:11:56  simsek
-;;; o Added support for Allegro CL
-;;;
-;;; Revision 1.5  2000/05/11 18:28:10  rtoy
-;;; After the great standard-matrix renaming, row-vector-p and
-;;; col-vector-p were swapped.
-;;;
-;;; Revision 1.4  2000/05/11 18:02:55  rtoy
-;;; o After the great standard-matrix renaming, I forgot a few initargs
-;;;   that needed to be changed
-;;; o MAKE-REAL-MATRIX and MAKE-COMPLEX-MATRIX didn't properly handle
-;;;   things like #(1 2 3 4) and #((1 2 3 4)).  Make them accept these.
-;;;
-;;; Revision 1.3  2000/05/08 17:19:18  rtoy
-;;; Changes to the STANDARD-MATRIX class:
-;;; o The slots N, M, and NXM have changed names.
-;;; o The accessors of these slots have changed:
-;;;      NROWS, NCOLS, NUMBER-OF-ELEMENTS
-;;;   The old names aren't available anymore.
-;;; o The initargs of these slots have changed:
-;;;      :nrows, :ncols, :nels
-;;;
-;;; Revision 1.2  2000/05/05 21:35:16  simsek
-;;; o Fixed row-vector-p and col-vector-p
-;;;
-;;; Revision 1.1  2000/04/14 00:11:12  simsek
-;;; o This file is adapted from obsolete files 'matrix-float.lisp'
-;;;   'matrix-complex.lisp' and 'matrix-extra.lisp'
-;;; o Initial revision.
-;;;
-;;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+(in-package :matlisp)
 
-;;; Definitions of STANDARD-MATRIX, REAL-MATRIX, COMPLEX-MATRIX.
-
-(in-package "MATLISP")
-
-#+nil (export '(real-matrix
-	  complex-matrix
-	  standard-matrix
-	  real-matrix-element-type
-	  real-matrix-store-type
-	  complex-matrix-element-type
-	  complex-matrix-store-type
-	  #|
-	  n
-	  m
-	  nxm
-	  |#
-	  nrows
-	  ncols
-	  number-of-elements
-	  row-vector-p
-	  col-vector-p
-	  row-or-col-vector-p
-	  square-matrix-p
-	  size
-	  fortran-matrix-indexing
-	  fortran-complex-matrix-indexing
-	  complex-coerce
-	  fill-matrix
-	  make-real-matrix-dim
-	  make-real-matrix
-	  make-complex-matrix-dim
-	  make-complex-matrix))
-
-(eval-when (load eval compile)
-(deftype integer4-matrix-element-type ()
-  '(signed-byte 32))
-
-(deftype real-matrix-element-type () 
-  "The type of the elements stored in a REAL-MATRIX"
-  'double-float)
-
-(deftype real-matrix-store-type (size) 
-  "The type of the storage structure for a REAL-MATRIX"
-  `(simple-array double-float ,size))
-
-(deftype complex-matrix-element-type () 
-  "The type of the elements stored in a COMPLEX-MATRIX"
-  'double-float)
-
-(deftype complex-matrix-store-type (size) 
-  "The type of the storage structure for a COMPLEX-MATRIX"
-  `(simple-array double-float ,size))
-)
-
-(declaim (ftype (function (standard-matrix) fixnum)
-		n
-		m
-		nxm
-		store-size)
-	 (ftype (function (real-matrix) fixnum)
-		n
-		m
-		nxm
-		store-size)
-	 (ftype (function (complex-matrix) fixnum)
-		n
-		m
-		nxm
-		store-size)
-	 (ftype (function (real-matrix) (simple-array double-float (*)))
-		 store)
-	 (ftype (function (complex-matrix) (simple-array double-float (*)))
-		 store))
-
-
-#|
-(defgeneric n (matrix)
-  (:documentation 
-"
-  Syntax
-  ======
-  (N matrix)
-
-  Purpose
-  =======
-  Returns the number of rows of MATRIX.
-"))
-
-(defgeneric m (matrix)
-  (:documentation
-"
-  Syntax
-  ======
-  (M matrix)
-
-  Purpose
-  =======
-  Returns the number of columns of MATRIX.
-"))
-
-(defgeneric nxm (matrix)
-  (:documentation
-"
-  Syntax
-  ======
-  (NxM matrix)
-
-  Purpose
-  =======
-  Returns the number of elements of MATRIX;
-  which is number of rows * number of columns. 
-"))
-|#
-
-(defgeneric store-size (matrix)
-  (:documentation
-"
-  Syntax
-  ======
-  (STORE-SIZE matrix)
-
-  Purpose
-  ======= 
-  Total number of elements needed to store the matrix.  (Usually
-  the same as (NxM matrix), but not necessarily so!
-"))
-
-(defgeneric store (matrix)
-  (:documentation
-"
-  Syntax
-  ======
-  (STORE matrix)
-
-  Purpose
-  =======
-The actual storage for the matrix.  It is typically a one dimensional
-array but not necessarily so.  The float and complex matrices do use
-1-D arrays.  The complex matrix actually stores the real and imaginary
-parts in successive elements of the matrix because Fortran stores them
-that way.
-"))
+;;
+(defclass standard-matrix (standard-tensor)  
+  ((rank
+    :accessor rank
+    :type index-type
+    :initform 2
+    :documentation "For a matrix, rank = 2."))
+  (:documentation "Basic matrix class."))
 
-#+(and (or cmu sbcl) gerds-pcl)
-(declaim (ext:slots (slot-boundp real-matrix complex-matrix)
-		    (inline standard-matri...
 
[truncated message content]