diff --git a/basis-splines.lisp b/basis-splines.lisp
index 9bef5258cae4e5494d4f189d2693966d4d8eccb8..17bcf76569c33e44ec08eeecde4790b29df3efda 100644
--- a/basis-splines.lisp
+++ b/basis-splines.lisp
@@ -1,6 +1,6 @@
 ;; Basis splines.
 ;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
-;; Time-stamp: <2008-08-23 15:35:59EDT basis-splines.lisp>
+;; Time-stamp: <2008-11-30 23:43:58EST basis-splines.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -75,13 +75,13 @@
 	 (bw (basis-spline order nbreak))
 	 (mw (fit-workspace ndata ncoeffs))
 	 (rng (random-number-generator *mt19937* 0))
-	 (B (vector-double-float ncoeffs))
-	 (c (vector-double-float ncoeffs))
-	 (cov (matrix-double-float (list ncoeffs ncoeffs)))
-	 (w (vector-double-float ndata))
-	 (x (vector-double-float ndata))
-	 (y (vector-double-float ndata))
-	 (Xmatrix (matrix-double-float (list ndata ncoeffs)))
+	 (B (make-array* 'double-float :dimensions ncoeffs))
+	 (c (make-array* 'double-float :dimensions ncoeffs))
+	 (cov (make-array* 'double-float :dimensions (list ncoeffs ncoeffs)))
+	 (w (make-array* 'double-float :dimensions ndata))
+	 (x (make-array* 'double-float :dimensions ndata))
+	 (y (make-array* 'double-float :dimensions ndata))
+	 (Xmatrix (make-array* 'double-float :dimensions (list ndata ncoeffs)))
 	 (sigma 0.1d0))
     ;; The data to be fitted.
     (dotimes (i ndata)
diff --git a/data/array-tests.lisp b/data/array-tests.lisp
index 0baa88d36bffbab8380a8ff5e44d0148d294c97b..52a32e33511957da21dfc40b7e8c4f953c96ffce 100644
--- a/data/array-tests.lisp
+++ b/data/array-tests.lisp
@@ -1,6 +1,6 @@
 ;; Tests of array functions 
 ;; Liam Healy 2008-10-20 22:41:48EDT array-tests.lisp
-;; Time-stamp: <2008-11-16 10:28:21EST array-tests.lisp>
+;; Time-stamp: <2008-11-30 22:51:15EST array-tests.lisp>
 ;; $Id: $
 
 ;;; Generate each file with #'write-test-to-file, e.g.
@@ -24,8 +24,8 @@
 
 (generate-all-array-tests
  matrix-set-all-m+ :no-complex
- (letm ((m1 (array-default '(3 3) t))
-	(m2 (array-default '(3 3))))
+ (let ((m1 (array-default '(3 3) t))
+       (m2 (array-default '(3 3))))
    (set-all m1 (scalar-default))
    (cl-array (m+ m1 m2))))
 
diff --git a/data/combination.lisp b/data/combination.lisp
index a95a66c7feedfd59636fc5cd8d1dd2b7a0fb98db..2ab29c625b4ce514d2546b02cc9a960aee611687 100644
--- a/data/combination.lisp
+++ b/data/combination.lisp
@@ -1,6 +1,6 @@
 ;; Combinations
 ;; Liam Healy, Sun Mar 26 2006 - 11:51
-;; Time-stamp: <2008-11-29 15:09:07EST combination.lisp>
+;; Time-stamp: <2008-11-30 19:17:47EST combination.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -29,11 +29,9 @@
 
 (arglist-only combination "A combination." n k sync-array-on-exit)
 
-(defun make-data-combination (nk)
+(export 'make-combination)
+(defun make-combination (nk)
   "Make the combination object with the data array."
-  ;; The replaces the function of make-data-from-array and
-  ;; make-data-from-dimensions because the class needs different
-  ;; initialization than other kinds of data.
   (let ((k (second nk)))
     (make-instance
      'combination
diff --git a/data/data.lisp b/data/data.lisp
index 4bc3f1d02e4b82036fb56724a6b9dcdec8dfa6bc..a58a4f1df30a5bde6972a9a7c2468210565f9b56 100644
--- a/data/data.lisp
+++ b/data/data.lisp
@@ -1,7 +1,7 @@
 ;; "Data" is bulk arrayed data, like vectors, matrices, permutations,
 ;; combinations, or histograms.
 ;; Liam Healy 2008-04-06 21:23:41EDT
-;; Time-stamp: <2008-11-29 15:08:23EST data.lisp>
+;; Time-stamp: <2008-11-30 19:27:55EST data.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -80,7 +80,6 @@
 ;;;;****************************************************************************
 
 (defparameter *class-element-type* nil
-  ;;; This is needed for make-data-from-dimensions
   "The mapping between the class name and the CL element type.")
 
 (defun data-class-name (category element-type)
@@ -104,6 +103,7 @@
 		 ((element-type :initform ',element-type-cl
 				:allocation :class)))
 	       ;; Define the letm-expansion method
+	       ;;#+(or)			; REMOVE
 	       (defmethod letm-expansion
 		   (symbol (type (eql ',class-name)) args body)
 		 (expand-data symbol type (first args) (second args) body))
@@ -111,6 +111,7 @@
 	       (pushnew ',(cons class-name element-type-cl)
 		*class-element-type* :test #'equal)
 	       ;; Set up dummy function to provide arglist and documentation
+	       ;;#+(or)			; REMOVE
 	       (arglist-only
 		,class-name
 		,(format nil "A ~(~a~) of element type ~(~a~)."
@@ -123,29 +124,47 @@
 ;;;; Make data from either the dimensions provided or from the initial values
 ;;;;****************************************************************************
 
-(defun make-data-from-array (class array)
-  "Make the data object from the CL array make with make-array*."
-  (make-instance class
-		 :cl-array array
-		 :mpointer nil ; this will be set by :before method below.
-		 #-native :c-pointer #-native nil ; this will be set by defmfun
-		 :dimensions (array-dimensions array)
-		 :total-size (array-total-size array)))
-
-(defun make-data-from-dimensions (type dimensions)
-  (make-data-from-array
-   type
-   (make-array* dimensions (lookup-type type *class-element-type*))))
-
-(defun make-data (type array-or-dimensions)
-  (if (eq type 'combination)
-      (make-data-combination array-or-dimensions)
-      (if (typep array-or-dimensions 'array)
-	  (make-data-from-array type array-or-dimensions)
-	  (make-data-from-dimensions type array-or-dimensions))))
+(export 'make-array*)
+(defun make-array* (element-type &rest keys)
+  "Make a GSLL array with the given element type,
+   :dimensions, :initial-contents and/or :initial-element."
+  (let ((ffa (apply #'make-ffa element-type keys)))
+    (make-instance
+     (data-class-name
+      (if (eql (array-rank ffa) 2) 'matrix 'vector)
+      element-type)
+     :cl-array ffa
+     :mpointer nil
+     #-native :c-pointer #-native nil	; this will be set by defmfun
+     :dimensions (array-dimensions ffa)
+     :total-size (array-total-size ffa))))
+
+(defun hashm-numeric-code (n)
+  "Get the appropriate element type for the numeric code n"
+  (case n
+    ((nil 1) 'double-float)
+    (2 '(complex double-float))
+    (3 'single-float)
+    (4 '(complex single-float))
+    (7 '(signed-byte 8))
+    (8 '(unsigned-byte 8))
+    (15 '(signed-byte 16))
+    (16 '(unsigned-byte 16))
+    (31 '(signed-byte 32))
+    (32 '(unsigned-byte 32))
+    (63 '(signed-byte 64))
+    (64 '(unsigned-byte 64))))
+
+(set-dispatch-macro-character
+ #\# #\m
+ (lambda (stream subchar arg)
+   (declare (ignore subchar))
+   (read-char stream)
+   (let ((list (read-delimited-list #\) stream)))
+     `(make-array* ',(hashm-numeric-code arg) :initial-contents ',list))))
 
 ;;;;****************************************************************************
-;;;; Expand letm bindings
+;;;; Expand letm bindings  TO BE OBSOLETE
 ;;;;****************************************************************************
 
 (defun component-type (eltype)
@@ -199,6 +218,7 @@
   "Create an array (vector or matrix) from the literal contents."
   `(abody ,@(mapcar (lambda (e) `(quote ,e)) type-contents)))
 
+;;#+(or)			; REMOVE
 (defun expand-data (symbol type init-or-spec sync-exit body)
   "Expand the form for a gsl-data object.  The symbol is bound
    within the body to the object.  The argument
diff --git a/data/foreign-friendly.lisp b/data/foreign-friendly.lisp
index 47708326f2f996c63e3621745bfdcebe32a8cecc..13b7d0c210f0fa380c1c6e714ad91494afc816e2 100644
--- a/data/foreign-friendly.lisp
+++ b/data/foreign-friendly.lisp
@@ -1,6 +1,6 @@
 ;; Use the foreign-friendly arrays package.
 ;; Liam Healy 2008-03-22 15:40:08EDT
-;; Time-stamp: <2008-11-29 11:16:03EST foreign-friendly.lisp>
+;; Time-stamp: <2008-11-30 18:44:50EST foreign-friendly.lisp>
 ;; $Id$
 
 ;;; Foreign-friendly arrays (original implementation by Tamas Papp)
@@ -14,32 +14,58 @@
 ;;;; Make arrays for possible foreign use 
 ;;;;****************************************************************************
 
-(export '(make-array*))
-
-(defun make-array*
-    (dimensions element-type
-     &key (initial-element nil initial-element-p)
-     (initial-contents nil initial-contents-p))
+(defun make-ffa
+    (element-type
+     &key dimensions (initial-element nil initial-element-p)
+     (initial-contents nil initial-contents-p)
+     &allow-other-keys)
   "Make an array of one or two dimensions for possible use in foreign code.
    Syntax is similar to make-array, but note that element-type
    is mandatory and limited to certain types."
   (assert (member element-type *array-element-types* :test 'equal)
 	  (element-type)
 	  "Specified element-type must be one of *array-element-types*.")
-  (if initial-contents-p
-      (if (and (subtypep element-type 'complex)
-	       (= (length initial-contents)
-		  (* 2 (if (listp dimensions) (apply #'* dimensions) dimensions))))
-	  (make-ffa dimensions element-type
-			:initial-contents
-			(loop for (re im) on initial-contents by #'cddr
-			   collect (complex re im)))
-	  (make-ffa dimensions element-type :initial-contents initial-contents))
-      (if initial-element-p
-	  (make-ffa dimensions element-type :initial-element initial-element)
-	  (make-ffa dimensions element-type))))
-
-(defun make-ffa (dimensions element-type &key
+  (let* ((initial-matrix (listp (first initial-contents)))
+	 (complex-initial-real
+	  (and initial-contents-p
+	       (subtypep element-type 'complex)
+	       (typep
+		(if initial-matrix
+		    (caar initial-contents)
+		    (first initial-contents))
+		'real))))
+    (unless dimensions
+      (if initial-contents-p
+	  (setf dimensions
+		(if initial-matrix
+		    ;; matrix
+		    (list (length initial-contents)
+			  (ceiling (length (first initial-contents))
+				   (if complex-initial-real 2 1)))
+		    ;; vector
+		    (ceiling (length initial-contents)
+			     (if complex-initial-real 2 1))))
+	  (error "dimensions must be specified if contents are not")))
+    (when (and initial-matrix initial-contents-p)
+      ;; flatten matrix spec
+      (setf initial-contents (mapcar #'append initial-contents)))
+    (when complex-initial-real
+      (setf initial-contents (list-complex-from-real initial-contents)))
+    (apply 'make-ffa-1d
+	   dimensions element-type
+	   (append
+	    (when initial-contents-p
+		     (list :initial-contents initial-contents))
+	    (when initial-element-p
+		     (list :initial-element initial-element))))))
+
+(defun list-complex-from-real (list)
+  "Collect a list of complex numbers from real components."
+  (assert (evenp (length list)) (list) 
+	  "Cannot collect reals into complex from an odd number of reals.")
+  (loop for (re im) on list by #'cddr collect (complex re im)))
+
+(defun make-ffa-1d (dimensions element-type &key
 		 (initial-element 0 initial-element-p)
 		 (initial-contents nil initial-contents-p))
   "Make an array that is either one-dimensional or displaced to a
@@ -60,6 +86,9 @@
 							element-type)))
 		  ;; contents given, copy or coerce
 		  (initial-contents-p
+		   (when (listp (first initial-contents))
+		     ;; flatten list for matrices
+		     (setf initial-contents (apply 'append initial-contents)))
 		   (assert (= (length initial-contents) length))
 		   (if (typep initial-contents (list 'vector element-type))
 		       (copy-seq initial-contents)
diff --git a/data/permutation.lisp b/data/permutation.lisp
index 3607fe2b7fa12a3d8de9833e96f926c8f5f23931..28c6fae4b12e0d5db2b36738a539b6d4b765b914 100644
--- a/data/permutation.lisp
+++ b/data/permutation.lisp
@@ -1,6 +1,6 @@
 ;; Permutations
 ;; Liam Healy, Sun Mar 26 2006 - 11:51
-;; Time-stamp: <2008-10-25 18:52:27EDT permutation.lisp>
+;; Time-stamp: <2008-11-30 23:32:48EST permutation.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -287,7 +287,7 @@
    (swap-elements perm-1 1 3)
    (cl-array perm-1))
  (letm ((perm-1 (permutation 4 t))	;permute-vector
-	(intvec (vector-signed-byte-32 (a 11 22 33 44))))
+	(intvec #31m(11 22 33 44)))
    (set-identity perm-1)
    (swap-elements perm-1 1 3)
    (swap-elements perm-1 0 2)
diff --git a/documentation/documentation.html b/documentation/documentation.html
index e44c666d14d6b92d9ac2dec32eba3d76df151300..e5521903fd5b9d5cdd5c7a8fec5397182180543e 100644
--- a/documentation/documentation.html
+++ b/documentation/documentation.html
@@ -84,6 +84,7 @@ will be visible; however, the function should not be used outside a
 
  <h3>Arrays</h3>
  <div class="content">
+ <p>FIX: #m, #nm, <code>make-array*</code>
  <p>
  GSLL has many functions that work on vectors (one-dimensional arrays)
  and matrices (two-dimensional arrays).  GSLL supports all array element
@@ -204,7 +205,7 @@ gamma-randist, negative-binomial.
     <address><a href="mailto:gsll-devel@common-lisp.net">Liam Healy</a></address>
 <!-- hhmts start -->
     <small>
-      Time-stamp: <2008-11-15 21:53:03EST documentation.html>
+      Time-stamp: <2008-11-30 19:28:51EST documentation.html>
       </small>
 <!-- hhmts end -->
  </div>
diff --git a/eigensystems.lisp b/eigensystems.lisp
index 7a2ed5fa2b91fbf6576b837deba8772688b82482..a907e0696413cadae41c916c1f7d050811377951 100644
--- a/eigensystems.lisp
+++ b/eigensystems.lisp
@@ -1,6 +1,6 @@
 ;; Eigenvectors and eigenvalues
 ;; Liam Healy, Sun May 21 2006 - 19:52
-;; Time-stamp: <2008-10-25 11:34:13EDT eigensystems.lisp>
+;; Time-stamp: <2008-11-30 23:37:31EST eigensystems.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -158,13 +158,12 @@
 ;;;;****************************************************************************
  
 (defun eigenvalue-eigenvectors-example ()
-  (letm ((evecs (matrix-double-float '(3 3)))
-	 (evals (vector-double-float 3))
+  (letm ((evecs (make-array* 'double-float :dimensions '(3 3)))
+	 (evals (make-array* 'double-float :dimensions 3))
 	 (ws (eigen-symmv 3))
-	 (mat (matrix-double-float 
-	       (a (20.0d0 -10.0d0 0.0d0)
-		  (-10.0d0 30.0d0 0.0d0)
-		  (0.0d0 0.0d0 40.0d0)))))
+	 (mat #m((20.0d0 -10.0d0 0.0d0)
+		 (-10.0d0 30.0d0 0.0d0)
+		 (0.0d0 0.0d0 40.0d0))))
     (eigenvalues-eigenvectors mat evals evecs ws)
     (values (cl-array evals) (cl-array evecs))))
 
diff --git a/gsll.asd b/gsll.asd
index 11a224ed246e91473a94bf466e351ce3e6d9d79e..ef1a5a00cb571e77bcf3f08c6e94256ead3a7554 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2008-11-29 11:18:52EST gsll.asd>
+;; Time-stamp: <2008-11-30 23:17:52EST gsll.asd>
 ;; $Id$
 
 (asdf:defsystem "gsll"
@@ -94,7 +94,7 @@
    ;; Skip fft for now, I'm not sure how it works in C
    (:file "numerical-integration" :depends-on (init))
    (:module random
-	    :depends-on (init)
+	    :depends-on (init data)
 	    :components
 	    ((:file "rng-types")
 	     (:file "generators" :depends-on (rng-types))
diff --git a/init/generate-examples.lisp b/init/generate-examples.lisp
index d69e4b4fd3eeb7f1ea614ba7949d08e15318c9b0..33223a79ec1866afa6b60151d3a2cc2dcb817189 100644
--- a/init/generate-examples.lisp
+++ b/init/generate-examples.lisp
@@ -1,6 +1,6 @@
 ;; Define examples.
 ;; Liam Healy 2008-09-07 21:00:48EDT generate-tests.lisp
-;; Time-stamp: <2008-11-16 12:13:10EST generate-examples.lisp>
+;; Time-stamp: <2008-11-30 23:02:06EST generate-examples.lisp>
 ;; $Id: $
 
 ;;; Define examples that can be displayed by users with the function
@@ -102,27 +102,28 @@
 ;;;; Generate forms for all array types
 ;;;;****************************************************************************
 
-(defun array-default (spec &optional no-init sync-on-exit)
+(defun array-default (spec &optional no-init)
   "Make an array of the current type and intialize from the pool."
   (declare (special default-element-type starting-element))
   (let ((matrixp (listp spec)))
-    `(,(data-class-name (if matrixp 'matrix 'vector) default-element-type)
-       ,(if no-init
-	    `',spec		   ; no initial values, just dimension
-	    (cons		   ; macro to make initial-contents
-	     'a
-	     (if matrixp
-		 (loop repeat (first spec)
-		    collect
+    (if no-init
+	;; No initial values, just dimension
+	`(make-array* ',default-element-type :dimensions ',spec)
+	;; Initial contents
+	`(make-array*
+	  ',default-element-type
+	  :initial-contents
+	  ',(if matrixp
+		(loop repeat (first spec)
+		   collect
+		   (make-list-from-pool
+		    default-element-type (second spec) starting-element)
+		   do
+		   (incf starting-element (second spec)))
+		(prog1 
 		    (make-list-from-pool
-		     default-element-type (second spec) starting-element)
-		    do
-		    (incf starting-element (second spec)))
-		 (prog1 
-		     (make-list-from-pool
-		      default-element-type spec starting-element)
-		   (incf starting-element spec)))))
-       ,sync-on-exit)))
+		     default-element-type spec starting-element)
+		  (incf starting-element spec)))))))
 
 (defun scalar-default (&optional float-type)
   "Make a scalar of the current type from the pool.  For complex
@@ -162,7 +163,7 @@
        (setf starting-element 0)
        (stupid-code-walk-eval-some
 	test
-	'(array-default vector-default scalar-default)))))
+	'(array-default scalar-default)))))
 
 (defmacro generate-all-array-tests (name element-types test)
   `(save-test ,name ,@(generate-all-array-tests-body element-types test)))
diff --git a/interpolation/spline-example.lisp b/interpolation/spline-example.lisp
index d11b81646f9fe19c4feb18e60a93534fb8d87213..66caddc66b6084fdd87e7fcf39e9ef05c2a07bd8 100644
--- a/interpolation/spline-example.lisp
+++ b/interpolation/spline-example.lisp
@@ -1,6 +1,6 @@
 ;; Example spline
 ;; Liam Healy, Sat Nov 10 2007 - 21:18
-;; Time-stamp: <2008-03-09 19:21:47EDT spline-example.lisp>
+;; Time-stamp: <2008-11-30 23:45:40EST spline-example.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -20,8 +20,8 @@
   (multiple-value-bind (xarr yarr)
       (spline-example-arrays)
     (letm ((acc (acceleration))
-	   (cxarr (vector-double-float xarr))
-	   (cyarr (vector-double-float yarr))
+	   (cxarr (make-array* 'double-float :dimensions xarr))
+	   (cyarr (make-array* 'double-float :dimensions yarr))
 	   (spline (spline *cubic-spline-interpolation* cxarr cyarr)))
       (loop for xi from (aref xarr 0) below (aref xarr 9) by 0.01d0
 	    collect (list xi (evaluate-spline spline xi acc))))))
diff --git a/linear-algebra/exponential.lisp b/linear-algebra/exponential.lisp
index ecd54695ea27b38d208114e9d9849162f2282796..b9d56e5fd3a4e59287d3efa631f5d562611922fb 100644
--- a/linear-algebra/exponential.lisp
+++ b/linear-algebra/exponential.lisp
@@ -1,6 +1,6 @@
 ;; Exponential of a matrix
 ;; Liam Healy 2008-08-10 17:25:35EDT exponential.lisp
-;; Time-stamp: <2008-08-10 18:27:04EDT exponential.lisp>
+;; Time-stamp: <2008-11-30 23:35:44EST exponential.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -25,8 +25,8 @@
 #|
 ;;; A matrix of the form [[0, 1], [-1, 0]]
 ;;; when exponentiated gives [[cos x,  sin x], [-sin x,  cos x]]
-(letm ((mat (matrix-double-float (a (0.0d0 1.0d0) (-1.0d0 0.0d0))))
-       (exp (matrix-double-float '(2 2))))
+(letm ((mat #m((0.0d0 1.0d0) (-1.0d0 0.0d0)))
+       (exp (make-array* 'double-float :dimensions '(2 2))))
   (cl-array (matrix-exponential mat exp)))
 #2A((0.5403023058681384d0 0.841470984807895d0)
     (-0.841470984807895d0 0.5403023058681384d0))
diff --git a/linear-algebra/lu.lisp b/linear-algebra/lu.lisp
index 49bf7e47b2050f09b98234af805eee63028461a9..c2aa8564dc97ab77d1666b9c2a1def2f01756b62 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: <2008-11-29 15:30:45EST lu.lisp>
+;; Time-stamp: <2008-11-30 23:44:35EST lu.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -136,10 +136,10 @@
 (export 'invert-matrix)
 (defun invert-matrix (mat)
   "Invert the matrix."
-  (letm ((mmat (matrix-double-float mat))
+  (letm ((mmat mat)
 	 (dim (array-dimension mat 0))
 	 (per (permutation dim))
-	 (inv (matrix-double-float (list dim dim))))
+	 (inv (make-array* 'double-float :dimensions (list dim dim))))
     (LU-decomposition mmat per)
     (lu-invert mmat per inv)
     (cl-array inv)))
diff --git a/monte-carlo.lisp b/monte-carlo.lisp
index 201432eeef5474f9cf89b8c4c375dd0982c9c734..a55bc64423316e742d2ca511f1a161880e5a68c6 100644
--- a/monte-carlo.lisp
+++ b/monte-carlo.lisp
@@ -1,6 +1,6 @@
 ;; Monte Carlo Integration
 ;; Liam Healy Sat Feb  3 2007 - 17:42
-;; Time-stamp: <2008-10-25 11:44:35EDT monte-carlo.lisp>
+;; Time-stamp: <2008-11-30 23:04:24EST monte-carlo.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -294,22 +294,22 @@
 
 (defun random-walk-plain-example (&optional (nsamples 500000))
   (letm ((ws (monte-carlo-plain 3))
-	 (lower (vector-double-float (a 0.0d0 0.0d0 0.0d0)))
-	 (upper (vector-double-float (a* pi pi pi)))
+	 (lower #m(0.0d0 0.0d0 0.0d0))
+	 (upper (make-array* 'double-float :initial-contents (list pi pi pi)))
 	 (rng (random-number-generator *mt19937* 0)))
     (monte-carlo-integrate-plain monte-carlo-g lower upper nsamples rng ws)))
 
 (defun random-walk-miser-example (&optional (nsamples 500000))
   (letm ((ws (monte-carlo-miser 3))
-	 (lower (vector-double-float (a 0.0d0 0.0d0 0.0d0)))
-	 (upper (vector-double-float (a* pi pi pi)))
+	 (lower #m(0.0d0 0.0d0 0.0d0))
+	 (upper (make-array* 'double-float :initial-contents (list pi pi pi)))
 	 (rng (random-number-generator *mt19937* 0)))
     (monte-carlo-integrate-miser monte-carlo-g lower upper nsamples rng ws)))
 
 (defun random-walk-vegas-example (&optional (nsamples 500000))
   (letm ((ws (monte-carlo-vegas 3))
-	 (lower (vector-double-float (a 0.0d0 0.0d0 0.0d0)))
-	 (upper (vector-double-float (a* pi pi pi)))
+	 (lower #m(0.0d0 0.0d0 0.0d0))
+	 (upper (make-array* 'double-float :initial-contents (list pi pi pi)))
 	 (rng (random-number-generator *mt19937* 0)))
     (monte-carlo-integrate-vegas monte-carlo-g lower upper nsamples rng ws)))
 
diff --git a/polynomial.lisp b/polynomial.lisp
index 97e93abae009eb29c769441914de84569555c6e9..0b4d2c7984535c5d6dd312c368ab8892b7be0b2b 100644
--- a/polynomial.lisp
+++ b/polynomial.lisp
@@ -1,6 +1,6 @@
 ;; Polynomials
 ;; Liam Healy, Tue Mar 21 2006 - 18:33
-;; Time-stamp: <2008-10-25 11:47:01EDT polynomial.lisp>
+;; Time-stamp: <2008-11-30 23:38:39EST polynomial.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -146,11 +146,11 @@
   w of the appropriate size.  The n-1 roots are returned in
   the packed complex array z of length 2(n-1), alternating
   real and imaginary parts."
-  (let ((len (length coefficients)))
-    (letm ((coef (vector-double-float coefficients))
-	   (answer (vector-double-float (* 2 (1- len))))
+  (let ((len (total-size coefficients)))
+    ;; Should this be making a complex array?
+    (letm ((answer (make-array* 'double-float :dimensions (* 2 (1- len))))
 	   (ws (complex-workspace len)))
-      (values-list (polynomial-solve-ws coef ws answer)))))
+      (values-list (polynomial-solve-ws coefficients ws answer)))))
 
 (defmfun polynomial-solve-ws (coefficients workspace answer-pd)
   "gsl_poly_complex_solve"
@@ -172,16 +172,16 @@
 ;;;;****************************************************************************
 
 (save-test polynomial
- (letm ((xa (vector-double-float (a 0.0d0 1.0d0 2.0d0 3.0d0)))
-	(ya (vector-double-float (a 2.5d0 7.2d0 32.7d0 91.0d0)))
-	(dd (vector-double-float 4)))
+ (letm ((xa #m(0.0d0 1.0d0 2.0d0 3.0d0))
+	(ya (#m(2.5d0 7.2d0 32.7d0 91.0d0) ))
+	(dd (make-array* 'double-float :dimensions 4)))
    (divided-difference dd xa ya)
    (list
     (polynomial-eval-divided-difference dd xa 0.0d0)
     (polynomial-eval-divided-difference dd xa 1.0d0)
     (polynomial-eval-divided-difference dd xa 2.0d0)
     (polynomial-eval-divided-difference dd xa 3.0d0)))
- (letm ((vec (vector-double-float (a 1.0d0 2.0d0 3.0d0))))
+ (letm ((vec #m(1.0d0 2.0d0 3.0d0)))
    (polynomial-eval vec -1.0d0))
  (solve-quadratic 1.0d0 0.0d0 1.0d0)
  (solve-quadratic 1.0d0 -2.0d0 1.0d0)
@@ -189,5 +189,5 @@
  (solve-cubic -6.0d0 -13.0d0 42.0d0)
  (solve-cubic-complex -1.0d0 1.0d0 -1.0d0)
  ;; Example from GSL manual
- (polynomial-solve (a double-float -1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
+ (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
 
diff --git a/random/dirichlet.lisp b/random/dirichlet.lisp
index c77e6b547ce36857e100cd73a4df84b70fea5f6a..1f8af5c3099fdcc199eb7ebd9770cf993ad2da78 100644
--- a/random/dirichlet.lisp
+++ b/random/dirichlet.lisp
@@ -1,6 +1,6 @@
 ;; Dirichlet distribution
 ;; Liam Healy, Sun Oct 29 2006
-;; Time-stamp: <2008-10-25 12:13:42EDT dirichlet.lisp>
+;; Time-stamp: <2008-11-30 23:17:06EST dirichlet.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -57,15 +57,15 @@
 ;;; Examples and unit test
 (save-test dirichlet
   (letm ((rng (random-number-generator *mt19937* 0))
-	   (alpha (vector-double-float (a 1.0d0 2.0d0 3.0d0 4.0d0)))
-	   (theta (vector-double-float 4)))
+	   (alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
+	   (theta (make-array* 'double-float :dimensions 4)))
       (dirichlet rng alpha theta)
       (cl-array theta))
-  (letm ((alpha (vector-double-float (a 1.0d0 2.0d0 3.0d0 4.0d0)))
-	  (theta (vector-double-float (a 0.1d0 0.3d0 0.4d0 0.2d0))))
+  (letm ((alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
+	  (theta #m(1.0d0 2.0d0 3.0d0 4.0d0)))
      (dirichlet-pdf alpha theta))
-  (letm ((alpha (vector-double-float (a 1.0d0 2.0d0 3.0d0 4.0d0)))
-	  (theta (vector-double-float (a 0.1d0 0.3d0 0.4d0 0.2d0))))
+  (letm ((alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
+	  (theta #m(1.0d0 2.0d0 3.0d0 4.0d0)))
      (dirichlet-log-pdf alpha theta)))
 
 
diff --git a/random/discrete.lisp b/random/discrete.lisp
index 08602820148d87c8d6a07a14869e693d742d2ff0..d9b6dc94034d4943cb5b1b562f3c545f35a9ea84 100644
--- a/random/discrete.lisp
+++ b/random/discrete.lisp
@@ -1,6 +1,6 @@
 ;; Discrete random variables
 ;; Liam Healy, Sat Nov 11 2006 - 21:51
-;; Time-stamp: <2008-10-25 12:14:54EDT discrete.lisp>
+;; Time-stamp: <2008-11-30 23:13:24EST discrete.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -53,15 +53,12 @@
 (save-test discrete
    ;; Must have two letms because the vector value is not set until
    ;; the body, but the discrete-random needs that set value.
-   (letm ((probabilities (vector-double-float (a 0.25d0 0.5d0 0.25d0))))
+   (letm ((probabilities #m(0.25d0 0.5d0 0.25d0)))
      (letm ((table (discrete-random probabilities))
 	    (rng (random-number-generator *mt19937* 0)))
        (loop for i from 0 to 10
 	     collect
 	     (discrete rng table))))
-   (letm ((probabilities (vector-double-float (a 0.25d0 0.5d0 0.25d0))))
+   (letm ((probabilities #m(0.25d0 0.5d0 0.25d0)))
       (letm ((table (discrete-random probabilities)))
 	(discrete-pdf 1 table))))
-
-
-
diff --git a/random/multinomial.lisp b/random/multinomial.lisp
index f26579dd767cfa1aa6d14cada8c2bb6388d22759..712624cc12da976877aa55ccbd07d28f03d809e9 100644
--- a/random/multinomial.lisp
+++ b/random/multinomial.lisp
@@ -1,6 +1,6 @@
 ;; Multinomial distribution
 ;; Liam Healy, Sat Nov 25 2006 - 16:00
-;; Time-stamp: <2008-10-25 18:04:19EDT multinomial.lisp>
+;; Time-stamp: <2008-11-30 23:19:48EST multinomial.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -50,13 +50,13 @@
 ;;; Examples and unit test
 (save-test multinomial
  (letm ((rng (random-number-generator *mt19937* 0))
-	(p (vector-double-float (a 0.1d0 0.2d0 0.3d0 0.4d0)))
-	(n (vector-signed-byte-32 4)))
+	(p #m(0.1d0 0.2d0 0.3d0 0.4d0))
+	(n (make-array* '(signed-byte 32) :dimensions 4)))
    (multinomial rng 8 p n)
    (cl-array n))
- (letm ((p (vector-double-float (a 0.1d0 0.2d0 0.3d0 0.4d0)))
-	(n (vector-signed-byte-32 (a 5 0 1 2))))
+ (letm ((p #m(0.1d0 0.2d0 0.3d0 0.4d0))
+	(n #31m(5 0 1 2)))
    (multinomial-pdf p N))
- (letm ((p (vector-double-float (a 0.1d0 0.2d0 0.3d0 0.4d0)))
-	(n (vector-signed-byte-32 (a 5 0 1 2))))
+ (letm ((p #m(0.1d0 0.2d0 0.3d0 0.4d0))
+	(n #31m(5 0 1 2)))
    (multinomial-log-pdf p n)))
diff --git a/random/quasi.lisp b/random/quasi.lisp
index 39eef30ea557497fbb3411a9d1e6562a93f0ac76..7144b6f688d4ed5e1e3cfbaf36e6e7751a2f912f 100644
--- a/random/quasi.lisp
+++ b/random/quasi.lisp
@@ -1,6 +1,6 @@
 ;; Quasi-random sequences in arbitrary dimensions.
 ;; Liam Healy, Sun Jul 16 2006 - 15:54
-;; Time-stamp: <2008-10-25 18:09:21EDT quasi.lisp>
+;; Time-stamp: <2008-11-30 23:38:59EST quasi.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -128,7 +128,7 @@
 (save-test quasi-random-number-generators
   ;; This example is given in the GSL documentation
   (letm ((gen (quasi-random-number-generator 2 *sobol*))
-	  (vec (vector-double-float 2)))
+	  (vec (make-array* 'double-float :dimensions 2)))
      (loop repeat 5
 	   do (qrng-get gen vec)
 	   append (coerce (cl-array vec) 'list))))
diff --git a/random/shuffling-sampling.lisp b/random/shuffling-sampling.lisp
index 1fe33b9b2f5f9fbc5c60ff2a09ac1a834c510539..8073d4e7d29eae8c0f7eb94b9c475ca732448e03 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: <2008-10-25 18:16:45EDT shuffling-sampling.lisp>
+;; Time-stamp: <2008-11-30 23:15:32EST shuffling-sampling.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -58,16 +58,16 @@
 ;;; Examples and unit test
 (save-test shuffling-sampling
  (letm ((rng (random-number-generator *mt19937* 0))
-	(v1 (vector-signed-byte-32 (a 1 2 3 4 5 6 7 8))))
+	(v1 #31m(1 2 3 4 5 6 7 8)))
    (shuffle rng v1)
    (cl-array v1))
  (letm ((rng (random-number-generator *mt19937* 0))
-	(v1 (vector-signed-byte-32 (a 1 2 3 4 5 6 7 8)))
-	(v2 (vector-signed-byte-32 4)))
+	(v1 #31m(1 2 3 4 5 6 7 8))
+	(v2 (make-array* '(signed-byte 32) :dimensions 4)))
    (choose-random rng v2 v1)
    (cl-array v2))
  (letm ((rng (random-number-generator *mt19937* 0))
-	(v1 (vector-signed-byte-32 (a 1 2 3 4 5 6 7 8)))
-	(v2 (vector-signed-byte-32 10)))
+	(v1 #31m(1 2 3 4 5 6 7 8))
+	(v2 (make-array* '(signed-byte 32) :dimensions 10)))
    (sample rng v2 v1)
    (cl-array v2)))
diff --git a/series-acceleration.lisp b/series-acceleration.lisp
index 070093808fd838190e67d07506756609b62edc17..bd7dbedd6622bcd1584269c82e2e075f86c69c0e 100644
--- a/series-acceleration.lisp
+++ b/series-acceleration.lisp
@@ -1,6 +1,6 @@
 ;; Series acceleration.
 ;; Liam Healy, Wed Nov 21 2007 - 18:41
-;; Time-stamp: <2008-08-23 11:30:37EDT series-acceleration.lisp>
+;; Time-stamp: <2008-11-30 23:45:16EST series-acceleration.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -110,7 +110,7 @@
 	(sum 0.0d0)
 	(zeta2 (/ (expt pi 2) 6)))
     (letm ((levin (levin maxterms))
-	   (array (vector-double-float maxterms)))
+	   (array (make-array* 'double-float :dimensions maxterms)))
       (dotimes (n maxterms)
 	(setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
 	(incf sum (maref array n)))
diff --git a/solve-minimize-fit/linear-least-squares.lisp b/solve-minimize-fit/linear-least-squares.lisp
index 6c4fdc91c37d39d6498d9aa01febf22a92aca46e..f8e11c19e9977cd6c05fcb549078e9c37818e7d7 100644
--- a/solve-minimize-fit/linear-least-squares.lisp
+++ b/solve-minimize-fit/linear-least-squares.lisp
@@ -1,6 +1,6 @@
 ;; Linear least squares, or linear regression
 ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
-;; Time-stamp: <2008-08-31 12:44:33EDT linear-least-squares.lisp>
+;; Time-stamp: <2008-11-30 23:28:48EST linear-least-squares.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -258,9 +258,9 @@
 (defun univariate-linear-least-squares-example ()
   "First example in Section 36.5 of the GSL manual."
   ;; Results not given in manual so not verified yet.
-  (letm ((x (vector-double-float (a 1970.0d0 1980.0d0 1990.0d0 2000.0d0)))
-	 (y (vector-double-float (a 12.0d0 11.0d0 14.0d0 13.0d0)))
-	 (w (vector-double-float (a 0.1d0 0.2d0 0.3d0 0.4d0))))
+  (letm ((x #m(1970.0d0 1980.0d0 1990.0d0 2000.0d0))
+	 (y #m(12.0d0 11.0d0 14.0d0 13.0d0))
+	 (w #m(0.1d0 0.2d0 0.3d0 0.4d0)))
 	(multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq)
 	    (weighted-linear-fit x w y)
 	  (format t "~&Best fit: Y = ~8,5f + ~8,5f X" c0 c1)
@@ -301,18 +301,18 @@
 (defun mv-linear-least-squares-example (data)
   "Second example in Section 36.5 of the GSL manual."
   (letm ((n (length data)) chisq
-	 (x (matrix-double-float (list n 3)))
-	 (cov (matrix-double-float '(3 3)))
-	 (y (vector-double-float n))
-	 (w (vector-double-float n))
-	 (c (vector-double-float 3)))
+	 (x (make-array* 'double-float :dimensions (list n 3)))
+	 (cov (make-array* 'double-float :dimensions '(3 3)))
+	 (y (make-array* 'double-float :dimensions n))
+	 (w (make-array* 'double-float :dimensions n))
+	 (c (make-array* 'double-float :dimensions 3)))
     (loop for i from 0
-	  for row in data do
-	  (setf (maref X i 0) 1.0d0
-		(maref X i 1) (first row)
-		(maref X i 2) (expt (first row) 2)
-		(maref y i) (second row)
-		(maref w i) (/ (expt (third row) 2))))
+       for row in data do
+       (setf (maref X i 0) 1.0d0
+	     (maref X i 1) (first row)
+	     (maref X i 2) (expt (first row) 2)
+	     (maref y i) (second row)
+	     (maref w i) (/ (expt (third row) 2))))
     (letm ((ws (fit-workspace n 3)))
       (setf chisq
 	    (weighted-linear-mfit X w y c cov ws)))
diff --git a/solve-minimize-fit/minimization-multi.lisp b/solve-minimize-fit/minimization-multi.lisp
index 29f52c50dfd7ae0c996cad42eb73b63de5d2b971..c6d99a5c2db7947b47e2f776f3fcec6209075dc9 100644
--- a/solve-minimize-fit/minimization-multi.lisp
+++ b/solve-minimize-fit/minimization-multi.lisp
@@ -1,6 +1,6 @@
 ;; Multivariate minimization.
 ;; Liam Healy  <Tue Jan  8 2008 - 21:28>
-;; Time-stamp: <2008-08-31 12:25:59EDT minimization-multi.lisp>
+;; Time-stamp: <2008-11-30 23:29:45EST minimization-multi.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -345,7 +345,7 @@
     parabaloid 2 parabaloid-derivative parabaloid-and-derivative)
 
 (defun multimin-example-fletcher-reeves ()
-  (letm ((initial (vector-double-float (a 5.0d0 7.0d0)))
+  (letm ((initial #m(5.0d0 7.0d0))
 	 (minimizer
 	  (mfdfminimizer *conjugate-fletcher-reeves* 2 parabaloid
 			 initial 0.01d0 1.0d-4)))
@@ -375,8 +375,8 @@
 (def-minimization-functions parabaloid-f 2)
 
 (defun multimin-example-nelder-mead ()
-  (letm ((initial (vector-double-float (a 5.0d0 7.0d0)))
-	 (step-size (vector-double-float 2)))
+  (letm ((initial #m(5.0d0 7.0d0))
+	 (step-size (make-array* 'double-float :dimensions 2)))
     (set-all step-size 1.0d0)
     (letm ((minimizer
 	    (mfminimizer *simplex-nelder-mead* 2 parabaloid-f initial step-size)))
diff --git a/solve-minimize-fit/nonlinear-least-squares.lisp b/solve-minimize-fit/nonlinear-least-squares.lisp
index dcab38ab2036c25a192682e1737ed57d8de908c5..f3f0a68842a94b6f452b561c5a27fba1672e555b 100644
--- a/solve-minimize-fit/nonlinear-least-squares.lisp
+++ b/solve-minimize-fit/nonlinear-least-squares.lisp
@@ -1,6 +1,6 @@
 ;; Nonlinear least squares fitting.
 ;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp
-;; Time-stamp: <2008-11-16 14:49:15EST nonlinear-least-squares.lisp>
+;; Time-stamp: <2008-11-30 23:31:40EST nonlinear-least-squares.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -387,13 +387,15 @@
 
 (defun norm-f (fit)
   "Find the norm of the fit function f."
-  (letm ((arr (vector-double-float (cl-array (fdffit-slot fit 'f)))))
-    (euclidean-norm arr)))
+  ;; Fix this
+  (euclidean-norm (fdffit-slot fit 'f)))
 
 (defun solve-nonlinear-least-squares-example ()
-  (letm ((init (vector-double-float (a 1.0d0 0.0d0 0.0d0)))
+  (letm ((init #m(1.0d0 0.0d0 0.0d0))
 	 (covariance
-	  (matrix-double-float (list *number-of-parameters* *number-of-parameters*)))
+	  (make-array* 'double-float
+		       :dimensions
+		       (list *number-of-parameters* *number-of-parameters*)))
 	 (fit (nonlinear-fdffit
 	       *levenberg-marquardt*
 	       *number-of-observations*
diff --git a/solve-minimize-fit/roots-multi.lisp b/solve-minimize-fit/roots-multi.lisp
index e91af7be1c3398947d9075ec8475e1c39426a6f0..0dc853a2e95714cbe11dc1449ac8ea0c79f9883d 100644
--- a/solve-minimize-fit/roots-multi.lisp
+++ b/solve-minimize-fit/roots-multi.lisp
@@ -1,6 +1,6 @@
 ;;; Multivariate roots.                
 ;;; Liam Healy 2008-01-12 12:49:08
-;;; Time-stamp: <2008-11-16 14:44:38EST roots-multi.lisp>
+;;; Time-stamp: <2008-11-30 23:27:04EST roots-multi.lisp>
 ;;; $Id$
 
 (in-package :gsl)
@@ -423,7 +423,7 @@
 (defun roots-multi-example ()
   "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual."
   (let ((max-iter 1000))
-    (letm ((vect (vector-double-float (a -10.0d0 -5.0d0)))
+    (letm ((vect #m(a -10.0d0 -5.0d0))
 	   (solver (mfsolver *hybrid-scaled* rosenbrock vect)))
       (loop for iter from 0
 	 with fnval and argval
@@ -476,7 +476,7 @@
 		   (maref fnval 0)
 		   (maref fnval 1))))
     (let ((max-iter 1000))
-      (letm ((vect (vector-double-float (a -10.0d0 -5.0d0))))
+      (letm ((vect #m(a -10.0d0 -5.0d0)))
 	(letm
 	    ((solver (mfdfsolver *gnewton-mfdfsolver* rosenbrock-f vect)))
 	  (loop for iter from 0
diff --git a/special-functions/bessel.lisp b/special-functions/bessel.lisp
index 391ca0ecc41e096caab8de12542cd210486f04b7..0fa9447c24e1472a9d04a0709b462dbb1b2e5dc7 100644
--- a/special-functions/bessel.lisp
+++ b/special-functions/bessel.lisp
@@ -1,6 +1,6 @@
 ;; Bessel functions
 ;; Liam Healy, Fri Mar 17 2006 - 18:42
-;; Time-stamp: <2008-10-23 22:37:15EDT bessel.lisp>
+;; Time-stamp: <2008-11-30 23:22:00EST bessel.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -460,61 +460,61 @@
 	   (cylindrical-bessel-J0 4.0d0)
 	   (cylindrical-bessel-J1 4.0d0)
 	   (cylindrical-bessel-Jn 2 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (cylindrical-bessel-Jn-array 2.0d0 besarr 2)
 	     (cl-array besarr))
 	   (cylindrical-bessel-Y0 4.0d0)
 	   (cylindrical-bessel-Y1 4.0d0)
 	   (cylindrical-bessel-Yn 3 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (cylindrical-bessel-Yn-array 2.0d0 besarr 2)
 	     (cl-array besarr))
 	   (cylindrical-bessel-I0 4.0d0)
 	   (cylindrical-bessel-I1 4.0d0)
 	   (cylindrical-bessel-In 3 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (cylindrical-bessel-In-array 2.0d0 besarr 2)
 	     (cl-array besarr))
 	   (cylindrical-bessel-I0-scaled 4.0d0)
 	   (cylindrical-bessel-I1-scaled 4.0d0)
 	   (cylindrical-bessel-In-scaled 3 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (cylindrical-bessel-In-scaled-array 2.0d0 besarr 2)
 	     (cl-array besarr))
 	   (cylindrical-bessel-K0 4.0d0)
 	   (cylindrical-bessel-K1 4.0d0)
 	   (cylindrical-bessel-Kn 2 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (cylindrical-bessel-Kn-array 2.0d0 besarr 2)
 	     (cl-array besarr))
 	   (cylindrical-bessel-K0-scaled 4.0d0)
 	   (cylindrical-bessel-K1-scaled 4.0d0)
 	   (cylindrical-bessel-Kn-scaled 2 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (cylindrical-bessel-Kn-array 2.0d0 besarr 2)
 	     (cl-array besarr))
 	   (spherical-bessel-j0 4.0d0)
 	   (spherical-bessel-j1 4.0d0)
 	   (spherical-bessel-j2 4.0d0)
 	   (spherical-bessel-jl 3 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (spherical-bessel-jl-array 4.0d0 besarr)
 	     (cl-array besarr))
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (spherical-bessel-jl-steed-array 4.0d0 besarr)
 	     (cl-array besarr))
 	   (spherical-bessel-y0 4.0d0)
 	   (spherical-bessel-y1 4.0d0)
 	   (spherical-bessel-y2 4.0d0)
 	   (spherical-bessel-yl 2 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (spherical-bessel-yl-array 4.0d0 besarr)
 	     (cl-array besarr))
 	   (spherical-bessel-i0-scaled 4.0d0)
 	   (spherical-bessel-i1-scaled 4.0d0)
 	   (spherical-bessel-i2-scaled 4.0d0)
 	   (spherical-bessel-il-scaled 3 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (spherical-bessel-il-scaled-array 4.0d0 besarr)
 	     (cl-array besarr))
 	   (spherical-bessel-k0-scaled 4.0d0)
@@ -522,11 +522,11 @@
 
 	   (spherical-bessel-k2-scaled 4.0d0)
 	   (spherical-bessel-kl-scaled 3 4.0d0)
-	   (letm ((besarr (vector-double-float 4)))
+	   (letm ((besarr (make-array* 'double-float :dimensions 4)))
 	     (spherical-bessel-kl-scaled-array 4.0d0 besarr)
 	     (cl-array besarr))
 	   (bessel-jnu 3.0d0 4.0d0)
-	   (letm ((besarr (vector-double-float (a 1.0d0 2.0d0 3.0d0))))
+	   (letm ((besarr #m(1.0d0 2.0d0 3.0d0)))
 	     (spherical-Jnu-array 0.5d0 besarr)
 	     (cl-array besarr))
 	   (bessel-Ynu 3.0d0 4.0d0)
diff --git a/special-functions/coulomb.lisp b/special-functions/coulomb.lisp
index 5001b60068aeac214ef906c07fc3cda4061f8d71..1e47356f06c56d20ba23d1851bf2381eb52a0c82 100644
--- a/special-functions/coulomb.lisp
+++ b/special-functions/coulomb.lisp
@@ -1,6 +1,6 @@
 ;; Coulumb functions
 ;; Liam Healy, Sat Mar 18 2006 - 23:23
-;; Time-stamp: <2008-11-16 14:39:08EST coulomb.lisp>
+;; Time-stamp: <2008-11-30 23:44:05EST coulomb.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -128,15 +128,16 @@
   (hydrogenicr-1 1.0d0 2.5d0)
   (hydrogenicr 3 1 1.0d0 2.5d0)
   (coulomb-wave-FG 0.0d0 1.0d0 2.0d0 0)
-  (letm ((arr (vector-double-float 3)))
+  (letm ((arr (make-array* 'double-float :dimensions 3)))
     (coulomb-wave-F-array 0.0d0 1.0d0 2.0d0 arr)
     (cl-array arr))
   (coulomb-wave-fg 1.0d0 2.0d0 2.5d0 1)
-  (letm ((Farr (vector-double-float 3)) (Garr (vector-double-float 3)))
+  (letm ((Farr (make-array* 'double-float :dimensions 3))
+	 (Garr (make-array* 'double-float :dimensions 3)))
     (coulomb-wave-FG-array 1.5d0 1.0d0 1.0d0 Farr Garr)
     (append (coerce (cl-array Farr) 'list) (coerce (cl-array Garr) 'list)))
-  (letm ((arr (vector-double-float 3)))
+  (letm ((arr (make-array* 'double-float :dimensions 3)))
     (coulomb-wave-sphF-array  0.0d0 1.0d0 2.0d0 arr) (cl-array arr))
   (coulomb-cl 1.0d0 2.5d0)
-  (letm ((cl (vector-double-float 3)))
+  (letm ((cl (make-array* 'double-float :dimensions 3)))
     (coulomb-CL-array 0.0d0 1.0d0 cl) (cl-array cl)))
diff --git a/special-functions/gegenbauer.lisp b/special-functions/gegenbauer.lisp
index ee65117c6b19ba66e327c3987d9785aac25abe4b..ec205bbda6c58347c9c6e9e0ca1eefbb82b47ba7 100644
--- a/special-functions/gegenbauer.lisp
+++ b/special-functions/gegenbauer.lisp
@@ -1,6 +1,6 @@
 ;; Gegenbauer polynomials
 ;; Liam Healy, Fri Apr 28 2006 - 20:40
-;; Time-stamp: <2008-10-23 23:09:12EDT gegenbauer.lisp>
+;; Time-stamp: <2008-11-30 23:44:04EST gegenbauer.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -48,5 +48,5 @@
   (gegenbauer-2 1.0d0 3.0d0)
   (gegenbauer-3 1.0d0 3.0d0)
   (gegenbauer 4 1.0d0 3.0d0)
-  (letm ((arr (vector-double-float 4)))
+  (letm ((arr (make-array* 'double-float :dimensions 4)))
       (gegenbauer-array 1.0d0 3.0d0 arr) (cl-array arr)))
diff --git a/special-functions/legendre.lisp b/special-functions/legendre.lisp
index 8cd252e8fdab1bc8ecd8893fcf92deb1b95c01da..81e23b6b317038338d82ff256a10da850c123c0d 100644
--- a/special-functions/legendre.lisp
+++ b/special-functions/legendre.lisp
@@ -1,6 +1,6 @@
 ;; Legendre functions
 ;; Liam Healy, Sat Apr 29 2006 - 19:16
-;; Time-stamp: <2008-11-16 14:39:40EST legendre.lisp>
+;; Time-stamp: <2008-11-30 23:42:24EST legendre.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -258,26 +258,26 @@
   (legendre-Pl -4 0.3d0)
   (legendre-Pl 4 3.0d0)
   (legendre-Pl 4 0.3d0)
-  (letm ((arr (vector-double-float 4)))
+  (letm ((arr (make-array* 'double-float :dimensions 4)))
       (legendre-Pl-array 0.5d0 arr)
       (cl-array arr))
   (legendre-Q0 3.3d0)
   (legendre-Q1 3.3d0)
   (legendre-Ql 2 3.3d0)
   (legendre-Plm 4 3 0.5d0)
-  (letm ((arr (vector-double-float 4)))
+  (letm ((arr (make-array* 'double-float :dimensions 4)))
       (legendre-Plm-array 2 0.5d0 arr)
       (cl-array arr))
-  (letm ((val (vector-double-float 4))
-	   (deriv (vector-double-float 4)))
+  (letm ((val (make-array* 'double-float :dimensions 4))
+	   (deriv (make-array* 'double-float :dimensions 4)))
       (legendre-Plm-deriv-array 2 0.5d0 val deriv)
       (cl-array deriv))
   (legendre-sphplm 1200 1100 0.3d0)
-  (letm ((arr (vector-double-float 4)))
+  (letm ((arr (make-array* 'double-float :dimensions 4)))
       (legendre-sphPlm-array 4 0.5d0 arr)
       (cl-array arr))
-  (letm ((val (vector-double-float 4))
-	   (deriv (vector-double-float 4)))
+  (letm ((val (make-array* 'double-float :dimensions 4))
+	   (deriv (make-array* 'double-float :dimensions 4)))
 	(legendre-sphPlm-deriv-array 4 0.5d0 val deriv)
 	(cl-array deriv))
   (legendre-conicalp-half 3.5d0 10.0d0)
@@ -289,7 +289,7 @@
   (legendre-h3d-0 1.0d0 0.5d0)
   (legendre-h3d-1 1.0d0 0.5d0)
   (legendre-h3d 4 1.0d0 0.5d0)
-  (letm ((arr (vector-double-float 4)))
+  (letm ((arr (make-array* 'double-float :dimensions 4)))
       (legendre-h3d-array 1.0d0 0.5d0 arr)
       (cl-array arr)))
 
diff --git a/statistics/absolute-deviation.lisp b/statistics/absolute-deviation.lisp
index a9d35a28f558c46b8076f57168de9e9eeb0e19de..76922878d7bed99b1fb33d5f2792e2914bfcf745 100644
--- a/statistics/absolute-deviation.lisp
+++ b/statistics/absolute-deviation.lisp
@@ -1,6 +1,6 @@
 ;; Absolute deviation
 ;; Liam Healy, Sun Dec 31 2006 - 13:19
-;; Time-stamp: <2008-10-25 18:46:33EDT absolute-deviation.lisp>
+;; Time-stamp: <2008-11-30 23:35:46EST absolute-deviation.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -49,8 +49,8 @@
 ;;; Examples and unit test
 
 (save-test absolute-deviation
-  (letm ((vec (vector-double-float (a -3.21d0 1.0d0 12.8d0)))
-	 (weights (vector-double-float (a 3.0d0 1.0d0 2.0d0))))
+  (letm ((vec #m(-3.21d0 1.0d0 12.8d0))
+	 (weights #m(3.0d0 1.0d0 2.0d0)))
       (let ((mean (mean vec)))
 	(list
 	 (absolute-deviation vec)
diff --git a/statistics/autocorrelation.lisp b/statistics/autocorrelation.lisp
index 5e8d8ac113f16ab72df5dae103e8b45f748542a6..f4800d578227d2752d83c240dc3006bc4999fd02 100644
--- a/statistics/autocorrelation.lisp
+++ b/statistics/autocorrelation.lisp
@@ -1,6 +1,6 @@
 ;; Autocorrelation
 ;; Liam Healy, Sun Dec 31 2006 - 13:19
-;; Time-stamp: <2008-10-25 18:47:48EDT autocorrelation.lisp>
+;; Time-stamp: <2008-11-30 23:34:37EST autocorrelation.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -26,7 +26,7 @@
 ;;; Examples and unit test
 
 (save-test autocorrelation
-  (letm ((vec (vector-double-float (a -3.21d0 1.0d0 12.8d0))))
+  (letm ((vec #m(-3.21d0 1.0d0 12.8d0)))
       (let ((mean (mean vec)))
 	(list
 	 (autocorrelation vec)
diff --git a/statistics/covariance.lisp b/statistics/covariance.lisp
index 98719f995bc082523f62881ea0685e7e37a1eeb3..fc64e4e12c9bbd47184171709770783ff7827840 100644
--- a/statistics/covariance.lisp
+++ b/statistics/covariance.lisp
@@ -1,6 +1,6 @@
 ;; Covariance
 ;; Liam Healy, Sun Dec 31 2006 - 13:19
-;; Time-stamp: <2008-10-25 18:50:08EDT covariance.lisp>
+;; Time-stamp: <2008-11-30 23:34:37EST covariance.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -29,8 +29,8 @@
 ;;; Examples and unit test
 
 (save-test covariance
-  (letm ((vec1 (vector-double-float (a -3.21d0 1.0d0 12.8d0)))
-	   (vec2 (vector-double-float (a 1.15d0 -1.0d0 0.5d0))))
+  (letm ((vec1 #m(-3.21d0 1.0d0 12.8d0))
+	 (vec2 #m(1.15d0 -1.0d0 0.5d0)))
       (let ((mean1 (mean vec1))
 	    (mean2 (mean vec2)))
 	(list
diff --git a/statistics/higher-moments.lisp b/statistics/higher-moments.lisp
index def1e3f05afea8345b44ddd4c5c22e2a663ed1f1..6abdcead57d84544cb904484158917ec2f229df8 100644
--- a/statistics/higher-moments.lisp
+++ b/statistics/higher-moments.lisp
@@ -1,6 +1,6 @@
 ;; Skewness and kurtosis.
 ;; Liam Healy, Sun Dec 31 2006 - 14:20
-;; Time-stamp: <2008-10-25 18:46:34EDT higher-moments.lisp>
+;; Time-stamp: <2008-11-30 23:32:16EST higher-moments.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -81,7 +81,7 @@
 ;;; Examples and unit test
 
 (save-test higher-moments
-  (letm ((vec (vector-double-float (a -3.21d0 1.0d0 12.8d0))))
+  (letm ((vec #m(-3.21d0 1.0d0 12.8d0)))
       (let* ((mean (mean vec))
 	     (sd (standard-deviation vec mean)))
 	(list
diff --git a/statistics/median-percentile.lisp b/statistics/median-percentile.lisp
index c279f21be3a877260807a3a56ae827a3dd2afcc5..6fd26d02494a5eeb9088051aa6e194d27d2284a8 100644
--- a/statistics/median-percentile.lisp
+++ b/statistics/median-percentile.lisp
@@ -1,6 +1,6 @@
 ;; Median and percentile
 ;; Liam Healy, Sun Dec 31 2006 - 13:19
-;; Time-stamp: <2008-10-25 18:50:07EDT median-percentile.lisp>
+;; Time-stamp: <2008-11-30 23:33:15EST median-percentile.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -53,8 +53,7 @@
 ;;; Examples and unit test
 
 (save-test median-percentile
-  (letm ((vec (vector-double-float (a -3.21d0 1.0d0 12.8d0))))
+  (letm ((vec #m(-3.21d0 1.0d0 12.8d0)))
      (median vec))
-  (letm ((vec (vector-double-float
-	       (a -18.0d0 -12.0d0 -3.21d0 0.5d0 1.0d0 2.7d0 12.8d0))))
+  (letm ((vec #m(-18.0d0 -12.0d0 -3.21d0 0.5d0 1.0d0 2.7d0 12.8d0)))
      (quantile vec 0.75d0)))
diff --git a/wavelet.lisp b/wavelet.lisp
index 98f3373fff67e783ff46ea8ce8f4d9c0f1c95cf4..c5e02ab00a592a03e305cbbf701062205c7e8586 100644
--- a/wavelet.lisp
+++ b/wavelet.lisp
@@ -1,6 +1,6 @@
 ;; Wavelet transforms.
 ;; Liam Healy, Mon Nov 26 2007 - 20:43
-;; Time-stamp: <2008-11-16 13:46:02EST wavelet.lisp>
+;; Time-stamp: <2008-11-30 23:47:59EST wavelet.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -312,7 +312,7 @@
 
 ;;; See GSL manual Section 30.4.
 (defparameter *wavelet-sample*
- (a 0.0462458471760794d0 0.0462458471760794d0 0.0512458471760794d0 0.0712458471760795d0
+ #m(0.0462458471760794d0 0.0462458471760794d0 0.0512458471760794d0 0.0712458471760795d0
  0.0712458471760795d0 0.0662458471760795d0 0.0962458471760795d0 0.101245847176079d0
  0.116245847176079d0 0.121245847176079d0 0.116245847176079d0 0.106245847176079d0
  0.0912458471760794d0 0.101245847176079d0 0.0962458471760795d0 0.0962458471760795d0
@@ -387,11 +387,11 @@
    256) using the 20 largest components of the wavelet transform, while
    setting the others to zero.  See GSL manual Section 30.4."
   (letm ((n (length cl-data))
-	 (vector (vector-double-float cl-data))
+	 (vector cl-data)
 	 (wavelet (wavelet *daubechies-wavelet* 4))
 	 (workspace (wavelet-workspace n)))
     (wavelet-transform-forward wavelet vector 1 workspace)
-    (letm ((absvector (vector-double-float n))
+    (letm ((absvector (make-array* 'double-float :dimensions n))
 	   (permutation (permutation n)))
       (dotimes (i n)
 	(setf (maref absvector i) (abs (maref vector i))))
@@ -407,7 +407,7 @@
 (defun wavelet-forward-example (&optional (cl-data *wavelet-sample*))
   "Simpler example, with only a Daubechies wavelet forward transformation."
   (letm ((n (length cl-data))
-	 (vector (vector-double-float cl-data))
+	 (vector cl-data)
 	 (wavelet (wavelet *daubechies-wavelet* 4))
 	 (workspace (wavelet-workspace n)))
     (wavelet-transform-forward wavelet vector 1 workspace)