Commit 44b9d868 authored by Liam Healy's avatar Liam Healy
Browse files

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.
parent 75f04669
;; 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")
......
;; 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)))
;; 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))
;; 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*))
;; 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*))
......@@ -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)))))))
;; 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*))))