diff --git a/gsll-tests.asd b/gsll-tests.asd
index 7d081ac4e142cd0981596a1113ab7fca149cb20d..82fff6ee9a545a285aef33df5569b1adf9e57108 100644
--- a/gsll-tests.asd
+++ b/gsll-tests.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2009-09-26 12:46:47EDT gsll-tests.asd>
+;; Time-stamp: <2009-09-26 16:44:38EDT gsll-tests.asd>
 
 (asdf:defsystem "gsll-tests"
   :name "gsll-tests"
@@ -137,6 +137,7 @@
 	     (:file "power")
 	     (:file "psi")
 	     (:file "qr")
+	     (:file "qrpt")
 	     (:file "quasi-random-number-generators")
 	     (:file "random-number-generators")
 	     (:file "rayleigh")
diff --git a/linear-algebra/matrix-generation.lisp b/linear-algebra/matrix-generation.lisp
index 0c3e306712914b5a4afa396dbbc6bc238df9bc70..c8f717469c148985a16cd568878bbfe7abcdab96 100644
--- a/linear-algebra/matrix-generation.lisp
+++ b/linear-algebra/matrix-generation.lisp
@@ -1,6 +1,6 @@
 ;; Generate matrices used in tests of linear algebra functions
 ;; Liam Healy 2009-09-19 18:28:31EDT matrix-generation.lisp
-;; Time-stamp: <2009-09-20 23:13:47EDT matrix-generation.lisp>
+;; Time-stamp: <2009-09-26 16:42:11EDT matrix-generation.lisp>
 
 (in-package :gsl)
 
@@ -59,6 +59,9 @@
 (defun create-general-matrix (dim0 dim1)
   (create-matrix (lambda (i j) (/ (+ 1 i j))) dim0 dim1))
 
+(defun create-singular-matrix (dim0 dim1)
+  (create-matrix (lambda (i j) (if (zerop i) 0 (/ (+ 1 i j)))) dim0 dim1))
+
 (defun create-hilbert-matrix (dim)
   "Make Hilbert matrix used to test linear algebra functions."
   (create-general-matrix dim dim))
@@ -125,3 +128,5 @@
 
 (defparameter *m35* (create-general-matrix 3 5))
 (defparameter *m53* (create-general-matrix 5 3))
+(defparameter *s35* (create-singular-matrix 3 5))
+(defparameter *s53* (create-singular-matrix 5 3))
diff --git a/linear-algebra/qrpt.lisp b/linear-algebra/qrpt.lisp
index 4f58c2527ff41b71f0011f11b818e9c7897617db..edada31125ee79c59bd15e68f335b06317664a60 100644
--- a/linear-algebra/qrpt.lisp
+++ b/linear-algebra/qrpt.lisp
@@ -1,6 +1,6 @@
 ;; QR with column pivoting
 ;; Liam Healy, Fri Apr 28 2006 - 16:53
-;; Time-stamp: <2009-04-26 23:05:14EDT qrpt.lisp>
+;; Time-stamp: <2009-09-26 16:46:01EDT qrpt.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,20 +8,24 @@
 ;;; /usr/include/gsl/gsl_linalg.h
 
 (defmfun QRPT-decomposition
-    (A tau p &optional (norm (make-marray 'double-float :dimensions (dim1 A))))
+    (A
+     &optional
+     (tau (make-marray 'double-float :dimensions (min (dim0 A) (dim1 A))))
+     (permutation (make-permutation (dim1 A)))
+     (norm (make-marray 'double-float :dimensions (dim1 A))))
   "gsl_linalg_QRPT_decomp"
   (((mpointer A) :pointer) ((mpointer tau) :pointer)
-   ((mpointer p) :pointer)
+   ((mpointer permutation) :pointer)
    (signum (:pointer :int)) ((mpointer norm) :pointer))
   :inputs (A)
-  :outputs (A tau p norm)
-  :return (A tau signum p)
+  :outputs (A tau permutation norm)
+  :return (A tau permutation signum)
   :documentation			; FDL
   "Factorizes the M-by-N matrix A into
    the QRP^T decomposition A = Q R P^T.  On output the
    diagonal and upper triangular part of the input matrix contain the
    matrix R.  The permutation matrix P is stored in the
-   permutation p.  The sign of the permutation is given by
+   permutation.  The sign of the permutation is given by
    signum. It has the value (-1)^n, where n is the
    number of interchanges in the permutation. The vector tau and the
    columns of the lower triangular part of the matrix A contain the
@@ -38,37 +42,44 @@
    5.4.1).")
 
 (defmfun QRPT-decomposition*
-    (A q r tau p &optional (norm (make-marray 'double-float :dimensions (dim1 A))))
+    (A
+     &optional
+     (q (make-marray 'double-float :dimensions (list (dim0 A) (dim0 A))))
+     (r (make-marray 'double-float :dimensions (dimensions A)))
+     (tau (make-marray 'double-float :dimensions (min (dim0 A) (dim1 A))))
+     (permutation (make-permutation (dim1 A)))
+     (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 :int))
+   ((mpointer permutation) :pointer) (signum (:pointer :int))
    ((mpointer norm) :pointer))
   :inputs (A)
   :outputs (q r norm)
-  :return (q r p signum)
+  :return (q r permutation signum)
   :documentation			; FDL
   "Factorize the matrix A into the decomposition
   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 &optional x-spec
+    (QR tau permutation 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 QR) :pointer) ((mpointer tau) :pointer)
+    ((mpointer permutation) :pointer) ((mpointer b) :pointer))
+   (((mpointer QR) :pointer) ((mpointer tau) :pointer)
+    ((mpointer permutation) :pointer)
     ((mpointer b) :pointer) ((mpointer x) :pointer)))
-  :inputs (QR tau p b x)
+  :inputs (QR tau permutation b x)
   :outputs (x)
   :return ((or x b))
   :documentation			; FDL
   "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
+   into (QR, tau, permutation) 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
@@ -76,23 +87,25 @@
    modified.  The solution is returned from the function call.")
 
 (defmfun QRPT-QRsolve
-    (QR p b &optional (x (make-marray 'double-float :dimensions (dimensions b))))
+    (Q R permutation b
+       &optional (x (make-marray 'double-float :dimensions (dim0 b))))
   "gsl_linalg_QRPT_QRsolve"
-  (((mpointer QR) :pointer) ((mpointer p) :pointer)
+  (((mpointer Q) :pointer) ((mpointer R) :pointer)
+   ((mpointer permutation) :pointer)
    ((mpointer b) :pointer) ((mpointer x) :pointer))
-  :inputs (QR p b)
+  :inputs (Q R permutation b)
   :outputs (x)
   :documentation			; FDL
   "Solve the square system R P^T x = Q^T b for
    x. It can be used when the QR decomposition of a matrix is
    available in unpacked form as (Q, R).")
 
-(defmfun QRPT-update (Q R p w v)
+(defmfun QRPT-update (Q R permutation w v)
   "gsl_linalg_QRPT_update"
   (((mpointer Q) :pointer) ((mpointer R) :pointer)
-   ((mpointer p) :pointer)
+   ((mpointer permutation) :pointer)
    ((mpointer w) :pointer) ((mpointer v) :pointer))
-  :inputs (Q R p w v)
+  :inputs (Q R permutation w v)
   :outputs (w Q R)
   :return (Q R)
   :documentation			; FDL
@@ -100,19 +113,20 @@
    decomposition (Q, R, p). The update is given by
    Q'R' = Q R + w v^T where the output matrices Q' and
    R' are also orthogonal and right triangular. Note that w is
-   destroyed by the update. The permutation p is not changed.")
+   destroyed by the update. The permutation is not changed.")
 
 (defmfun QRPT-Rsolve
-    (QR p b &optional x-spec
+    (QR permutation 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 QR) :pointer) ((mpointer permutation) :pointer)
+    ((mpointer b) :pointer))
+   (((mpointer QR) :pointer) ((mpointer permutation) :pointer)
     ((mpointer b) :pointer) ((mpointer x) :pointer)))
-  :inputs (QR p b x)
+  :inputs (QR permutation b x)
   :outputs (x b)
   :return ((or x b))
   :documentation			; FDL
@@ -124,3 +138,67 @@
   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.")
+
+;;; Examples and unit test, from linalg/test.c
+
+(defun test-qrpt-solve-dim (matrix)
+  "Solve the linear equation using QRPT with the supplied matrix and
+   a right-hand side vector which is the reciprocal of one more than
+   the index."
+  (let ((dim (dim0 matrix)))
+    (multiple-value-bind (QRPT tau permutation)
+	(QRPT-decomposition (copy matrix))
+      (QRPT-solve QRPT tau permutation (create-rhs-vector dim) T))))
+
+(defun test-qrpt-qrsolve-dim (matrix)
+  "Solve the linear equation using QRPT with the supplied matrix and
+   a right-hand side vector which is the reciprocal of one more than
+   the index."
+  (let ((dim (dim0 matrix)))
+    (multiple-value-bind (Q R permutation)
+	(QRPT-decomposition* (copy matrix))
+      (QRPT-QRsolve Q R permutation (create-rhs-vector dim)))))
+
+(defun test-qrpt-decomp-dim (matrix)
+  "Solve the QRPT decomposition with the supplied
+   matrix and a right-hand side vector which is the reciprocal of one
+   more than the index."
+  (multiple-value-bind (QRPT tau permutation)
+      (QRPT-decomposition (copy matrix))
+    (multiple-value-bind (Q R)
+	(QR-unpack QRPT tau)
+      (let* ((qr (matrix-product Q R))
+	     (ans (make-marray 'double-float :dimensions (dimensions qr))))
+	(dotimes (i (dim0 qr) ans)
+	  (setf (row ans i) (permute-inverse permutation (row qr i))))))))
+
+(save-test qrpt
+ ;; test_QRPT_solve
+ (test-qrpt-solve-dim *hilb2*)
+ (test-qrpt-solve-dim *hilb3*)
+ (test-qrpt-solve-dim *hilb4*)
+ (test-qrpt-solve-dim *hilb12*)
+ (test-qrpt-solve-dim *vander2*)
+ (test-qrpt-solve-dim *vander3*)
+ (test-qrpt-solve-dim *vander4*)
+ (test-qrpt-solve-dim *vander12*)
+ ;; test_QRPT_QRsolve
+ (test-qrpt-qrsolve-dim *hilb2*)
+ (test-qrpt-qrsolve-dim *hilb3*)
+ (test-qrpt-qrsolve-dim *hilb4*)
+ (test-qrpt-qrsolve-dim *hilb12*)
+ (test-qrpt-qrsolve-dim *vander2*)
+ (test-qrpt-qrsolve-dim *vander3*)
+ (test-qrpt-qrsolve-dim *vander4*)
+ (test-qrpt-qrsolve-dim *vander12*)
+ ;; test_QRPT_decomp
+ (test-qrpt-decomp-dim *m35*)
+ (test-qrpt-decomp-dim *m53*)
+ (test-qrpt-decomp-dim *hilb2*)
+ (test-qrpt-decomp-dim *hilb3*)
+ (test-qrpt-decomp-dim *hilb4*)
+ (test-qrpt-decomp-dim *hilb12*)
+ (test-qrpt-decomp-dim *vander2*)
+ (test-qrpt-decomp-dim *vander3*)
+ (test-qrpt-decomp-dim *vander4*)
+ (test-qrpt-decomp-dim *vander12*))
diff --git a/tests/qrpt.lisp b/tests/qrpt.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..7be269167e79657958347126c120af10e739d448
--- /dev/null
+++ b/tests/qrpt.lisp
@@ -0,0 +1,123 @@
+;; Regression test QR for GSLL, automatically generated
+
+(in-package :gsl)
+
+;;; Answers inserted from linalg/test.c
+;;; GSL has #define GSL_DBL_EPSILON        2.2204460492503131e-16
+;;; which is 2x what double-float-epsilon is.
+(LISP-UNIT:DEFINE-TEST qrpt
+  ;; QRPT solve
+  (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (list *hilb2-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *HILB2*))))
+  (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *hilb3-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *HILB3*))))
+  (let ((lisp-unit:*epsilon* (* 2 4096 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *hilb4-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *HILB4*))))
+  (let ((lisp-unit:*epsilon* 0.5d0))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *hilb12-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *HILB12*))))
+  (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander2-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *VANDER2*))))
+  (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander3-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *vander3*))))
+  (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander4-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *vander4*))))
+  (let ((lisp-unit:*epsilon* 0.05d0))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander12-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-SOLVE-DIM *vander12*))))
+  ;; QRPT QRsolve
+  (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (list *hilb2-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *HILB2*))))
+  (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (list *hilb3-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *HILB3*))))
+  (let ((lisp-unit:*epsilon* (* 2 4096 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (list *hilb4-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *HILB4*))))
+  (let ((lisp-unit:*epsilon* 0.5d0))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (list *hilb12-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *HILB12*))))
+  (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander2-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *VANDER2*))))
+  (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander3-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *vander3*))))
+  (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander4-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *vander4*))))
+  (let ((lisp-unit:*epsilon* 0.05d0))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander12-soln*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-QRSOLVE-DIM *vander12*))))
+  ;; QRPT decomp
+  (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *m35*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *M35*))))
+  (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *m53*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *M53*))))
+  (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *s35*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *s35*))))
+  (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *s53*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *s53*))))
+  (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *hilb2*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *HILB2*))))
+  (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *hilb3*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *HILB3*))))
+  (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *hilb4*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *HILB4*))))
+  (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *hilb12*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *HILB12*))))
+  (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander2*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *vander2*))))
+  (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander3*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *vander3*))))
+  (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander4*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *vander4*))))
+  (let ((lisp-unit:*epsilon* 0.0005d0))	; "FIXME: bad accuracy"
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST *vander12*)
+     (MULTIPLE-VALUE-LIST (TEST-QRPT-DECOMP-DIM *VANDER12*)))))