From 44b9d868a51985da43fd0c3ea82a7b06d0a66c68 Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Sun, 20 Sep 2009 23:23:18 -0400 Subject: [PATCH] QR: More optional arguments, add tests form GSL Changed required to optional arguments for several of the QR linear algebra functions, for returned quantities. 45 new tests added for QR that are translated from GSL's tests. Matrices and solutions that are used by more than one of the linear algebra techniques are defined in matrix-generation.lisp. --- gsll-tests.asd | 3 +- linear-algebra/lu.lisp | 18 +-- linear-algebra/matrix-generation.lisp | 65 +++++++- linear-algebra/qr.lisp | 142 +++++++++++++++- linear-algebra/svd.lisp | 18 +-- tests/lu.lisp | 59 ++----- tests/qr.lisp | 223 ++++++++++++++++++++++++++ tests/svd.lisp | 50 ++---- 8 files changed, 468 insertions(+), 110 deletions(-) create mode 100644 tests/qr.lisp diff --git a/gsll-tests.asd b/gsll-tests.asd index 47f189ed..b7ea1509 100644 --- a/gsll-tests.asd +++ b/gsll-tests.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2009-09-18 16:26:18EDT gsll-tests.asd> +;; Time-stamp: <2009-09-20 21:52:19EDT gsll-tests.asd> (asdf:defsystem "gsll-tests" :name "gsll-tests" @@ -135,6 +135,7 @@ (:file "polynomial") (:file "power") (:file "psi") + (:file "qr") (:file "quasi-random-number-generators") (:file "random-number-generators") (:file "rayleigh") diff --git a/linear-algebra/lu.lisp b/linear-algebra/lu.lisp index 36d3d1bb..beb503ac 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: <2009-09-20 14:55:07EDT lu.lisp> +;; Time-stamp: <2009-09-20 19:35:33EDT lu.lisp> ;; $Id$ (in-package :gsl) @@ -185,12 +185,12 @@ (LU-refine matrix upper permutation rhs initial-solution))))) (save-test lu - (test-lu-solve-dim (create-hilbert-matrix 2)) - (test-lu-solve-dim (create-hilbert-matrix 3)) - (test-lu-solve-dim (create-hilbert-matrix 4)) - (test-lu-solve-dim (create-hilbert-matrix 12)) - (test-lu-solve-dim (create-vandermonde-matrix 2)) - (test-lu-solve-dim (create-vandermonde-matrix 3)) - (test-lu-solve-dim (create-vandermonde-matrix 4)) - (test-lu-solve-dim (create-vandermonde-matrix 12)) + (test-lu-solve-dim *hilb2*) + (test-lu-solve-dim *hilb3*) + (test-lu-solve-dim *hilb4*) + (test-lu-solve-dim *hilb12*) + (test-lu-solve-dim *vander2*) + (test-lu-solve-dim *vander3*) + (test-lu-solve-dim *vander4*) + (test-lu-solve-dim *vander12*) (test-lu-solve-dim (create-complex-matrix 7))) diff --git a/linear-algebra/matrix-generation.lisp b/linear-algebra/matrix-generation.lisp index 6292de65..0c3e3067 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-19 22:39:15EDT matrix-generation.lisp> +;; Time-stamp: <2009-09-20 23:13:47EDT matrix-generation.lisp> (in-package :gsl) @@ -14,18 +14,21 @@ ;;;; General array creation from indices ;;;;**************************************************************************** -;;; Maybe this should be exported. Come to think of it, didn't Glen +;;; These should be exported. Come to think of it, didn't Glen ;;; have something more general than this? +(defun set-matrix (matrix function) + (dotimes (i (dim0 matrix) matrix) + (dotimes (j (dim1 matrix)) + (setf (maref matrix i j) + (coerce (funcall function i j) (element-type matrix)))))) + (defun create-matrix (function dim0 &optional (dim1 dim0) (element-type 'double-float)) "Make a matrix of the specified dimensions, with contents based on a function of the element indices i, j." - (let ((matrix - (make-marray (cl-single element-type) :dimensions (list dim0 dim1)))) - (dotimes (i dim0 matrix) - (dotimes (j dim1) - (setf (maref matrix i j) - (coerce (funcall function i j) element-type)))))) + (set-matrix + (make-marray (cl-single element-type) :dimensions (list dim0 dim1)) + function)) (defun create-vector (function dim &optional (element-type 'double-float)) @@ -43,6 +46,12 @@ (lambda (i j) (if (= i j) (maref vector i) 0)) (dim0 vector))) +(defun constant-matrix + (constant dim0 &optional (dim1 dim0) (element-type 'double-float)) + (let ((cst (coerce constant element-type))) + (create-matrix (lambda (i j) (declare (ignore i j)) cst) + dim0 dim1 element-type))) + ;;;;**************************************************************************** ;;;; Specific arrays used in linear algebra tests ;;;;**************************************************************************** @@ -76,3 +85,43 @@ (create-vector (lambda (i) (complex (1+ (* 2 i)) (+ 2 (* 2 i)))) 7 element-type) (create-vector '1+ dim element-type))) + +(defparameter *hilb2* (create-hilbert-matrix 2)) +(defparameter *hilb3* (create-hilbert-matrix 3)) +(defparameter *hilb4* (create-hilbert-matrix 4)) +(defparameter *hilb12* (create-hilbert-matrix 12)) + +(defparameter *hilb2-soln* + (make-marray 'double-float :initial-contents '(-8.0d0 18.0d0))) +(defparameter *hilb3-soln* + (make-marray 'double-float :initial-contents '(27.0d0 -192.0d0 210.0d0))) +(defparameter *hilb4-soln* + (make-marray 'double-float + :initial-contents '(-64.0d0 900.0d0 -2520.0d0 1820.0d0))) +(defparameter *hilb12-soln* + (make-marray 'double-float :initial-contents + '(-1728.0d0 245388.0d0 -8528520.0d0 + 127026900.0d0 -1009008000.0d0 4768571808.0d0 + -14202796608.0d0 27336497760.0d0 -33921201600.0d0 + 26189163000.0d0 -11437874448.0d0 2157916488.0d0))) + +(defparameter *vander2* (create-vandermonde-matrix 2)) +(defparameter *vander3* (create-vandermonde-matrix 3)) +(defparameter *vander4* (create-vandermonde-matrix 4)) +(defparameter *vander12* (create-vandermonde-matrix 12)) + +(defparameter *vander2-soln* + (make-marray 'double-float :initial-contents '(1.0d0 0.0d0))) +(defparameter *vander3-soln* + (make-marray 'double-float :initial-contents + '(0.0d0 1.0d0 0.0d0))) +(defparameter *vander4-soln* + (make-marray 'double-float :initial-contents + '(0.0d0 0.0d0 1.0d0 0.0d0))) +(defparameter *vander12-soln* + (make-marray 'double-float :initial-contents + '(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 + 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0))) + +(defparameter *m35* (create-general-matrix 3 5)) +(defparameter *m53* (create-general-matrix 5 3)) diff --git a/linear-algebra/qr.lisp b/linear-algebra/qr.lisp index 268e8b1d..bf463a1f 100644 --- a/linear-algebra/qr.lisp +++ b/linear-algebra/qr.lisp @@ -1,10 +1,12 @@ ;; QR decomposition ;; Liam Healy 2008-02-17 11:05:20EST qr.lisp -;; Time-stamp: <2009-02-23 21:42:14EST qr.lisp> +;; Time-stamp: <2009-09-20 22:26:14EDT qr.lisp> ;; $Id$ (in-package :gsl) +;;; /usr/include/gsl/gsl_linalg.h + ;;; FDL ;;; A general rectangular M-by-N matrix A has a ;;; QR decomposition into the product of an orthogonal @@ -18,7 +20,10 @@ ;;; columns of Q form an orthonormal basis for the range of A, ;;; ran(A), when A has full column rank. -(defmfun QR-decomposition (A tau) +(defmfun QR-decomposition + (A + &optional + (tau (make-marray 'double-float :dimensions (min (dim0 A) (dim1 A))))) "gsl_linalg_QR_decomp" (((mpointer A) :pointer) ((mpointer tau) :pointer)) :inputs (A) @@ -64,8 +69,10 @@ 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)))) + (QR tau b + &optional + (x (make-marray 'double-float :dimensions (dim1 QR))) + (residual (make-marray 'double-float :dimensions (dim0 QR)))) "gsl_linalg_QR_lssolve" (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer) @@ -127,7 +134,11 @@ 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) +(defmfun QR-unpack + (QR tau + &optional + (Q (make-marray 'double-float :dimensions (list (dim0 QR) (dim0 QR)))) + (R (make-marray 'double-float :dimensions (dimensions QR)))) "gsl_linalg_QR_unpack" (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer Q) :pointer) ((mpointer R) :pointer)) @@ -138,7 +149,8 @@ (QR, tau) into the matrices Q and R where Q is M-by-M and R is M-by-N.") -(defmfun QR-QRsolve (Q R b x) +(defmfun QR-QRsolve + (Q R b &optional (x (make-marray 'double-float :dimensions (dim0 b)))) "gsl_linalg_QR_QRsolve" (((mpointer Q) :pointer) ((mpointer R) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)) @@ -184,3 +196,121 @@ 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-qr-solve-dim (matrix) + "Solve the linear equation using QR 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 (QR tau) + (QR-decomposition (copy matrix)) + (QR-solve QR tau (create-rhs-vector dim) T)))) + +(defun test-qr-qrsolve-dim (matrix) + "Solve the linear equation using QR 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 (QR tau) + (QR-decomposition (copy matrix)) + (multiple-value-bind (Q R) + (QR-unpack QR tau) + (QR-QRsolve Q R (create-rhs-vector dim)))))) + +(defun test-qr-lssolve-dim (matrix) + "Solve the linear equation using QR least squares with the supplied + matrix and a right-hand side vector which is the reciprocal of one + more than the index. Returns the solution and the residual." + (let ((dim (dim0 matrix))) + (multiple-value-bind (QR tau) + (QR-decomposition (copy matrix)) + ;; Residual not checked. + (QR-solve-least-squares QR tau (create-rhs-vector dim))))) + +(defun test-qr-decomp-dim (matrix) + "Solve the QR decomposition with the supplied + matrix and a right-hand side vector which is the reciprocal of one + more than the index." + (multiple-value-bind (QR tau) + (QR-decomposition (copy matrix)) + (multiple-value-bind (Q R) + (QR-unpack QR tau) + (matrix-product Q R)))) + +(defun test-qr-update-dim (matrix) + "Test QR rank-1 update; this should return a matrix with all + elements near zero." + (let* ((dim0 (dim0 matrix)) (dim1 (dim1 matrix)) + (u (create-vector (lambda (i) (sin (1+ i))) dim0)) + (v (create-vector + (lambda (i) (+ (cos (+ 2 i)) (sin (+ 3 (expt i 2))))) + dim1)) + (qr1 + (create-matrix + (lambda (i j) (+ (maref matrix i j) (* (maref u i) (maref v j)))) + dim0 dim1)) + (qr2 (copy matrix)) + (w (make-marray 'double-float :dimensions dim0))) + (multiple-value-bind (QR2 tau) + (QR-decomposition qr2) + (multiple-value-bind (Q2 R2) + (QR-unpack QR2 tau) + ;; compute w = Q^T u + (matrix-product Q2 u w 1.0d0 0.0d0 :trans) + (QR-update Q2 R2 w v) + (matrix-product Q2 R2 qr2 1.0d0 0.0d0) + (elt- qr1 qr2))))) + +(save-test qr + ;; test_QR_solve + (test-qr-solve-dim *hilb2*) + (test-qr-solve-dim *hilb3*) + (test-qr-solve-dim *hilb4*) + (test-qr-solve-dim *hilb12*) + (test-qr-solve-dim *vander2*) + (test-qr-solve-dim *vander3*) + (test-qr-solve-dim *vander4*) + (test-qr-solve-dim *vander12*) + ;; test_QR_QRsolve + (test-qr-qrsolve-dim *hilb2*) + (test-qr-qrsolve-dim *hilb3*) + (test-qr-qrsolve-dim *hilb4*) + (test-qr-qrsolve-dim *hilb12*) + (test-qr-qrsolve-dim *vander2*) + (test-qr-qrsolve-dim *vander3*) + (test-qr-qrsolve-dim *vander4*) + (test-qr-qrsolve-dim *vander12*) + ;; test_QR_lssolve + (test-qr-lssolve-dim *m53*) + (test-qr-lssolve-dim *hilb2*) + (test-qr-lssolve-dim *hilb3*) + (test-qr-lssolve-dim *hilb4*) + (test-qr-lssolve-dim *hilb12*) + (test-qr-lssolve-dim *vander2*) + (test-qr-lssolve-dim *vander3*) + (test-qr-lssolve-dim *vander4*) + (test-qr-lssolve-dim *vander12*) + ;; test_QR_decomp + (test-qr-decomp-dim *m35*) + (test-qr-decomp-dim *m53*) + (test-qr-decomp-dim *hilb2*) + (test-qr-decomp-dim *hilb3*) + (test-qr-decomp-dim *hilb4*) + (test-qr-decomp-dim *hilb12*) + (test-qr-decomp-dim *vander2*) + (test-qr-decomp-dim *vander3*) + (test-qr-decomp-dim *vander4*) + (test-qr-decomp-dim *vander12*) + ;; test_QR_update + (test-qr-update-dim *m35*) + (test-qr-update-dim *m53*) + (test-qr-update-dim *hilb2*) + (test-qr-update-dim *hilb3*) + (test-qr-update-dim *hilb4*) + (test-qr-update-dim *hilb12*) + (test-qr-update-dim *vander2*) + (test-qr-update-dim *vander3*) + (test-qr-update-dim *vander4*) + (test-qr-update-dim *vander12*)) diff --git a/linear-algebra/svd.lisp b/linear-algebra/svd.lisp index d2daaf1c..131e8e61 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: <2009-09-19 22:24:15EDT svd.lisp> +;; Time-stamp: <2009-09-20 19:35:07EDT svd.lisp> (in-package :gsl) @@ -118,11 +118,11 @@ (SV-solve u q d (create-rhs-vector dim))))) (save-test svd - (test-sv-solve-dim (create-hilbert-matrix 2)) - (test-sv-solve-dim (create-hilbert-matrix 3)) - (test-sv-solve-dim (create-hilbert-matrix 4)) - (test-sv-solve-dim (create-hilbert-matrix 12)) - (test-sv-solve-dim (create-vandermonde-matrix 2)) - (test-sv-solve-dim (create-vandermonde-matrix 3)) - (test-sv-solve-dim (create-vandermonde-matrix 4)) - (test-sv-solve-dim (create-vandermonde-matrix 12))) + (test-sv-solve-dim *hilb2*) + (test-sv-solve-dim *hilb3*) + (test-sv-solve-dim *hilb4*) + (test-sv-solve-dim *hilb12*) + (test-sv-solve-dim *vander2*) + (test-sv-solve-dim *vander3*) + (test-sv-solve-dim *vander4*) + (test-sv-solve-dim *vander12*)) diff --git a/tests/lu.lisp b/tests/lu.lisp index 09af1c18..5135f6c8 100644 --- a/tests/lu.lisp +++ b/tests/lu.lisp @@ -9,62 +9,36 @@ (LISP-UNIT:DEFINE-TEST LU (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(-8.0d0 18.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-HILBERT-MATRIX 2))))) + (LIST *hilb2-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *hilb2*)))) (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(27.0d0 -192.0d0 210.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-HILBERT-MATRIX 3))))) + (LIST *hilb3-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *hilb3*)))) (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT - :INITIAL-CONTENTS '(-64.0d0 900.0d0 -2520.0d0 1820.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-HILBERT-MATRIX 4))))) + (LIST *hilb4-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *hilb4*)))) (let ((lisp-unit:*epsilon* 0.5d0)) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(-1728.0d0 245388.0d0 -8528520.0d0 - 127026900.0d0 -1009008000.0d0 4768571808.0d0 - -14202796608.0d0 27336497760.0d0 -33921201600.0d0 - 26189163000.0d0 -11437874448.0d0 2157916488.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-HILBERT-MATRIX 12))))) + (LIST *hilb12-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *hilb12*)))) (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 2))))) + (LIST *vander2-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *vander2*)))) (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(0.0d0 1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 3))))) + (LIST *vander3-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *vander3*)))) (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(0.0d0 0.0d0 1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 4))))) + (LIST *vander4-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *vander4*)))) (let ((lisp-unit:*epsilon* 0.05d0)) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 - 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-LU-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 12))))) + (LIST *vander12-soln*) + (MULTIPLE-VALUE-LIST (TEST-LU-SOLVE-DIM *vander12*)))) (let ((lisp-unit:*epsilon* (* 2 1024 1024 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL (LIST @@ -158,4 +132,3 @@ (MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS '(2 2) :INITIAL-CONTENTS '(1.0d0 2.0d0 3.0d0 4.0d0))))))) - diff --git a/tests/qr.lisp b/tests/qr.lisp new file mode 100644 index 00000000..79384ab3 --- /dev/null +++ b/tests/qr.lisp @@ -0,0 +1,223 @@ +;; 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 QR + ;; QR solve + (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (list *hilb2-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *HILB2*)))) + (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *hilb3-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *HILB3*)))) + (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *hilb4-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *HILB4*)))) + (let ((lisp-unit:*epsilon* 0.5d0)) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *hilb12-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *HILB12*)))) + (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander2-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *VANDER2*)))) + (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander3-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *vander3*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander4-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *vander4*)))) + (let ((lisp-unit:*epsilon* 0.05d0)) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander12-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-SOLVE-DIM *vander12*)))) + ;; QR QRsolve + (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (list *hilb2-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *HILB2*)))) + (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (list *hilb3-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *HILB3*)))) + (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (list *hilb4-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *HILB4*)))) + (let ((lisp-unit:*epsilon* 0.5d0)) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (list *hilb12-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *HILB12*)))) + (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander2-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *VANDER2*)))) + (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander3-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *vander3*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander4-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *vander4*)))) + (let ((lisp-unit:*epsilon* 0.05d0)) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander12-soln*) + (MULTIPLE-VALUE-LIST (TEST-QR-QRSOLVE-DIM *vander12*)))) + ;; QR LSsolve + (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(52.5992295702070d0 -337.7263113752073d0 351.8823436427604)) + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(-0.03018843019022285d0 + 0.30523642927546163d0 + -0.4779674081526012d0 + -0.25160239130369944d0 + 0.4989664138130865d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *M53*)))) + (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST + *hilb2-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *HILB2*)))) + (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST + *hilb3-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *HILB3*)))) + (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST + *hilb4-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0 0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *HILB4*)))) + (let ((lisp-unit:*epsilon* 0.5d0)) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST + *hilb12-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 + 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *HILB12*)))) + (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander2-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *VANDER2*)))) + (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander3-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *vander3*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander4-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0 0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *vander4*)))) + (let ((lisp-unit:*epsilon* 0.05d0)) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander12-soln* + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 + 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0))) + (MULTIPLE-VALUE-LIST (TEST-QR-LSSOLVE-DIM *vander12*)))) + ;; QR decomp + (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *m35*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *M35*)))) + (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *m53*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *M53*)))) + (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *hilb2*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *HILB2*)))) + (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *hilb3*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *HILB3*)))) + (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *hilb4*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *HILB4*)))) + (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *hilb12*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *HILB12*)))) + (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander2*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *vander2*)))) + (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander3*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *vander3*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander4*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *vander4*)))) + (let ((lisp-unit:*epsilon* 0.0005d0)) ; "FIXME: bad accuracy" + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST *vander12*) + (MULTIPLE-VALUE-LIST (TEST-QR-DECOMP-DIM *VANDER12*)))) + ;; QR update + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 3 5)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *M35*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 5 3)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *M53*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 2)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *HILB2*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 3)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *HILB3*)))) + (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 4)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *HILB4*)))) + (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 12)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *HILB12*)))) + (let ((lisp-unit:*epsilon* (* 2 8 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 2)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *vander2*)))) + (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 3)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *vander3*)))) + (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 4)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *vander4*)))) + (let ((lisp-unit:*epsilon* 0.0005d0)) ; "FIXME: bad accuracy" + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL + (LIST (constant-matrix 0 12)) + (MULTIPLE-VALUE-LIST (TEST-QR-UPDATE-DIM *VANDER12*))))) diff --git a/tests/svd.lisp b/tests/svd.lisp index 554ed20e..4230897d 100644 --- a/tests/svd.lisp +++ b/tests/svd.lisp @@ -8,53 +8,35 @@ (LISP-UNIT:DEFINE-TEST SVD (let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(-8.0d0 18.0d0))) + (LIST *hilb2-soln*) (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 2))))) + (TEST-SV-SOLVE-DIM *HILB2*)))) (let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(27.0d0 -192.0d0 210.0d0))) + (LIST *hilb3-soln*) (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 3))))) + (TEST-SV-SOLVE-DIM *HILB3*)))) (let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT - :INITIAL-CONTENTS '(-64.0d0 900.0d0 -2520.0d0 1820.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 4))))) + (LIST *hilb4-soln*) + (MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *hilb4*)))) (let ((lisp-unit:*epsilon* 0.5d0)) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(-1728.0d0 245388.0d0 -8528520.0d0 - 127026900.0d0 -1009008000.0d0 4768571808.0d0 - -14202796608.0d0 27336497760.0d0 -33921201600.0d0 - 26189163000.0d0 -11437874448.0d0 2157916488.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 12))))) + (LIST *hilb12-soln*) + (MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *hilb12*)))) (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 2))))) + (LIST *vander2-soln*) + (MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *vander2*)))) (let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(0.0d0 1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 3))))) + (LIST *vander3-soln*) + (MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *vander3*)))) (let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon))) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(0.0d0 0.0d0 1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 4))))) + (LIST *vander4-soln*) + (MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *vander4*)))) (let ((lisp-unit:*epsilon* 0.05d0)) (LISP-UNIT:ASSERT-NUMERICAL-EQUAL - (LIST - (MAKE-MARRAY - 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0))) - (MULTIPLE-VALUE-LIST - (TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 12)))))) + (LIST *vander12-soln*) + (MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *vander12*))))) -- GitLab