Skip to content
Snippets Groups Projects
Commit a42b3a3a authored by Liam Healy's avatar Liam Healy
Browse files

QRPT: add tests, make arguments optional, fix bugs

Added the translated GSL tests to the linear algebra for QR
decomposition with column pivoting.  More arguments corresponding to
values returned by the functions are now optional.  Fixed the order of
return values from QRPT-decomposition so that permutation is before
signum.  Fixed argument lists to QRPT-QRsolve.
parent bcb8e0b5
No related branches found
No related tags found
No related merge requests found
;; Definition of GSLL system ;; Definition of GSLL system
;; Liam Healy ;; 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" (asdf:defsystem "gsll-tests"
:name "gsll-tests" :name "gsll-tests"
...@@ -137,6 +137,7 @@ ...@@ -137,6 +137,7 @@
(:file "power") (:file "power")
(:file "psi") (:file "psi")
(:file "qr") (:file "qr")
(:file "qrpt")
(:file "quasi-random-number-generators") (:file "quasi-random-number-generators")
(:file "random-number-generators") (:file "random-number-generators")
(:file "rayleigh") (:file "rayleigh")
......
;; Generate matrices used in tests of linear algebra functions ;; Generate matrices used in tests of linear algebra functions
;; Liam Healy 2009-09-19 18:28:31EDT matrix-generation.lisp ;; 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) (in-package :gsl)
...@@ -59,6 +59,9 @@ ...@@ -59,6 +59,9 @@
(defun create-general-matrix (dim0 dim1) (defun create-general-matrix (dim0 dim1)
(create-matrix (lambda (i j) (/ (+ 1 i j))) 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) (defun create-hilbert-matrix (dim)
"Make Hilbert matrix used to test linear algebra functions." "Make Hilbert matrix used to test linear algebra functions."
(create-general-matrix dim dim)) (create-general-matrix dim dim))
...@@ -125,3 +128,5 @@ ...@@ -125,3 +128,5 @@
(defparameter *m35* (create-general-matrix 3 5)) (defparameter *m35* (create-general-matrix 3 5))
(defparameter *m53* (create-general-matrix 5 3)) (defparameter *m53* (create-general-matrix 5 3))
(defparameter *s35* (create-singular-matrix 3 5))
(defparameter *s53* (create-singular-matrix 5 3))
;; QR with column pivoting ;; QR with column pivoting
;; Liam Healy, Fri Apr 28 2006 - 16:53 ;; 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$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -8,20 +8,24 @@ ...@@ -8,20 +8,24 @@
;;; /usr/include/gsl/gsl_linalg.h ;;; /usr/include/gsl/gsl_linalg.h
(defmfun QRPT-decomposition (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" "gsl_linalg_QRPT_decomp"
(((mpointer A) :pointer) ((mpointer tau) :pointer) (((mpointer A) :pointer) ((mpointer tau) :pointer)
((mpointer p) :pointer) ((mpointer permutation) :pointer)
(signum (:pointer :int)) ((mpointer norm) :pointer)) (signum (:pointer :int)) ((mpointer norm) :pointer))
:inputs (A) :inputs (A)
:outputs (A tau p norm) :outputs (A tau permutation norm)
:return (A tau signum p) :return (A tau permutation signum)
:documentation ; FDL :documentation ; FDL
"Factorizes the M-by-N matrix A into "Factorizes the M-by-N matrix A into
the QRP^T decomposition A = Q R P^T. On output the the QRP^T decomposition A = Q R P^T. On output the
diagonal and upper triangular part of the input matrix contain the diagonal and upper triangular part of the input matrix contain the
matrix R. The permutation matrix P is stored in 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 signum. It has the value (-1)^n, where n is the
number of interchanges in the permutation. The vector tau and the number of interchanges in the permutation. The vector tau and the
columns of the lower triangular part of the matrix A contain the columns of the lower triangular part of the matrix A contain the
...@@ -38,37 +42,44 @@ ...@@ -38,37 +42,44 @@
5.4.1).") 5.4.1).")
(defmfun QRPT-decomposition* (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" "gsl_linalg_QRPT_decomp2"
(((mpointer A) :pointer) ((mpointer q) :pointer) (((mpointer A) :pointer) ((mpointer q) :pointer)
((mpointer r) :pointer) ((mpointer tau) :pointer) ((mpointer r) :pointer) ((mpointer tau) :pointer)
((mpointer p) :pointer) (signum (:pointer :int)) ((mpointer permutation) :pointer) (signum (:pointer :int))
((mpointer norm) :pointer)) ((mpointer norm) :pointer))
:inputs (A) :inputs (A)
:outputs (q r norm) :outputs (q r norm)
:return (q r p signum) :return (q r permutation signum)
:documentation ; FDL :documentation ; FDL
"Factorize the matrix A into the decomposition "Factorize the matrix A into the decomposition
A = Q R P^T without modifying A itself and storing the A = Q R P^T without modifying A itself and storing the
output in the separate matrices q and r.") output in the separate matrices q and r.")
(defmfun QRPT-solve (defmfun QRPT-solve
(QR tau p b &optional x-spec (QR tau permutation b &optional x-spec
&aux &aux
(x (if (eq x-spec t) (x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b)) (make-marray 'double-float :dimensions (dimensions b))
x-spec))) x-spec)))
("gsl_linalg_QRPT_svx" "gsl_linalg_QRPT_solve") ("gsl_linalg_QRPT_svx" "gsl_linalg_QRPT_solve")
((((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer p) :pointer) ((((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer b) :pointer)) ((mpointer permutation) :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 x) :pointer))) ((mpointer b) :pointer) ((mpointer x) :pointer)))
:inputs (QR tau p b x) :inputs (QR tau permutation b x)
:outputs (x) :outputs (x)
:return ((or x b)) :return ((or x b))
:documentation ; FDL :documentation ; FDL
"Solve the square system A x = b using the QRP^T decomposition of A "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 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 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 x-spec is a marray, the solution will be returned in it. If x-spec
...@@ -76,23 +87,25 @@ ...@@ -76,23 +87,25 @@
modified. The solution is returned from the function call.") modified. The solution is returned from the function call.")
(defmfun QRPT-QRsolve (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" "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)) ((mpointer b) :pointer) ((mpointer x) :pointer))
:inputs (QR p b) :inputs (Q R permutation b)
:outputs (x) :outputs (x)
:documentation ; FDL :documentation ; FDL
"Solve the square system R P^T x = Q^T b for "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 x. It can be used when the QR decomposition of a matrix is
available in unpacked form as (Q, R).") 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" "gsl_linalg_QRPT_update"
(((mpointer Q) :pointer) ((mpointer R) :pointer) (((mpointer Q) :pointer) ((mpointer R) :pointer)
((mpointer p) :pointer) ((mpointer permutation) :pointer)
((mpointer w) :pointer) ((mpointer v) :pointer)) ((mpointer w) :pointer) ((mpointer v) :pointer))
:inputs (Q R p w v) :inputs (Q R permutation w v)
:outputs (w Q R) :outputs (w Q R)
:return (Q R) :return (Q R)
:documentation ; FDL :documentation ; FDL
...@@ -100,19 +113,20 @@ ...@@ -100,19 +113,20 @@
decomposition (Q, R, p). The update is given by decomposition (Q, R, p). The update is given by
Q'R' = Q R + w v^T where the output matrices Q' and Q'R' = Q R + w v^T where the output matrices Q' and
R' are also orthogonal and right triangular. Note that w is 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 (defmfun QRPT-Rsolve
(QR p b &optional x-spec (QR permutation b &optional x-spec
&aux &aux
(x (if (eq x-spec t) (x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b)) (make-marray 'double-float :dimensions (dimensions b))
x-spec))) x-spec)))
("gsl_linalg_QRPT_Rsvx" "gsl_linalg_QRPT_Rsolve") ("gsl_linalg_QRPT_Rsvx" "gsl_linalg_QRPT_Rsolve")
((((mpointer QR) :pointer) ((mpointer p) :pointer) ((mpointer b) :pointer)) ((((mpointer QR) :pointer) ((mpointer permutation) :pointer)
(((mpointer QR) :pointer) ((mpointer p) :pointer) ((mpointer b) :pointer))
(((mpointer QR) :pointer) ((mpointer permutation) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer))) ((mpointer b) :pointer) ((mpointer x) :pointer)))
:inputs (QR p b x) :inputs (QR permutation b x)
:outputs (x b) :outputs (x b)
:return ((or x b)) :return ((or x b))
:documentation ; FDL :documentation ; FDL
...@@ -124,3 +138,67 @@ ...@@ -124,3 +138,67 @@
a marray, the solution will be 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 non-NIL, on output the solution is stored in x and b is not
modified. The solution is returned from the function call.") 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*))
;; 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*)))))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment