diff --git a/init/defmfun-single.lisp b/init/defmfun-single.lisp
index b60b1e1d4f96529670ca4a76c29ce3fddce1052a..b1c2350aa59bd1d689091b2262c604b15ee01259 100644
--- a/init/defmfun-single.lisp
+++ b/init/defmfun-single.lisp
@@ -1,6 +1,6 @@
 ;; Helpers that define a single GSL function interface
 ;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
-;; Time-stamp: <2009-02-18 18:33:05EST defmfun-single.lisp>
+;; Time-stamp: <2009-02-23 22:06:29EST defmfun-single.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -166,7 +166,14 @@
 		   (error "Could not find ~a among the arguments" s)))
 	     allocated))
 	   (clret (or			; better as a symbol macro
-		   (substitute cret-name :c-return return)
+		   (substitute
+		    cret-name :c-return
+		    (mapcar (lambda (sym)
+			      (let ((it (find sym allocated-decl :key 'st-symbol)))
+				(if it
+				    (first (cl-convert-form it))
+				    sym)))
+			    return))
 		   (mapcan #'cl-convert-form
 			   (callback-remove-arg allocated-decl 'st-symbol))
 		   outputs
diff --git a/linear-algebra/householder.lisp b/linear-algebra/householder.lisp
index 1bfd10fd72416407784f32bb6c285f4644947d13..2bd60bc9607151c5ec8e0ea8c963bf7acfd198c7 100644
--- a/linear-algebra/householder.lisp
+++ b/linear-algebra/householder.lisp
@@ -1,6 +1,6 @@
 ;; Householder Transformations
 ;; Liam Healy, Wed May 10 2006 - 10:03
-;; Time-stamp: <2008-12-07 18:31:26EST householder.lisp>
+;; Time-stamp: <2009-02-23 21:09:50EST householder.lisp>
 ;; $Id$
 
 ;;; For householder-transform, it would be nice to be able to pick the
@@ -69,27 +69,25 @@
 ;;;; Householder solver for linear systems
 ;;;;****************************************************************************
 
-(defmfun householder-solve (A b x)
-  "gsl_linalg_HH_solve"
-  (((mpointer A) :pointer) ((mpointer b) :pointer)
-   ((mpointer x) :pointer))
-  :inputs (A b x)
-  :outputs (x A)
-  :return (x)
+(defmfun householder-solve
+    (A b &optional x-spec
+       &aux
+       (x (if (eq x-spec t)
+	      (make-marray 'double-float :dimensions (dimensions b))
+	      x-spec)))
+  ("gsl_linalg_HH_svx" "gsl_linalg_HH_solve")
+  ((((mpointer A) :pointer) ((mpointer b) :pointer))
+   (((mpointer A) :pointer) ((mpointer b) :pointer)
+    ((mpointer x) :pointer)))
+  :inputs (A b)
+  :outputs (x)
+  :return ((or x b))
   :documentation			; FDL
-  "Solve the system A x = b directly using
-   Householder transformations. On output the solution is stored in x
-   and b is not modified. The matrix A is destroyed by the
-   Householder transformations.")
-
-(defmfun householder-svx (A x)
-  "gsl_linalg_HH_svx"
-  (((mpointer A) :pointer) ((mpointer x) :pointer))
-  :inputs (A x)
-  :outputs (x A)
-  :return (x)
-  :documentation			; FDL
-  "Solve the system A x = b in-place using
-  Householder transformations.  On input x should contain the
-  right-hand side b, which is replaced by the solution on output.  The
-  matrix A is destroyed by the Householder transformations.")
+  "Solve the system A x = b directly using Householder
+   transformations. If x-spec is NIL (default), the solution will
+   replace b.  If x-spec is T, then an array will be created and the
+   solution returned in it.  If x-spec is a marray, the solution will
+   be returned in it.  If x-spec is non-NIL, on output the solution is
+   stored in x and b is not modified.  The matrix A is destroyed by
+   the Householder transformations.  The solution is returned from the
+   function call.")
diff --git a/linear-algebra/qr.lisp b/linear-algebra/qr.lisp
index 3a719db17965d60b022480342f3182cafe7075cf..268e8b1de8e7c9e2315bca243d978b4ec259516d 100644
--- a/linear-algebra/qr.lisp
+++ b/linear-algebra/qr.lisp
@@ -1,6 +1,6 @@
 ;; QR decomposition
 ;; Liam Healy 2008-02-17 11:05:20EST qr.lisp
-;; Time-stamp: <2008-12-07 18:32:41EST qr.lisp>
+;; Time-stamp: <2009-02-23 21:42:14EST qr.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -39,31 +39,33 @@
    The algorithm used to perform the decomposition is Householder QR (Golub
    & Van Loan, Matrix Computations, Algorithm 5.2.1).")
 
-(defmfun QR-solve (QR tau b x)
-    "gsl_linalg_QR_solve"
-  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
-   ((mpointer b) :pointer) ((mpointer x) :pointer))
-  :inputs (QR tau b)
-  :outputs (x)
-  :documentation			; FDL
-  "Solve the square system A x = b using the QR
-   decomposition of A into (QR, tau) given by
-   QR-decomp. The least-squares solution for rectangular systems can
-   be found using QR-lssolve.")
-
-(defmfun QR-solvex (QR tau x)
-  "gsl_linalg_QR_svx"
-  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
-   ((mpointer x) :pointer))
-  :inputs (QR tau x)
+(defmfun QR-solve
+    (QR tau b &optional x-spec
+       &aux
+       (x (if (eq x-spec t)
+	      (make-marray 'double-float :dimensions (dimensions b))
+	      x-spec)))
+  ("gsl_linalg_QR_svx" "gsl_linalg_QR_solve")
+  ((((mpointer QR) :pointer) ((mpointer tau) :pointer)
+    ((mpointer b) :pointer))
+   (((mpointer QR) :pointer) ((mpointer tau) :pointer)
+    ((mpointer b) :pointer) ((mpointer x) :pointer)))
+  :inputs (QR tau b x)
   :outputs (x)
+  :return ((or x b))
   :documentation			; FDL
-  "Solves the square system A x = b in-place using the
-  QR decomposition of A into (QR, tau) given by
-  QR-decomp.  On input x should contain the
-  right-hand side b, which is replaced by the solution on output.")
-
-(defmfun QR-solve-least-squares (QR tau b x residual)
+  "Solve the square system A x = b using the QR decomposition of A
+   into (QR, tau) given by QR-decomp. The least-squares solution for
+   rectangular systems can be found using QR-lssolve.  If x-spec is
+   NIL (default), the solution will replace b.  If x-spec is T, then
+   an array will be created and the solution returned in it.  If
+   x-spec is a marray, the solution will be returned in it.  If x-spec
+   is non-NIL, on output the solution is stored in x and b is not
+   modified.  The solution is returned from the function call.")
+
+(defmfun QR-solve-least-squares
+    (QR tau b x
+	&optional (residual (make-marray 'double-float :dimensions (dimensions b))))
   "gsl_linalg_QR_lssolve"
   (((mpointer QR) :pointer) ((mpointer tau) :pointer)
    ((mpointer b) :pointer) ((mpointer x) :pointer)
@@ -104,26 +106,26 @@
    the encoding of the Householder vectors without needing to form the full
    matrix Q.")
 
-(defmfun QR-Rsolve (QR b x)
-  "gsl_linalg_QR_Rsolve"
-  (((mpointer QR) :pointer) ((mpointer b) :pointer)
-   ((mpointer x) :pointer))
-  :inputs (QR b)
-  :outputs (x)
+(defmfun QR-Rsolve
+    (QR b &optional x-spec
+       &aux
+       (x (if (eq x-spec t)
+	      (make-marray 'double-float :dimensions (dimensions b))
+	      x-spec)))
+  ("gsl_linalg_QR_Rsvx" "gsl_linalg_QR_Rsolve")
+  ((((mpointer QR) :pointer) ((mpointer b) :pointer))
+   (((mpointer QR) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)))
+  :inputs (QR b x)
+  :outputs (b x)
+  :return ((or x b))
   :documentation			; FDL
   "Solve the triangular system R x = b for x.  It may be useful if the
-   product b' = Q^T b has already been computed using QR-QTvec}.")
-
-(defmfun QR-Rsolvex (QR x)
-  "gsl_linalg_QR_Rsvx"
-  (((mpointer QR) :pointer) ((mpointer x) :pointer))
-  :inputs (QR x)
-  :outputs (x)
-  :documentation			; FDL
-  "Solve the triangular system R x = b for x in-place. On input x
-  should contain the right-hand side b and is replaced by the solution
-  on output. This function may be useful if the product b' = Q^T b has
-  already been computed using QR-QTvec}.")
+   product b' = Q^T b has already been computed using QR-QTvec.  If
+   x-spec is NIL (default), the solution will replace b.  If x-spec is
+   T, then an array will be created and the solution returned in it.
+   If x-spec is a marray, the solution will be returned in it.  If
+   x-spec is non-NIL, on output the solution is stored in x and b is
+   not modified.  The solution is returned from the function call.")
 
 (defmfun QR-unpack (QR tau Q R)
   "gsl_linalg_QR_unpack"
@@ -161,21 +163,24 @@
   orthogonal and right triangular. Note that w is destroyed by the
   update.")
 
-(defmfun R-solve (R b x)
-  "gsl_linalg_R_solve"
-  (((mpointer R) :pointer) ((mpointer b) :pointer)
-   ((mpointer x) :pointer))
-  :inputs (R b)
-  :outputs (x)
-  :documentation			; FDL
-  "Solves the triangular system R x = b for the N-by-N matrix R.")
-
-(defmfun R-solvex (R x)
-  "gsl_linalg_R_svx"
-  (((mpointer R) :pointer) ((mpointer x) :pointer))
-  :inputs (R x)
+(defmfun R-solve
+    (R b &optional x-spec
+       &aux
+       (x (if (eq x-spec t)
+	      (make-marray 'double-float :dimensions (dimensions b))
+	      x-spec)))
+  ("gsl_linalg_R_svx" "gsl_linalg_R_solve")
+  ((((mpointer R) :pointer) ((mpointer b) :pointer))
+   (((mpointer R) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)))
+  :inputs (R b x)
   :outputs (x)
+  :return ((or x b))
   :documentation			; FDL
-  "Solve the triangular system R x = b in-place. On
-  input x should contain the right-hand side b, which is
-  replaced by the solution on output.")
+  "Solve the triangular system R x = b in-place. On input x should
+  contain the right-hand side b, which is replaced by the solution on
+  output.  If x-spec is NIL (default), the solution will replace b.
+  If x-spec is T, then an array will be created and the solution
+  returned in it.  If x-spec is a marray, the solution will be
+  returned in it.  If x-spec is non-NIL, on output the solution is
+  stored in x and b is not modified.  The solution is returned from
+  the function call.")
diff --git a/linear-algebra/qrpt.lisp b/linear-algebra/qrpt.lisp
index 2f006466d536bbab6167edf94e9edbe7ffe74377..b2bf31a570ba2d2aa15b311936d77ccd3e91f021 100644
--- a/linear-algebra/qrpt.lisp
+++ b/linear-algebra/qrpt.lisp
@@ -1,15 +1,16 @@
 ;; QR with column pivoting
 ;; Liam Healy, Fri Apr 28 2006 - 16:53
-;; Time-stamp: <2008-12-07 18:33:01EST qrpt.lisp>
+;; Time-stamp: <2009-02-23 22:11:10EST qrpt.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
-(defmfun QRPT-decomposition (A tau p signum norm)
+(defmfun QRPT-decomposition
+    (A tau p &optional (norm (make-marray 'double-float :dimensions (dim1 A))))
   "gsl_linalg_QRPT_decomp"
   (((mpointer A) :pointer) ((mpointer tau) :pointer)
    ((mpointer p) :pointer)
-   (signum :pointer) ((mpointer norm) :pointer))
+   (signum :int) ((mpointer norm) :pointer))
   :inputs (A)
   :outputs (A tau p norm)
   :return (A tau signum p)
@@ -34,11 +35,12 @@
    column pivoting (Golub & Van Loan, Matrix Computations, Algorithm
    5.4.1).")
 
-(defmfun QRPT-decomposition* (A q r tau p signum norm)
+(defmfun QRPT-decomposition*
+    (A q r tau p &optional (norm (make-marray 'double-float :dimensions (dim1 A))))
   "gsl_linalg_QRPT_decomp2"
   (((mpointer A) :pointer) ((mpointer q) :pointer)
    ((mpointer r) :pointer) ((mpointer tau) :pointer)
-   ((mpointer p) :pointer) (signum :pointer)
+   ((mpointer p) :pointer) (signum :int)
    ((mpointer norm) :pointer))
   :inputs (A)
   :outputs (q r norm)
@@ -48,29 +50,31 @@
   A = Q R P^T without modifying A itself and storing the
   output in the separate matrices q and r.")
 
-(defmfun QRPT-solve (QR tau p b x)
-  "gsl_linalg_QRPT_solve"
-  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
-   ((mpointer p) :pointer)
-   ((mpointer b) :pointer) ((mpointer x) :pointer))
-  :inputs (QR tau p b)
-  :outputs (x)
-  :documentation 			; FDL
-  "Solve the square system A x = b using the QRP^T
-   decomposition of A into (QR, tau, p) given by #'QRPT-decomposition.")
-
-(defmfun QRPT-solvex (QR tau p x)
-  "gsl_linalg_QRPT_svx"
-  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
-   ((mpointer p) :pointer) ((mpointer x) :pointer))
-  :inputs (QR tau p x)
+(defmfun QRPT-solve
+    (QR tau p b &optional x-spec
+       &aux
+       (x (if (eq x-spec t)
+	      (make-marray 'double-float :dimensions (dimensions b))
+	      x-spec)))
+  ("gsl_linalg_QRPT_svx" "gsl_linalg_QRPT_solve")
+  ((((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer p) :pointer)
+    ((mpointer b) :pointer))
+   (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer p) :pointer)
+    ((mpointer b) :pointer) ((mpointer x) :pointer)))
+  :inputs (QR tau p b x)
   :outputs (x)
+  :return ((or x b))
   :documentation			; FDL
-  "Solve the square system A x = b in-place using the
-   QRP^T decomposition of A into (QR, tau, p). On input x should contain the
-   right-hand side b, which is replaced by the solution on output.")
+  "Solve the square system A x = b using the QRP^T decomposition of A
+   into (QR, tau, p) given by #'QRPT-decomposition.  If x-spec is
+   NIL (default), the solution will replace b.  If x-spec is T, then
+   an array will be created and the solution returned in it.  If
+   x-spec is a marray, the solution will be returned in it.  If x-spec
+   is non-NIL, on output the solution is stored in x and b is not
+   modified.  The solution is returned from the function call.")
 
-(defmfun QRPT-QRsolve (QR p b x)
+(defmfun QRPT-QRsolve
+    (QR p b &optional (x (make-marray 'double-float :dimensions (dimensions b))))
   "gsl_linalg_QRPT_QRsolve"
   (((mpointer QR) :pointer) ((mpointer p) :pointer)
    ((mpointer b) :pointer) ((mpointer x) :pointer))
@@ -96,24 +100,25 @@
    R' are also orthogonal and right triangular. Note that w is
    destroyed by the update. The permutation p is not changed.")
 
-(defmfun QRPT-Rsolve (QR p b x)
-  "gsl_linalg_QRPT_Rsolve"
-  (((mpointer QR) :pointer) ((mpointer p) :pointer)
-   ((mpointer b) :pointer) ((mpointer x) :pointer))
-  :inputs (QR p b)
-  :outputs (x)
-  :documentation			; FDL
-  "Solve the triangular system R P^T x = b for the
-   N-by-N matrix R contained in QR.")
-
-(defmfun QRPT-Rsolvex (QR p x)
-  "gsl_linalg_QRPT_Rsvx"
-  (((mpointer QR) :pointer) ((mpointer p) :pointer)
-   ((mpointer x) :pointer))
-  :inputs (QR p x)
-  :outputs (x)
+(defmfun QRPT-Rsolve
+    (QR p b &optional x-spec
+       &aux
+       (x (if (eq x-spec t)
+	      (make-marray 'double-float :dimensions (dimensions b))
+	      x-spec)))
+  ("gsl_linalg_QRPT_Rsvx" "gsl_linalg_QRPT_Rsolve")
+  ((((mpointer QR) :pointer) ((mpointer p) :pointer) ((mpointer b) :pointer))
+   (((mpointer QR) :pointer) ((mpointer p) :pointer)
+    ((mpointer b) :pointer) ((mpointer x) :pointer)))
+  :inputs (QR p b x)
+  :outputs (x b)
+  :return ((or x b))
   :documentation			; FDL
-  "Solve the triangular system R P^T x = b in-place
-  for the N-by-N matrix R contained in QR. On
-  input x should contain the right-hand side b, which is
-  replaced by the solution on output.")
+  "Solve the triangular system R P^T x = b in-place for the N-by-N
+  matrix R contained in QR. On input x should contain the right-hand
+  side b, which is replaced by the solution on output.  If x-spec is
+  NIL (default), the solution will replace b.  If x-spec is T, then an
+  array will be created and the solution returned in it.  If x-spec is
+  a marray, the solution will be returned in it.  If x-spec is
+  non-NIL, on output the solution is stored in x and b is not
+  modified.  The solution is returned from the function call.")
diff --git a/linear-algebra/svd.lisp b/linear-algebra/svd.lisp
index 6799e363a657f2b41b82c701e26190583456f025..797da70f22036ebee9bfae345733ee1b824870a3 100644
--- a/linear-algebra/svd.lisp
+++ b/linear-algebra/svd.lisp
@@ -1,6 +1,6 @@
 ;; Singular Value Decomposition
 ;; Liam Healy, Tue May  2 2006 - 12:15
-;; Time-stamp: <2008-12-07 18:33:23EST svd.lisp>
+;; Time-stamp: <2009-02-23 20:52:59EST svd.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -25,12 +25,14 @@
 ;;; precision. Small singular values should be edited by choosing a suitable
 ;;; tolerance.
 
-(defmfun SV-decomposition (A S V work)
+(defmfun SV-decomposition
+    (A S V &optional (work (make-marray 'double-float :dimensions (dim1 A))))
   "gsl_linalg_SV_decomp"
   (((mpointer A) :pointer) ((mpointer V) :pointer)
    ((mpointer S) :pointer) ((mpointer work) :pointer))
   :inputs (A)
   :outputs (A S V)
+  :return (A S V)
   :documentation			; FDL
   "Factorize the M-by-N matrix A into
   the singular value decomposition A = U S V^T for M >= N.
@@ -43,14 +45,18 @@
   transpose of V.  A workspace of length N is required in work.
   This routine uses the Golub-Reinsch SVD algorithm.")
 
-(defmfun SV-modified-decomposition (A S V X work)
+(defmfun SV-modified-decomposition
+    (A S V
+       &optional (X (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A))))
+       (work (make-marray 'double-float :dimensions (dim1 A))))
   "gsl_linalg_SV_decomp_mod"
   (((mpointer A) :pointer) ((mpointer X) :pointer)
    ((mpointer V) :pointer)
    ((mpointer S) :pointer) ((mpointer work) :pointer))
   :inputs (A)
   :outputs (A S V)
-  :documentation			; FDL
+  :return (S V)			    ; is A modified by this algorithm?
+  :documentation		    ; FDL
   "The SVD using the modified Golub-Reinsch algorithm, which is
    faster for M >> N.  It requires the vector work of length N and the
    N-by-N matrix X as additional working space.")
@@ -61,22 +67,25 @@
    ((mpointer S) :pointer))
   :inputs (A)
   :outputs (A S V)
-  :documentation			; FDL
+  :return (S V)			    ; is A modified by this algorithm?
+  :documentation		    ; FDL
   "The SVD of the M-by-N matrix A using one-sided Jacobi
    orthogonalization for M >= N.  The Jacobi method can compute singular
    values to higher relative accuracy than Golub-Reinsch algorithms (see
    references for details).")
 
-(defmfun SV-solve (U S V b x)
+(defmfun SV-solve
+    (U S V b &optional (x (make-marray 'double-float :dimensions (dim0 U))))
   "gsl_linalg_SV_solve"
   (((mpointer U) :pointer) ((mpointer V) :pointer)
    ((mpointer S) :pointer)
    ((mpointer b) :pointer) ((mpointer x) :pointer))
   :inputs (U S V b)
   :outputs (x)
+  :return (x)
   :documentation			; FDL
   "Solve the system A x = b using the singular value
-   decomposition (U, S, V) of A given by #'SV-decomp.
+   decomposition (U, S, V) of A given by #'SV-decomposition.
 
    Only non-zero singular values are used in computing the solution. The
    parts of the solution corresponding to singular values of zero are