diff --git a/gsll-tests.asd b/gsll-tests.asd index b3889103aa37b24d1197b61fe72d7673b2d32a18..47f189edc6ba4bc0110fd5a516375ccff8fbaec4 100644 --- a/gsll-tests.asd +++ b/gsll-tests.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2009-08-26 21:04:08EDT gsll-tests.asd> +;; Time-stamp: <2009-09-18 16:26:18EDT gsll-tests.asd> (asdf:defsystem "gsll-tests" :name "gsll-tests" @@ -159,6 +159,7 @@ (:file "sort-vector-smallest-index") (:file "sort-vector-smallest") (:file "spherical-vector") + (:file "svd") (:file "swap-columns") (:file "swap-elements") (:file "swap-row-column") diff --git a/linear-algebra/svd.lisp b/linear-algebra/svd.lisp index 797da70f22036ebee9bfae345733ee1b824870a3..692dbf30f76e707b3c2ffd7d3bccc90a5cdff8ef 100644 --- a/linear-algebra/svd.lisp +++ b/linear-algebra/svd.lisp @@ -1,14 +1,16 @@ ;; Singular Value Decomposition ;; Liam Healy, Tue May 2 2006 - 12:15 -;; Time-stamp: <2009-02-23 20:52:59EST svd.lisp> +;; Time-stamp: <2009-09-18 16:07:28EDT svd.lisp> ;; $Id$ (in-package :gsl) +;;; /usr/include/gsl/gsl_linalg.h + ;;; FDL ;;; A general rectangular M-by-N matrix A has a ;;; singular value decomposition (svd) into the product of an -;;; M-by-N orthogonal matrix U, an N-by- diagonal matrix of +;;; M-by-N orthogonal matrix U, an N-by-N diagonal matrix of ;;; singular values S and the transpose of an ;;; N-by-N orthogonal square matrix V, A = U S V^T ;;; The singular values sigma_i = S_{ii} are all non-negative and are @@ -26,7 +28,10 @@ ;;; tolerance. (defmfun SV-decomposition - (A S V &optional (work (make-marray 'double-float :dimensions (dim1 A)))) + (A &optional + (S (make-marray 'double-float :dimensions (dim1 A))) + (V (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A)))) + (work (make-marray 'double-float :dimensions (dim1 A)))) "gsl_linalg_SV_decomp" (((mpointer A) :pointer) ((mpointer V) :pointer) ((mpointer S) :pointer) ((mpointer work) :pointer)) @@ -46,9 +51,12 @@ This routine uses the Golub-Reinsch SVD algorithm.") (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)))) + (A + &optional + (S (make-marray 'double-float :dimensions (dim1 A))) + (V (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A)))) + (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) @@ -61,7 +69,10 @@ faster for M >> N. It requires the vector work of length N and the N-by-N matrix X as additional working space.") -(defmfun SV-jacobi-decomposition (A S V) +(defmfun SV-jacobi-decomposition + (A &optional + (S (make-marray 'double-float :dimensions (dim1 A))) + (V (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A))))) "gsl_linalg_SV_decomp_jacobi" (((mpointer A) :pointer) ((mpointer V) :pointer) ((mpointer S) :pointer)) @@ -95,3 +106,44 @@ In the over-determined case where A has more rows than columns the system is solved in the least squares sense, returning the solution x which minimizes ||A x - b||_2.") + +;;; Examples and unit test, from linalg/test.c +;;; These are general to all the linear solver techniques, so +;;; more tests need to be made. + +(defun create-hilbert-matrix (dim) + "Make Hilbert matrix used to test linear algebra functions." + (let ((matrix (make-marray 'double-float :dimensions (list dim dim)))) + (dotimes (i dim matrix) + (dotimes (j dim) + (setf (maref matrix i j) (coerce (/ (+ 1 i j)) 'double-float)))))) + +(defun create-vandermonde-matrix (dim) + "Make Van der Monde matrix used to test linear algebra functions." + (let ((matrix (make-marray 'double-float :dimensions (list dim dim)))) + (dotimes (i dim matrix) + (dotimes (j dim) + (setf (maref matrix i j) + (coerce (expt (1+ i) (- dim j 1)) 'double-float)))))) + +(defun test-sv-solve-dim (matrix) + "Solve the linear equation using SVD with the supplied matrix and + a right-hand side vector which is the reciprocal of one more than + the index." + (let* ((dim (dim0 matrix)) + (rhs (make-marray 'double-float :dimensions dim))) + (dotimes (i dim) + (setf (maref rhs i) (coerce (1+ i) 'double-float))) + (multiple-value-bind (u q d) + (SV-decomposition (copy matrix)) + (SV-solve u q d rhs)))) + +(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))) diff --git a/tests/svd.lisp b/tests/svd.lisp new file mode 100644 index 0000000000000000000000000000000000000000..73883f218cd84450ee35ef460db0488973d4d6a9 --- /dev/null +++ b/tests/svd.lisp @@ -0,0 +1,60 @@ +;; Regression test SVD for GSLL, aputomatically 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 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))) + (MULTIPLE-VALUE-LIST + (TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 2))))) + (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))) + (MULTIPLE-VALUE-LIST + (TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 3))))) + (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))))) + (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))))) + (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))))) + (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))))) + (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))))) + (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))))))