diff --git a/data/both.lisp b/data/both.lisp
index 2771d97a33dfbcf9e08daa87de0ddbe0df128997..70fa2de84ffd3e0160907965ec3f41c3ff54ad90 100644
--- a/data/both.lisp
+++ b/data/both.lisp
@@ -1,6 +1,6 @@
 ;; Functions for both vectors and matrices.
 ;; Liam Healy 2008-04-26 20:48:44EDT both.lisp
-;; Time-stamp: <2011-01-10 18:17:02EST both.lisp>
+;; Time-stamp: <2011-01-12 00:38:30EST both.lisp>
 ;;
 ;; Copyright 2008, 2009, 2010, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -35,7 +35,7 @@
   :export nil
   :documentation "Allocate memory for the GSL struct given a block pointer.")
 
-(defmfun alloc-from-block ((object matrix) blockptr)
+(defmfun alloc-from-block ((object grid:matrix) blockptr)
   ("gsl_" :category :type "_alloc_from_block")
   ((blockptr :pointer)
    (0 sizet)				; offset
@@ -113,7 +113,7 @@
     `(cffi:foreign-funcall
       ,(actual-gsl-function-name
 	`("gsl_" :category :type ,(if value "_set" "_get"))
-	(if matrixp 'matrix 'vector)
+	(if matrixp 'grid:matrix 'vector)
 	element-type)
       :pointer ,mpointer
       sizet ,(first indices)
@@ -192,7 +192,7 @@
   :documentation			; FDL
   "Add the scalar complex x to all the elements of array a.")
 
-(defmethod elt+ ((x float) (a foreign-array))
+(defmethod elt+ ((x float) (a grid:foreign-array))
   (elt+ a x))
   
 (defmfun elt- ((a both) (b both))
@@ -207,7 +207,7 @@
   "Subtract the elements of b from the elements of a.
    The two must have the same dimensions.")
 
-(defmethod elt- ((a foreign-array) (x float))
+(defmethod elt- ((a grid:foreign-array) (x float))
   (elt+ a (- x)))
 
 (defmfun elt* ((a vector) (b vector))
@@ -222,7 +222,7 @@
   "Multiply the elements of a by the elements of b.
    The two must have the same dimensions.")
 
-(defmfun elt* ((a matrix) (b matrix))
+(defmfun elt* ((a grid:matrix) (b grid:matrix))
   ("gsl_" :category :type "_mul_elements")
   (((mpointer a) :pointer) ((mpointer b) :pointer))
   :definition :methods
@@ -243,7 +243,7 @@
   "Divide the elements of a by the elements of b.
    The two must have the same dimensions.")
 
-(defmfun elt/ ((a matrix) (b matrix))
+(defmfun elt/ ((a grid:matrix) (b grid:matrix))
   ("gsl_" :category :type "_div_elements")
   (((mpointer a) :pointer) ((mpointer b) :pointer))
   :definition :methods
@@ -252,7 +252,7 @@
   :outputs (a)
   :return (a))
 
-(defmethod elt/ ((a foreign-array) (x number))
+(defmethod elt/ ((a grid:foreign-array) (x number))
   (elt* a (/ x)))
 
 (defmfun elt* ((a both) (x float))
@@ -279,7 +279,7 @@
   :documentation			; FDL
   "Multiply the elements of a by the scalar complex factor x.")
 
-(defmethod elt* ((x float) (a foreign-array))
+(defmethod elt* ((x float) (a grid:foreign-array))
   (elt* a x))
 
 ;;;;****************************************************************************
@@ -329,7 +329,7 @@
   "The index of the minimum value in a.  When there are several
   equal minimum elements, then the lowest index is returned.")
 
-(defmfun min-index ((a matrix))
+(defmfun min-index ((a grid:matrix))
   ("gsl_" :category :type "_min_index")
   (((mpointer a) :pointer) (imin (:pointer sizet)) (jmin (:pointer sizet)))
   :definition :methods
@@ -348,7 +348,7 @@
   "The index of the maximum value in a.  When there are several
   equal maximum elements, then the lowest index is returned.")
 
-(defmfun max-index ((a matrix))
+(defmfun max-index ((a grid:matrix))
   ("gsl_" :category :type "_max_index")
   (((mpointer a) :pointer) (imin (:pointer sizet)) (jmin (:pointer sizet)))
   :definition :methods
@@ -369,7 +369,7 @@
   returned.  Returned indices are minimum, maximum; for matrices
   imin, jmin, imax, jmax.")
 
-(defmfun minmax-index ((a matrix))
+(defmfun minmax-index ((a grid:matrix))
   ("gsl_" :category :type "_minmax_index")
   (((mpointer a) :pointer)
    (imin (:pointer sizet)) (jmin (:pointer sizet))
diff --git a/data/foreign-array.lisp b/data/foreign-array.lisp
index e7b28dab452881f42fcc66941fb2aa64fea7fec1..8ebf9e681c340889213c6d815c982e0ef82ebe97 100644
--- a/data/foreign-array.lisp
+++ b/data/foreign-array.lisp
@@ -1,6 +1,6 @@
 ;; A grid:foreign-array with added metadata for GSL.
 ;; Liam Healy 2008-04-06 21:23:41EDT
-;; Time-stamp: <2011-01-10 17:54:07EST foreign-array.lisp>
+;; Time-stamp: <2011-01-12 00:34:42EST foreign-array.lisp>
 ;;
 ;; Copyright 2008, 2009, 2010, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -53,7 +53,7 @@
 	 ;; from the block to the vector/matrix; we must do that manually here
 	 (cffi:foreign-slot-value
 	  array-struct
-	  (if (typep object 'matrix) 'gsl-matrix-c 'gsl-vector-c) 'data)
+	  (if (typep object 'grid:matrix) 'gsl-matrix-c 'gsl-vector-c) 'data)
 	 (grid:foreign-pointer object))
 	(tg:finalize object
 		     (lambda ()
@@ -82,20 +82,20 @@
   (let* ((category
 	  (case category-or-rank
 	    ((vector :vector 1) 'vector)
-	    ((matrix :matrix 2) 'matrix)
+	    ((matrix grid:matrix :matrix 2) 'grid:matrix)
 	    (t (error "Unrecognized category ~a" category-or-rank))))
 	 (cstruct
 	  (case category
 	    (vector 'gsl-vector-c)
-	    (matrix 'gsl-matrix-c)))
+	    (grid:matrix 'gsl-matrix-c)))
 	 (fa
 	  (grid:make-grid
 	   (case category
 	     (vector
-		`((foreign-array ,(cffi:foreign-slot-value mpointer cstruct 'size))
+		`((grid:foreign-array ,(cffi:foreign-slot-value mpointer cstruct 'size))
 		  ,element-type))
-	     (matrix
-		`((foreign-array
+	     (grid:matrix
+		`((grid:foreign-array
 		   ,(cffi:foreign-slot-value mpointer cstruct 'size0)
 		   ,(cffi:foreign-slot-value mpointer cstruct 'size1))
 		  ,element-type)))
diff --git a/data/matrix.lisp b/data/matrix.lisp
index eac02936bf7a365a7da5075af26f333e28c72646..22df33a2171efefae0a62f92421c7a6dda3fa4d3 100644
--- a/data/matrix.lisp
+++ b/data/matrix.lisp
@@ -1,6 +1,6 @@
 ;; Matrices
 ;; Liam Healy 2008-04-15 21:57:52EDT matrix.lisp
-;; Time-stamp: <2011-01-10 18:17:17EST matrix.lisp>
+;; Time-stamp: <2011-01-12 00:34:45EST matrix.lisp>
 ;;
 ;; Copyright 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -26,7 +26,7 @@
 ;;;; Mathematical
 ;;;;****************************************************************************
 
-(defmfun set-identity ((matrix matrix))
+(defmfun set-identity ((matrix grid:matrix))
   ("gsl_matrix" :type "_set_identity")
   (((mpointer matrix) :pointer))
   :definition :generic
@@ -43,8 +43,7 @@
 ;;;; Copying rows and columns
 ;;;;****************************************************************************
 
-(defmfun row
-    ((matrix matrix) i
+(defmfun row ((matrix grid:matrix) i
      &optional (vector (grid:make-foreign-array element-type :dimensions (dim1 matrix))))
   ("gsl_matrix" :type "_get_row")
   (((mpointer vector) :pointer) ((mpointer matrix) :pointer) (i sizet))
@@ -57,7 +56,7 @@
    into the vector.  The length of the vector must be the
    same as the length of the row.")
 
-(defmfun (setf row) ((vector vector) (matrix matrix) i)
+(defmfun (setf row) ((vector vector) (matrix grid:matrix) i)
   ("gsl_matrix" :type "_set_row")
   (((mpointer matrix) :pointer) (i sizet) ((mpointer vector) :pointer))
   :definition :generic
@@ -68,8 +67,7 @@
   "Copy the elements of the vector into the jth row of the matrix.
   The length of the vector must be the same as the length of the row.")
 
-(defmfun column
-    ((matrix matrix) i
+(defmfun column ((matrix grid:matrix) i
      &optional (vector (grid:make-foreign-array element-type :dimensions (dim0 matrix))))
   ("gsl_matrix" :type "_get_col")
   (((mpointer vector) :pointer) ((mpointer matrix) :pointer) (i sizet))
@@ -82,7 +80,7 @@
    into the vector.  The length of the vector must be the
    same as the length of the column.")
 
-(defmfun (setf column) ((vector vector) (matrix matrix) i)
+(defmfun (setf column) ((vector vector) (matrix grid:matrix) i)
   ("gsl_matrix" :type "_set_col")
   (((mpointer matrix) :pointer) (i sizet) ((mpointer vector) :pointer))
   :definition :generic
@@ -97,7 +95,7 @@
 ;;;; Exchanging rows and columns
 ;;;;****************************************************************************
 
-(defmfun swap-rows ((matrix matrix) i j)
+(defmfun swap-rows ((matrix grid:matrix) i j)
   ("gsl_matrix" :type "_swap_rows")
   (((mpointer matrix) :pointer) (i sizet) (j sizet))
   :definition :generic
@@ -107,7 +105,7 @@
   :documentation 			; FDL
   "Exchange the ith and jth rows of the matrix in-place.")
 
-(defmfun swap-columns ((matrix matrix) i j)
+(defmfun swap-columns ((matrix grid:matrix) i j)
   ("gsl_matrix" :type "_swap_columns")
   (((mpointer matrix) :pointer) (i sizet) (j sizet))
   :definition :generic
@@ -117,7 +115,7 @@
   :documentation 			; FDL
   "Exchange the ith and jth columns of the matrix in-place.")
 
-(defmfun swap-row-column ((matrix matrix) i j)
+(defmfun swap-row-column ((matrix grid:matrix) i j)
   ("gsl_matrix" :type "_swap_rowcol")
   (((mpointer matrix) :pointer) (i sizet) (j sizet))
   :definition :generic
@@ -129,7 +127,7 @@
    matrix in-place.  The matrix must be square for this operation to
    be possible.")
 
-(defmfun matrix-transpose* ((matrix matrix))
+(defmfun matrix-transpose* ((matrix grid:matrix))
   ("gsl_matrix" :type "_transpose")
   (((mpointer matrix) :pointer))
   :definition :generic
@@ -142,7 +140,7 @@
    operation to be possible.")
 
 (defmfun matrix-transpose
-    ((source matrix)
+    ((source grid:matrix)
      &optional
      (destination
       (grid:make-foreign-array element-type :dimensions (reverse (grid:dimensions source)))))
diff --git a/eigensystems/generalized.lisp b/eigensystems/generalized.lisp
index c6778ef3ad204bf1fa5a86ced689c37f82810c60..466d48fc334291772a96fcaad63b5c128c117533 100644
--- a/eigensystems/generalized.lisp
+++ b/eigensystems/generalized.lisp
@@ -1,6 +1,6 @@
 ;; Real generalized eigensystems
 ;; Liam Healy 2009-02-16 12:55:04EST real-generalized.lisp
-;; Time-stamp: <2011-01-10 18:17:27EST generalized.lisp>
+;; Time-stamp: <2011-01-12 00:08:22EST generalized.lisp>
 ;;
 ;; Copyright 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -56,8 +56,7 @@
     n-by-n complex generalized hermitian-definite eigensystems. The
     size of the workspace is O(5n).")
 
-(defmfun eigenvalues-gensymm
-    ((A matrix) (B matrix)
+(defmfun eigenvalues-gensymm ((A grid:matrix) (B grid:matrix)
      &optional
      (eigenvalues (grid:make-foreign-array element-type :dimensions (dim0 A)))
      (ws (eltcase complex (make-eigen-genherm (dim0 A))
@@ -78,8 +77,7 @@
   outlined above. On output, B contains its Cholesky decomposition and
   A is destroyed.")
 
-(defmfun eigenvalues-eigenvectors-gensymm
-    ((A matrix) (B matrix)
+(defmfun eigenvalues-eigenvectors-gensymm ((A grid:matrix) (B grid:matrix)
      &optional
      (eigenvalues (grid:make-foreign-array element-type :dimensions (dim0 A)))
      (eigenvectors (grid:make-foreign-array element-type :dimensions (grid:dimensions A)))
diff --git a/eigensystems/symmetric-hermitian.lisp b/eigensystems/symmetric-hermitian.lisp
index abed4dd5eae13394783de824c027a45e58900307..01534f333262e235ce18e9cfe4fe811de1aa8f8e 100644
--- a/eigensystems/symmetric-hermitian.lisp
+++ b/eigensystems/symmetric-hermitian.lisp
@@ -1,6 +1,6 @@
 ;; Eigenvectors and eigenvalues
 ;; Liam Healy, Sun May 21 2006 - 19:52
-;; Time-stamp: <2011-01-10 18:19:13EST symmetric-hermitian.lisp>
+;; Time-stamp: <2011-01-12 00:09:43EST symmetric-hermitian.lisp>
 ;;
 ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -66,8 +66,7 @@
 ;;;; Eigenvalues and eigenvectors
 ;;;;****************************************************************************
 
-(defmfun eigenvalues
-    ((A matrix)
+(defmfun eigenvalues ((A grid:matrix)
      &optional
      (eigenvalues (grid:make-foreign-array element-type :dimensions (dim0 A)))
      (ws (eltcase complex (make-eigen-herm (dim0 A))
@@ -91,8 +90,7 @@
   referenced.  The eigenvalues are stored in the vector eigenvalues and
   are unordered.")
 
-(defmfun eigenvalues-eigenvectors
-    ((A matrix)
+(defmfun eigenvalues-eigenvectors ((A grid:matrix)
      &optional
      (eigenvalues (grid:make-foreign-array element-type :dimensions (dim0 A)))
      (eigenvectors (grid:make-foreign-array element-type :dimensions (grid:dimensions A)))
@@ -126,7 +124,7 @@
 ;;;;****************************************************************************
 
 (defmfun sort-eigenvalues-eigenvectors
-    ((eigenvalues vector) (eigenvectors matrix) sort-type)
+    ((eigenvalues vector) (eigenvectors grid:matrix) sort-type)
   (double-float "gsl_eigen_symmv_sort"
    complex-double-float "gsl_eigen_hermv_sort")
   (((mpointer eigenvalues) :pointer) ((mpointer eigenvectors) :pointer)
diff --git a/fast-fourier-transforms/backward.lisp b/fast-fourier-transforms/backward.lisp
index 637ca42d93fd46c3230ad1c514c86b1e9f511d29..a5bf3c7b4ff37353fe783d20a76e915c920535a8 100644
--- a/fast-fourier-transforms/backward.lisp
+++ b/fast-fourier-transforms/backward.lisp
@@ -1,6 +1,6 @@
 ;; Backward FFT
 ;; Sumant Oemrawsingh, Sat Oct 24 2009 - 12:55
-;; Time-stamp: <2011-01-10 17:59:22EST backward.lisp>
+;; Time-stamp: <2011-01-11 23:34:19EST backward.lisp>
 ;;
 ;; Copyright 2009, 2011 Sumant Oemrawsingh, Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -138,13 +138,13 @@
     (remf pass-on-args :decimation-in-frequency)
     (remf pass-on-args :non-radix-2)
     (if (and (not non-radix-2) (power-of-2-p (floor (size vector) stride)))
-	(if (subtypep (element-type vector) 'real)
+	(if (subtypep (grid:element-type vector) 'real)
 	    (apply 'backward-fourier-transform-halfcomplex-radix2
 		   vector pass-on-args)
 	    (if decimation-in-frequency
 		(apply 'backward-fourier-transform-dif-radix2 vector pass-on-args)
 		(apply 'backward-fourier-transform-radix2 vector pass-on-args)))
-	(if (subtypep (element-type vector) 'real)
+	(if (subtypep (grid:element-type vector) 'real)
 	    (apply 'backward-fourier-transform-halfcomplex-nonradix2
 		   vector pass-on-args)
 	    (apply 'backward-fourier-transform-nonradix2 vector pass-on-args)))))
diff --git a/fast-fourier-transforms/example.lisp b/fast-fourier-transforms/example.lisp
index a05d5d3bfccc723724e07c69572e8ea4adbd1342..d1ed0bf7f25058397a5207ed82789b07858b8611 100644
--- a/fast-fourier-transforms/example.lisp
+++ b/fast-fourier-transforms/example.lisp
@@ -1,8 +1,8 @@
 ;; Example FFT: transform a pulse (using the "clean" fft interface)
 ;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24
-;; Time-stamp: <2010-09-04 08:44:35EDT example.lisp>
+;; Time-stamp: <2011-01-11 23:32:17EST example.lisp>
 ;;
-;; Copyright 2009 Sumant Oemrawsingh, Liam M. Healy
+;; Copyright 2009, 2010, 2011 Sumant Oemrawsingh, Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
 ;;
 ;; This program is free software: you can redistribute it and/or modify
@@ -107,7 +107,7 @@
   "The real vector consisting of the real part of the complex vector."
   (let ((real-vector
           (make-and-init-vector
-            (grid:component-float-type (element-type complex-vector))
+            (grid:component-float-type (grid:element-type complex-vector))
             (size complex-vector)
             :init-offset init-offset)))
     (loop for i below (size complex-vector) by stride
@@ -119,19 +119,20 @@
 (defun copy-with-stride (vector &key (stride 1) init-offset)
   "Copy a vector and initialize it."
   (let ((vec
-          (make-and-init-vector (element-type vector)
-                                (size vector)
-                                :init-offset init-offset)))
+	 (make-and-init-vector
+	  (grid:element-type vector)
+	  (size vector)
+	  :init-offset init-offset)))
     (loop for i below (size vector) by stride
           do
-          (setf (grid:gref vec i) (grid:gref vector i)))
+       (setf (grid:gref vec i) (grid:gref vector i)))
     vec))
 
 (defun size-vector-scalar (vector &key (stride 1))
   "Return the size of a vector while taking the stride into account."
   (coerce (floor (size vector) stride)
-	  (if (subtypep (element-type vector) 'complex)
-	      (element-type vector)
+	  (if (subtypep (grid:element-type vector) 'complex)
+	      (grid:element-type vector)
 	      'double-float)))
 
 #+nil
@@ -139,7 +140,7 @@
   (elt/ vector (size-vector-scalar vector :stride stride)))
 
 (defun vector/length (vector &key (stride 1))
-  (let ((element-type (element-type vector)))
+  (let ((element-type (grid:element-type vector)))
     (loop with length = (size-vector-scalar vector :stride stride)
           for i from 0 below (size vector) by stride
           do
diff --git a/fast-fourier-transforms/inverse.lisp b/fast-fourier-transforms/inverse.lisp
index 590f2b748ed375b95be9103c49993796dca11a58..27c89e07d441758a2ab144def96ed7408c6a8464 100644
--- a/fast-fourier-transforms/inverse.lisp
+++ b/fast-fourier-transforms/inverse.lisp
@@ -1,6 +1,6 @@
 ;; Inverse FFT
 ;; Sumant Oemrawsingh, Sat Oct 24 2009 - 12:55
-;; Time-stamp: <2011-01-10 17:59:24EST inverse.lisp>
+;; Time-stamp: <2011-01-11 23:34:17EST inverse.lisp>
 ;;
 ;; Copyright 2009, 2011 Sumant Oemrawsingh, Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -138,13 +138,13 @@
     (remf pass-on-args :decimation-in-frequency)
     (remf pass-on-args :non-radix-2)
     (if (and (not non-radix-2) (power-of-2-p (floor (size vector) stride)))
-	(if (subtypep (element-type vector) 'real)
+	(if (subtypep (grid:element-type vector) 'real)
 	    (apply 'inverse-fourier-transform-halfcomplex-radix2
 		   vector pass-on-args)
 	    (if decimation-in-frequency
 		(apply 'inverse-fourier-transform-dif-radix2 vector pass-on-args)
 		(apply 'inverse-fourier-transform-radix2 vector pass-on-args)))
-	(if (subtypep (element-type vector) 'real)
+	(if (subtypep (grid:element-type vector) 'real)
 	    (apply 'inverse-fourier-transform-halfcomplex-nonradix2
 		   vector pass-on-args)
 	    (apply 'inverse-fourier-transform-nonradix2 vector pass-on-args)))))
diff --git a/histogram/ntuple-example.dat b/histogram/ntuple-example.dat
index be4cb5a5d29dc9072987e90d6886995088febcb2..88b5a4aa490e8f25b6f4e8c29786ceb205c796de 100644
Binary files a/histogram/ntuple-example.dat and b/histogram/ntuple-example.dat differ
diff --git a/init/defmfun-array.lisp b/init/defmfun-array.lisp
index 603319014986aeb2286c49efde178f8b64b1c4a5..a8c52a8d06f1090d2a6bd9ddad788a03b5f50ff8 100644
--- a/init/defmfun-array.lisp
+++ b/init/defmfun-array.lisp
@@ -1,8 +1,8 @@
 ;; Helpers for defining GSL functions on arrays
 ;; Liam Healy 2009-01-07 22:01:16EST defmfun-array.lisp
-;; Time-stamp: <2010-07-12 12:28:57EDT defmfun-array.lisp>
+;; Time-stamp: <2011-01-12 00:22:59EST defmfun-array.lisp>
 ;;
-;; Copyright 2009 Liam M. Healy
+;; Copyright 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
 ;;
 ;; This program is free software: you can redistribute it and/or modify
@@ -28,7 +28,7 @@
     (defn name arglist gsl-name c-arguments categories key-args)
   (setf categories
 	(remove-if-not
-	 (lambda (c) (member c '(both matrix vector)))
+	 (lambda (c) (member c '(both grid:matrix vector)))
 	 categories))
   (if (member 'both categories)
       (progn 
@@ -45,12 +45,12 @@
 	   (optional-args-to-switch-gsl-functions arglist gsl-name))
 	  (copy-list key-args) 'vector)
 	 (generate-methods
-	  defn 'matrix
+	  defn 'grid:matrix
 	  name arglist gsl-name
 	  (actual-array-c-type
-	   'matrix c-arguments
+	   'grid:matrix c-arguments
 	   (optional-args-to-switch-gsl-functions arglist gsl-name))
-	  (copy-list key-args) 'matrix)))
+	  (copy-list key-args) 'grid:matrix)))
       ;; Generate forms for one category
       (generate-methods
        defn (first categories)
@@ -149,7 +149,7 @@
 			base-name))))))))))
 
 (defun actual-array-class (category element-type &optional replace-both)
-  "From the category ('vector, 'matrix, or 'both) and element type,
+  "From the category ('vector, 'grid:matrix, or 'both) and element type,
    find the class name."
   (case category
     (:element-type (grid:number-class element-type))
diff --git a/init/init.lisp b/init/init.lisp
index d5999ccc7c6f899e969a895564f53d66633b7bc4..c9f8fe771e71597317c1b9148b69962b23f5faf0 100644
--- a/init/init.lisp
+++ b/init/init.lisp
@@ -1,6 +1,6 @@
 ;; Load GSL
 ;; Liam Healy Sat Mar  4 2006 - 18:53
-;; Time-stamp: <2011-01-11 19:06:42EST init.lisp>
+;; Time-stamp: <2011-01-12 00:35:16EST init.lisp>
 ;;
 ;; Copyright 2006, 2007, 2008, 2009, 2010, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -21,10 +21,8 @@
 (defpackage gsll
   (:nicknames :gsl)
   (:use :common-lisp :cffi)
-  (:import-from
-   :grid
-    #:element-type #:foreign-array #:matrix #:dim0 #:dim1 #:^ #:copy)
-  (:export #:element-type #:dim0 #:dim1  #:copy))
+  (:import-from :grid #:dim0 #:dim1 #:^ #:copy)
+  (:export #:dim0 #:dim1 #:copy))
 
 ;;; Where there is a symbol conflict, take the other one.
 (shadow '#:row :antik-user) ; conflict with grid:row; they are equivalent
diff --git a/linear-algebra/blas2.lisp b/linear-algebra/blas2.lisp
index b2281333e5961d817f4c5fd95d95760799c9c90b..5b9983d24e7939148e3b73293b5386f2cf64d237 100644
--- a/linear-algebra/blas2.lisp
+++ b/linear-algebra/blas2.lisp
@@ -1,6 +1,6 @@
 ;; BLAS level 2, Matrix-vector operations
 ;; Liam Healy, Wed Apr 26 2006 - 21:08
-;; Time-stamp: <2011-01-10 18:19:10EST blas2.lisp>
+;; Time-stamp: <2011-01-12 00:42:10EST blas2.lisp>
 ;;
 ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -41,13 +41,13 @@
 ;;;;****************************************************************************
 
 (defun matrix-product-dimensions (a b)
-  (if (typep b 'matrix)
+  (if (typep b 'grid:matrix)
       (list (first (grid:dimensions a))
 	    (second (grid:dimensions b)))
       (first (grid:dimensions a))))
 
 (defmfun matrix-product
-    ((A matrix) (x vector)
+    ((A grid:matrix) (x vector)
      &optional
      y
      (alpha 1) (beta 1) (TransA :notrans) TransB
@@ -74,7 +74,7 @@
   parameter TransB.")
 
 (defmfun matrix-product-triangular
-    ((A matrix) (x vector)
+    ((A grid:matrix) (x vector)
      &optional (alpha 1) (uplo :upper) (TransA :notrans)
      (diag :nonunit) (side :left))
   ("gsl_blas_" :type "trmv")
@@ -106,7 +106,7 @@
    are taken as unity and are not referenced.")
 
 (defmfun inverse-matrix-product
-    ((A matrix) (x vector)
+    ((A grid:matrix) (x vector)
      &optional (alpha 1) (uplo :upper) (TransA :notrans)
      (diag :nonunit) (side :left))
   ("gsl_blas_" :type "trsv")
@@ -137,7 +137,7 @@
    matrix A are taken as unity and are not referenced.")
 
 (defmfun matrix-product-symmetric
-    ((A matrix) (x vector)
+    ((A grid:matrix) (x vector)
      &optional y
      (alpha 1) (beta 1) (uplo :upper) (side :left)
      &aux
@@ -170,7 +170,7 @@
 
 #+fsbv
 (defmfun matrix-product-hermitian
-    ((A matrix) (x vector)
+    ((A grid:matrix) (x vector)
      &optional
      (y (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A x)
 		     :initial-element 0))
@@ -199,7 +199,7 @@
   are used. The imaginary elements of the diagonal are automatically
   set to zero.")
 
-(defmfun rank-1-update (alpha (x vector) (y vector) (A matrix))
+(defmfun rank-1-update (alpha (x vector) (y vector) (A grid:matrix))
   ("gsl_blas_" :type "ger" :suffix)
   ((alpha :element-c-type) ((mpointer x) :pointer)
    ((mpointer y) :pointer) ((mpointer A) :pointer))
@@ -211,7 +211,7 @@
    "The rank-1 update A = alpha x y^T + A of the matrix A.")
 
 #+fsbv
-(defmfun conjugate-rank-1-update (alpha (x vector) (y vector) (A matrix))
+(defmfun conjugate-rank-1-update (alpha (x vector) (y vector) (A grid:matrix))
   ("gsl_blas_" :type "gerc")
   ((alpha :element-c-type) ((mpointer x) :pointer)
    ((mpointer y) :pointer) ((mpointer A) :pointer))
@@ -223,7 +223,7 @@
   "The conjugate rank-1 update A = alpha x y^H + A of the matrix A.")
 
 (defmfun symmetric-rank-1-update
-    ((x vector) (A matrix)
+    ((x vector) (A grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "syr")
   ((uplo cblas-uplo) (alpha :element-c-type)
@@ -249,7 +249,7 @@
 
 #+fsbv
 (defmfun hermitian-rank-1-update
-    ((x vector) (A matrix)
+    ((x vector) (A grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "her")
   ((uplo cblas-uplo) (alpha :element-c-type)
@@ -276,7 +276,7 @@
   elements of the diagonal are automatically set to zero.")
 
 (defmfun symmetric-rank-2-update
-    ((x vector) (y vector) (A matrix)
+    ((x vector) (y vector) (A grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "syr2")
   ((uplo cblas-uplo) (alpha :element-c-type)
@@ -303,7 +303,7 @@
 
 #+fsbv
 (defmfun hermitian-rank-2-update
-    ((x vector) (y vector) (A matrix)
+    ((x vector) (y vector) (A grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "her2")
   ((uplo cblas-uplo) (alpha :element-c-type)
diff --git a/linear-algebra/blas3.lisp b/linear-algebra/blas3.lisp
index ab51b2a892fd69cf74bfe0099dcd332d24c98fcd..933448596d23930f85d2b1eb568f74ad87e1040d 100644
--- a/linear-algebra/blas3.lisp
+++ b/linear-algebra/blas3.lisp
@@ -1,8 +1,8 @@
 ;; BLAS level 3, Matrix-matrix operations
 ;; Liam Healy, Wed Apr 26 2006 - 21:08
-;; Time-stamp: <2010-07-07 14:25:00EDT blas3.lisp>
+;; Time-stamp: <2011-01-12 00:49:27EST blas3.lisp>
 ;;
-;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
+;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
 ;;
 ;; This program is free software: you can redistribute it and/or modify
@@ -34,7 +34,7 @@
 ;;;;****************************************************************************
 
 (defmfun matrix-product
-    ((A matrix) (B matrix)
+    ((A grid:matrix) (B grid:matrix)
      &optional
      C
      (alpha 1) (beta 1) (TransA :notrans) (TransB :notrans)
@@ -53,7 +53,7 @@
   :outputs (Carr))
 
 (defmfun matrix-product-symmetric
-    ((A matrix) (B matrix)
+    ((A grid:matrix) (B grid:matrix)
      &optional C (alpha 1) (beta 1) (uplo :upper) (side :left)
      &aux
      (Carr
@@ -71,7 +71,7 @@
 
 #+fsbv
 (defmfun matrix-product-hermitian
-    ((A matrix) (B matrix)
+    ((A grid:matrix) (B grid:matrix)
      &optional
      (C (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A B)
 		     :initial-element 0))
@@ -88,7 +88,7 @@
   :outputs (C))
 
 (defmfun matrix-product-triangular
-    ((A matrix) (B matrix)
+    ((A grid:matrix) (B grid:matrix)
      &optional (alpha 1) (uplo :upper) (TransA :notrans)
      (diag :nonunit) (side :left))
   ("gsl_blas_" :type "trmm")
@@ -101,7 +101,7 @@
   :outputs (B))
 
 (defmfun inverse-matrix-product
-    ((A matrix) (B matrix)
+    ((A grid:matrix) (B grid:matrix)
      &optional (alpha 1) (uplo :upper) (TransA :notrans)
      (diag :nonunit) (side :left))
   ;; This signals an error for complex arguments because you can't pass a
@@ -116,7 +116,7 @@
   :outputs (B))
 
 (defmfun symmetric-rank-1-update
-    ((A matrix) (C matrix)
+    ((A grid:matrix) (C grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "syrk")
   ((uplo cblas-uplo) (trans cblas-transpose)
@@ -129,7 +129,7 @@
 
 #+fsbv
 (defmfun hermitian-rank-1-update
-    ((A matrix) (C matrix)
+    ((A grid:matrix) (C grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "herk")
   ((uplo cblas-uplo) (trans cblas-transpose)
@@ -141,7 +141,7 @@
   :outputs (C))
 
 (defmfun symmetric-rank-2-update
-    ((A matrix) (B matrix) (C matrix)
+    ((A grid:matrix) (B grid:matrix) (C grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "syr2k")
   ((uplo cblas-uplo) (trans cblas-transpose)
@@ -155,7 +155,7 @@
 
 #+fsbv
 (defmfun hermitian-rank-2-update
-    ((A matrix) (B matrix) (C matrix)
+    ((A grid:matrix) (B grid:matrix) (C grid:matrix)
      &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans))
   ("gsl_blas_" :type "her2k")
   ((uplo cblas-uplo) (trans cblas-transpose)
diff --git a/linear-algebra/cholesky.lisp b/linear-algebra/cholesky.lisp
index d41918ea772d09f35362c8a249b55b416d117bae..253474c8c5da3b7419471f523e773030400c4a75 100644
--- a/linear-algebra/cholesky.lisp
+++ b/linear-algebra/cholesky.lisp
@@ -1,6 +1,6 @@
 ;; Cholesky Decomposition
 ;; Liam Healy, Wed May  3 2006 - 16:38
-;; Time-stamp: <2011-01-10 18:19:09EST cholesky.lisp>
+;; Time-stamp: <2011-01-12 00:29:25EST cholesky.lisp>
 ;;
 ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -34,7 +34,7 @@
 
 ;;; GSL version 1.10 introduced functions for complex matrices.
 
-(defmfun cholesky-decomposition ((A matrix))
+(defmfun cholesky-decomposition ((A grid:matrix))
   ("gsl_linalg" :complex "_cholesky_decomp")
   (((mpointer A) :pointer))
   :definition :generic
@@ -51,7 +51,7 @@
   returning the error input-domain.")
 
 (defmfun cholesky-solve
-    ((A matrix) (b vector) &optional x-spec
+    ((A grid:matrix) (b vector) &optional x-spec
      &aux
      (x (grid:make-foreign-array-or-default x-spec (grid:dimensions b) t)))
   (("gsl_linalg" :complex "_cholesky_svx")
diff --git a/linear-algebra/diagonal.lisp b/linear-algebra/diagonal.lisp
index cace5ac3e6b2ef4647f29c9fb9aedfc26464de4e..7df41dd5b8831acf6501cf43b1c279fea8fba142 100644
--- a/linear-algebra/diagonal.lisp
+++ b/linear-algebra/diagonal.lisp
@@ -1,8 +1,8 @@
 ;; Tridiagonal and Bidiagonal matrices
 ;; Liam Healy, Thu May  4 2006 - 15:43
-;; Time-stamp: <2010-06-30 19:57:28EDT diagonal.lisp>
+;; Time-stamp: <2011-01-12 00:49:45EST diagonal.lisp>
 ;;
-;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
+;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
 ;;
 ;; This program is free software: you can redistribute it and/or modify
@@ -35,7 +35,7 @@
 ;;; A = U T U^T where U is a unitary
 ;;; matrix and T is a real symmetric tridiagonal matrix.
 
-(defmfun tridiagonal-decomposition ((A matrix) tau)
+(defmfun tridiagonal-decomposition ((A grid:matrix) tau)
   (double-float "gsl_linalg_symmtd_decomp"
    complex-double-float "gsl_linalg_hermtd_decomp")
   (((mpointer A) :pointer) ((mpointer tau) :pointer))
@@ -53,7 +53,7 @@
    Q.  This storage scheme is the same as used by lapack.  The
    upper triangular part of A is not referenced.")
 
-(defmfun tridiagonal-unpack ((A matrix) tau Q diag subdiag)
+(defmfun tridiagonal-unpack ((A grid:matrix) tau Q diag subdiag)
   (double-float "gsl_linalg_symmtd_unpack"
    complex-double-float "gsl_linalg_hermtd_unpack")
   (((mpointer A) :pointer) ((mpointer tau) :pointer)
@@ -69,7 +69,7 @@
   orthogonal or unitary matrix Q, the vector of diagonal elements diag
   and the real vector of subdiagonal elements subdiag.")
 
-(defmfun tridiagonal-unpack-T ((A matrix) diag subdiag)
+(defmfun tridiagonal-unpack-T ((A grid:matrix) diag subdiag)
   (double-float "gsl_linalg_symmtd_unpack_T"
    complex-double-float "gsl_linalg_hermtd_unpack_T")
   (((mpointer A) :pointer) ((mpointer diag) :pointer)
diff --git a/linear-algebra/lu.lisp b/linear-algebra/lu.lisp
index 58392752621460744509a9eca8d71068f03d6378..02b0a5e8b3a8352cc277a42ea34895a74a4469c4 100644
--- a/linear-algebra/lu.lisp
+++ b/linear-algebra/lu.lisp
@@ -1,6 +1,6 @@
 ;; LU decomposition
 ;; Liam Healy, Thu Apr 27 2006 - 12:42
-;; Time-stamp: <2011-01-10 18:16:19EST lu.lisp>
+;; Time-stamp: <2011-01-12 00:44:43EST lu.lisp>
 ;;
 ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -23,7 +23,7 @@
 ;;; /usr/include/gsl/gsl_linalg.h
 
 (defmfun LU-decomposition
-    ((A matrix) &optional (permutation (make-permutation (dim0 A))))
+    ((A grid:matrix) &optional (permutation (make-permutation (dim0 A))))
   ("gsl_linalg" :complex "_LU_decomp")
   (((mpointer A) :pointer) ((mpointer permutation) :pointer)
    (signum (:pointer :int)))
@@ -53,7 +53,7 @@
   Algorithm 3.4.1).")
 
 (defmfun LU-solve
-    ((A matrix) (b vector) permutation &optional x-spec
+    ((A grid:matrix) (b vector) permutation &optional x-spec
      &aux
      (x (if (eq x-spec t)
 	    (grid:make-foreign-array element-type :dimensions (grid:dimensions b))
@@ -78,7 +78,7 @@
   will be computed there.  Otherwise it should be a supplied vector.")
 
 (defmfun LU-refine
-    ((A matrix) LU p (b vector) (x vector)
+    ((A grid:matrix) LU p (b vector) (x vector)
      &optional (residual (grid:make-foreign-array element-type :dimensions (dim0 A))))
   ("gsl_linalg" :complex "_LU_refine")
   (((mpointer A) :pointer) ((mpointer LU) :pointer)
@@ -95,7 +95,7 @@
   A x = b, using the LU decomposition of A into (LU,p). The initial
   residual r = A x - b is also computed and stored in residual. ")
 
-(defmfun LU-invert ((LU matrix) p inverse)
+(defmfun LU-invert ((LU grid:matrix) p inverse)
   ("gsl_linalg" :complex "_LU_invert")
   (((mpointer LU) :pointer) ((mpointer p) :pointer)
    ((mpointer inverse) :pointer))
@@ -112,7 +112,7 @@
    the same result more efficiently and reliably (consult any
    introductory textbook on numerical linear algebra for details).")
 
- (defmfun LU-determinant ((LU matrix) signum)
+ (defmfun LU-determinant ((LU grid:matrix) signum)
   ("gsl_linalg" :complex "_LU_det")
   (((mpointer LU) :pointer) (signum :int))
   :definition :generic
@@ -124,7 +124,7 @@
   decomposition, LU. The determinant is computed as the product of the
   diagonal elements of U and the sign of the row permutation signum.")
 
-(defmfun LU-log-determinant ((LU matrix))
+(defmfun LU-log-determinant ((LU grid:matrix))
   ("gsl_linalg" :complex "_LU_lndet")
   (((mpointer LU) :pointer))
   :definition :generic
@@ -137,7 +137,7 @@
    LU. This function may be useful if the direct computation of the
    determinant would overflow or underflow.")
 
-(defmfun LU-sgndet ((LU matrix) signum)
+(defmfun LU-sgndet ((LU grid:matrix) signum)
   ("gsl_linalg" :complex "_LU_sgndet")
   (((mpointer LU) :pointer) (signum :int))
   :definition :generic
@@ -170,7 +170,7 @@
    the index."
   (let* ((dim (dim0 matrix))
 	 (rhs (or vector
-		  (create-rhs-vector dim (element-type matrix)))))
+		  (create-rhs-vector dim (grid:element-type matrix)))))
     (multiple-value-bind (upper permutation signum)
 	(LU-decomposition (copy matrix))
       (declare (ignore signum))
diff --git a/linear-algebra/matrix-generation.lisp b/linear-algebra/matrix-generation.lisp
index 21645764859d726b5456943ddec4043a47f99330..902564b9e657505564269e965159dd85e71b2494 100644
--- a/linear-algebra/matrix-generation.lisp
+++ b/linear-algebra/matrix-generation.lisp
@@ -1,8 +1,8 @@
 ;; Generate matrices used in tests of linear algebra functions
 ;; Liam Healy 2009-09-19 18:28:31EDT matrix-generation.lisp
-;; Time-stamp: <2010-08-12 22:48:01EDT matrix-generation.lisp>
+;; Time-stamp: <2011-01-11 23:58:19EST matrix-generation.lisp>
 ;;
-;; Copyright 2009, 2010 Liam M. Healy
+;; Copyright 2009, 2010, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
 ;;
 ;; This program is free software: you can redistribute it and/or modify
@@ -40,13 +40,13 @@
    :source function
    :destination-specification
    (if dim1
-       `((foreign-array ,dim0 ,dim1) ,element-type)
-       `((foreign-array ,dim0) ,element-type))))
+       `((grid:foreign-array ,dim0 ,dim1) ,element-type)
+       `((grid:foreign-array ,dim0) ,element-type))))
 
 (defun constant-matrix
     (constant dim0 &optional (dim1 dim0) (element-type 'double-float))
   (grid:make-grid
-   `((foreign-array ,dim0 ,dim1) ,element-type)
+   `((grid:foreign-array ,dim0 ,dim1) ,element-type)
    :initial-element constant))
 
 ;;;;****************************************************************************
diff --git a/random/shuffling-sampling.lisp b/random/shuffling-sampling.lisp
index b41d27333553e80ae7e655df73c39bf43212c43a..47756947a673f62536eb43f4d99a2e2ee0e6c802 100644
--- a/random/shuffling-sampling.lisp
+++ b/random/shuffling-sampling.lisp
@@ -1,6 +1,6 @@
 ;; Shuffling and sampling
 ;; Liam Healy, Sat Dec  2 2006 - 18:40
-;; Time-stamp: <2011-01-10 17:59:29EST shuffling-sampling.lisp>
+;; Time-stamp: <2011-01-11 23:34:14EST shuffling-sampling.lisp>
 ;;
 ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
@@ -46,7 +46,7 @@
      &aux
      (destarr
       (if (integerp dest)
-	  (grid:make-foreign-array (element-type src) :dimensions dest)
+	  (grid:make-foreign-array (grid:element-type src) :dimensions dest)
 	  dest)))
   "gsl_ran_choose"
   (((mpointer generator) :pointer)
@@ -74,7 +74,7 @@
      &aux
      (destarr
       (if (integerp dest)
-	  (grid:make-foreign-array (element-type src) :dimensions dest)
+	  (grid:make-foreign-array (grid:element-type src) :dimensions dest)
 	  dest)))
   "gsl_ran_sample"
   (((mpointer generator) :pointer)
diff --git a/tests/fast-fourier-transform.lisp b/tests/fast-fourier-transform.lisp
index dac5cda25708d500cfc06727e4576a1520da12f2..72d880a467ce14a95efbf337b7bf114a779b5b49 100644
--- a/tests/fast-fourier-transform.lisp
+++ b/tests/fast-fourier-transform.lisp
@@ -1,8 +1,8 @@
 ;; Fast fourier transform tests
 ;; Liam Healy 2010-08-14 11:58:26EDT fast-fourier-transform.lisp
-;; Time-stamp: <2010-09-04 17:07:32EDT fast-fourier-transform.lisp>
+;; Time-stamp: <2011-01-12 00:50:11EST fast-fourier-transform.lisp>
 ;;
-;; Copyright 2010 Liam M. Healy
+;; Copyright 2010, 2011 Liam M. Healy
 ;; Distributed under the terms of the GNU General Public License
 ;;
 ;; This program is free software: you can redistribute it and/or modify
@@ -27,7 +27,7 @@
   "Make a grid that is all of the original array except that the
    stride points are excluded."
   (grid:make-foreign-array
-   (element-type array)
+   (grid:element-type array)
    :initial-contents
    (loop for i from 0 below (grid:total-size array)
       unless (zerop (mod i stride))