Commit 75f04669 authored by Liam Healy's avatar Liam Healy
Browse files

LU linear algebra tests from GSL

Added nine tests derived from the self-tests that GSL has for LU
decomposition.  All tests except the complex matrix pass within the
tolerances prescribed by GSL.
parent 72eaac6d
;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2009-04-26 23:04:02EDT lu.lisp>
;; Time-stamp: <2009-09-20 14:55:07EDT lu.lisp>
;; $Id$
(in-package :gsl)
......@@ -62,7 +62,9 @@
if it is T, an appropriate vector will be created and the solution
will be computed there. Otherwise it should be a supplied vector.")
(defmfun LU-refine ((A matrix) LU p (b vector) (x vector) residual)
(defmfun LU-refine
((A matrix) LU p (b vector) (x vector)
&optional (residual (make-marray element-type :dimensions (dim0 A))))
("gsl_linalg" :complex "_LU_refine")
(((mpointer A) :pointer) ((mpointer LU) :pointer)
((mpointer p) :pointer)
......@@ -131,8 +133,7 @@
"Compute the sign or phase factor of the determinant of a matrix A,
det(A)/|det(A)|, from its LU decomposition, LU.")
;;; Examples and unit test
;;; Invert a matrix using LU
(export 'invert-matrix)
(defun invert-matrix (mat)
"Invert the matrix."
......@@ -142,6 +143,10 @@
(LU-decomposition mat per)
(lu-invert mat per inv)))
;;; Examples and unit test
;;; These are direct tests of matrix inversion
(save-test
lu
(cl-array
......@@ -164,4 +169,28 @@
(matrix-product-triangular matrix x 1 :upper :notrans :nonunit)
1 :lower :notrans :unit)))))))
;;; From linalg/test.c
(defun test-lu-solve-dim (matrix &optional vector)
"Solve the linear equation using LU 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 (or vector
(create-rhs-vector dim (element-type matrix)))))
(multiple-value-bind (upper permutation signum)
(LU-decomposition (copy matrix))
(declare (ignore signum))
(let ((initial-solution (LU-solve upper rhs permutation T)))
(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 (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 21:51:22EDT matrix-generation.lisp>
;; Time-stamp: <2009-09-19 22:39:15EDT matrix-generation.lisp>
(in-package :gsl)
......@@ -10,18 +10,42 @@
;;; See linalg/test.c.
;;;;****************************************************************************
;;;; General array creation from indices
;;;;****************************************************************************
;;; Maybe this should be exported. Come to think of it, didn't Glen
;;; have something more general than this?
(defun create-matrix
(function dim0 &optional (dim1 dim0) (eltype 'double-float))
(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 eltype) :dimensions (list dim0 dim1))))
(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) eltype))))))
(coerce (funcall function i j) element-type))))))
(defun create-vector
(function dim &optional (element-type 'double-float))
"Make a vector of the specified dimension, with contents
based on a function of the element index."
(let ((vector
(make-marray (cl-single element-type) :dimensions dim)))
(dotimes (i dim vector)
(setf (maref vector i)
(coerce (funcall function i) element-type)))))
(defun create-diagonal-matrix (vector)
"Place the vector along the diagonal of square matrix."
(create-matrix
(lambda (i j) (if (= i j) (maref vector i) 0))
(dim0 vector)))
;;;;****************************************************************************
;;;; Specific arrays used in linear algebra tests
;;;;****************************************************************************
(defun create-general-matrix (dim0 dim1)
(create-matrix (lambda (i j) (/ (+ 1 i j))) dim0 dim1))
......@@ -41,21 +65,14 @@
;; This would be better named a column matrix, but they call it a row.
(create-matrix (lambda (i j) (if (zerop j) (/ (1+ i)) 0)) dim0 dim1))
;;; This should be exported too.
(defun create-diagonal-matrix (vector)
"Place the vector along the diagonal of square matrix."
(create-matrix
(lambda (i j) (if (= i j) (maref vector i) 0))
(dim0 vector)))
(defun create-complex-matrix (dim)
(create-matrix
(lambda (i j)
(complex (/ (+ 1 i j)) (+ 1/2 (expt i 2) (expt j 2))))
dim dim '(complex double-float)))
;;; Create a vector too
(defun create-vector (dim)
(let ((vec (make-marray 'double-float :dimensions dim)))
(dotimes (i dim vec)
(setf (maref vec i) (coerce (1+ i) 'double-float)))))
(defun create-rhs-vector (dim &optional (element-type 'double-float))
(if (subtypep element-type 'complex)
(create-vector
(lambda (i) (complex (1+ (* 2 i)) (+ 2 (* 2 i)))) 7 element-type)
(create-vector '1+ dim element-type)))
;; Singular Value Decomposition
;; Liam Healy, Tue May 2 2006 - 12:15
;; Time-stamp: <2009-09-19 19:15:31EDT svd.lisp>
;; Time-stamp: <2009-09-19 22:24:15EDT svd.lisp>
(in-package :gsl)
......@@ -115,7 +115,7 @@
(let ((dim (dim0 matrix)))
(multiple-value-bind (u q d)
(SV-decomposition (copy matrix))
(SV-solve u q d (create-vector dim)))))
(SV-solve u q d (create-rhs-vector dim)))))
(save-test svd
(test-sv-solve-dim (create-hilbert-matrix 2))
......
;; Regression test LU for GSLL, automatically generated
;; with some manual changes to the results
(in-package :gsl)
;;; Answers (except for invert-matrix)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 LU
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#(-39.65999999999999d0 -49.46000000000001d0
19.679999999999993d0 -5.549999999999997d0))
(MULTIPLE-VALUE-LIST
(LET ((MATRIX
(MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0)
(34.12d0 -6.15d0 49.27d0
-13.49d0)
(32.5d0 42.73d0 -17.24d0
43.31d0)
(-16.12d0 -8.25d0 21.44d0
-49.08d0))))
(VEC
(MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
'(-39.66d0 -49.46d0 19.68d0
-5.55d0))))
(MULTIPLE-VALUE-BIND
(MATRIX PERM)
(LU-DECOMPOSITION MATRIX)
(LET ((X (LU-SOLVE MATRIX VEC PERM)))
(CL-ARRAY
(PERMUTE-INVERSE PERM
(MATRIX-PRODUCT-TRIANGULAR
MATRIX
(MATRIX-PRODUCT-TRIANGULAR
MATRIX X 1 :UPPER :NOTRANS
:NONUNIT)
1 :LOWER :NOTRANS
:UNIT))))))))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(-39.65999999999999d0 -49.46000000000001d0)
#C(19.679999999999996d0 -5.549999999999995d0)
#C(-8.820000000000006d0 25.370000000000005d0)
#C(-30.580000000000002d0 31.67d0)))
(MULTIPLE-VALUE-LIST
(LET ((MATRIX
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0 49.27d0
-13.49d0)
(34.12d0 -6.15d0 49.27d0
-13.49d0 32.5d0 42.73d0
-17.24d0 43.31d0)
(32.5d0 42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0
-49.08d0)
(-16.12d0 -8.25d0 21.44d0
-49.08d0 -39.66d0 -49.46d0
19.68d0 -5.55d0))))
(VEC
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-39.66d0 -49.46d0 19.68d0
-5.55d0 -8.82d0 25.37d0 -30.58d0
31.67d0))))
(MULTIPLE-VALUE-BIND
(MATRIX PERM)
(LU-DECOMPOSITION MATRIX)
(LET ((X (LU-SOLVE MATRIX VEC PERM)))
(CL-ARRAY
(PERMUTE-INVERSE PERM
(MATRIX-PRODUCT-TRIANGULAR
MATRIX
(MATRIX-PRODUCT-TRIANGULAR
MATRIX X 1 :UPPER :NOTRANS
:NONUNIT)
1 :LOWER :NOTRANS
:UNIT))))))))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#2A((-1.9999999999999998d0 1.0d0)
(1.4999999999999998d0 -0.49999999999999994d0)))
(MULTIPLE-VALUE-LIST
(CL-ARRAY
(INVERT-MATRIX
(MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS '(2 2)
:INITIAL-CONTENTS
'(1.0d0 2.0d0 3.0d0 4.0d0)))))))
(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)))))
(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)))))
(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)))))
(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)))))
(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)))))
(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)))))
(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)))))
(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)))))
(let ((lisp-unit:*epsilon* (* 2 1024 1024 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT) :INITIAL-CONTENTS
'(2.40717272023734d+01 -9.84612797621247d+00
-2.69338853034031d+02 8.75455232472528d+01
2.96661356736296d+03 -1.02624473923993d+03
-1.82073812124749d+04 5.67384473042410d+03
5.57693879019068d+04 -1.61540963210502d+04
-7.88941207561151d+04 1.95053812987858d+04
3.95548551241728d+04 -7.76593696255317d+03)))
(MULTIPLE-VALUE-LIST
(TEST-LU-SOLVE-DIM (CREATE-COMPLEX-MATRIX 7)))))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#(-39.65999999999999d0 -49.46000000000001d0
19.679999999999993d0 -5.549999999999997d0))
(MULTIPLE-VALUE-LIST
(LET ((MATRIX
(MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0)
(34.12d0 -6.15d0 49.27d0
-13.49d0)
(32.5d0 42.73d0 -17.24d0
43.31d0)
(-16.12d0 -8.25d0 21.44d0
-49.08d0))))
(VEC
(MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
'(-39.66d0 -49.46d0 19.68d0
-5.55d0))))
(MULTIPLE-VALUE-BIND
(MATRIX PERM)
(LU-DECOMPOSITION MATRIX)
(LET ((X (LU-SOLVE MATRIX VEC PERM)))
(CL-ARRAY
(PERMUTE-INVERSE PERM
(MATRIX-PRODUCT-TRIANGULAR
MATRIX
(MATRIX-PRODUCT-TRIANGULAR
MATRIX X 1 :UPPER :NOTRANS
:NONUNIT)
1 :LOWER :NOTRANS
:UNIT))))))))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#(#C(-39.65999999999999d0 -49.46000000000001d0)
#C(19.679999999999996d0 -5.549999999999995d0)
#C(-8.820000000000006d0 25.370000000000005d0)
#C(-30.580000000000002d0 31.67d0)))
(MULTIPLE-VALUE-LIST
(LET ((MATRIX
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'((-34.5d0 8.24d0 3.29d0 -8.93d0
34.12d0 -6.15d0 49.27d0
-13.49d0)
(34.12d0 -6.15d0 49.27d0
-13.49d0 32.5d0 42.73d0
-17.24d0 43.31d0)
(32.5d0 42.73d0 -17.24d0 43.31d0
-16.12d0 -8.25d0 21.44d0
-49.08d0)
(-16.12d0 -8.25d0 21.44d0
-49.08d0 -39.66d0 -49.46d0
19.68d0 -5.55d0))))
(VEC
(MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
:INITIAL-CONTENTS
'(-39.66d0 -49.46d0 19.68d0
-5.55d0 -8.82d0 25.37d0 -30.58d0
31.67d0))))
(MULTIPLE-VALUE-BIND
(MATRIX PERM)
(LU-DECOMPOSITION MATRIX)
(LET ((X (LU-SOLVE MATRIX VEC PERM)))
(CL-ARRAY
(PERMUTE-INVERSE
PERM
(MATRIX-PRODUCT-TRIANGULAR
MATRIX
(MATRIX-PRODUCT-TRIANGULAR MATRIX X 1 :UPPER :NOTRANS :NONUNIT)
1 :LOWER :NOTRANS :UNIT))))))))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#2A((-1.9999999999999998d0 1.0d0)
(1.4999999999999998d0 -0.49999999999999994d0)))
(MULTIPLE-VALUE-LIST
(CL-ARRAY
(INVERT-MATRIX
(MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS '(2 2)
:INITIAL-CONTENTS
'(1.0d0 2.0d0 3.0d0 4.0d0)))))))
;; Regression test SVD for GSLL, aputomatically generated
;; Regression test SVD for GSLL, automatically generated
(in-package :gsl)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment