Commit 72eaac6d authored by Liam Healy's avatar Liam Healy
Browse files

Matrix generation for linear algebra tests

Matrices (and vectors) may be created from functions of the element
indices.  This is needed to create the matrices used in testing the
linear algebra functions.  It also can be of general use, especially
with further generalization.  For now, there are no exported symbols
from these definitions.
parent 03ad93e0
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-08-26 21:26:40EDT gsll.asd>
;; Time-stamp: <2009-09-19 18:54:29EDT gsll.asd>
(when (asdf:find-system :fsbv nil)
(pushnew :fsbv *features*))
......@@ -116,6 +116,7 @@
((:file "blas1")
(:file "blas2")
(:file "blas3" :depends-on ("blas2"))
(:file "matrix-generation")
(:file "exponential")
(:file "lu")
(:file "qr")
......
;; 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>
(in-package :gsl)
;;; These are general to all the linear solver techniques, so more
;;; tests need to be made. The symbols are not exported because it
;;; assumed they will only be used internally by the test functions.
;;; See linalg/test.c.
;;; 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))
"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))))
(dotimes (i dim0 matrix)
(dotimes (j dim1)
(setf (maref matrix i j)
(coerce (funcall function i j) eltype))))))
(defun create-general-matrix (dim0 dim1)
(create-matrix (lambda (i j) (/ (+ 1 i j))) dim0 dim1))
(defun create-hilbert-matrix (dim)
"Make Hilbert matrix used to test linear algebra functions."
(create-general-matrix dim dim))
(defun create-vandermonde-matrix (dim)
"Make Van der Monde matrix used to test linear algebra functions."
(create-matrix (lambda (i j) (expt (1+ i) (- dim j 1))) dim))
(defun create-moler-matrix (dim)
(create-matrix (lambda (i j) (- (min (1+ i) (1+ j)) 2)) dim))
(defun create-row-matrix (dim0 dim1)
;; 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)))))
;; Singular Value Decomposition
;; Liam Healy, Tue May 2 2006 - 12:15
;; Time-stamp: <2009-09-18 16:07:28EDT svd.lisp>
;; $Id$
;; Time-stamp: <2009-09-19 19:15:31EDT svd.lisp>
(in-package :gsl)
......@@ -108,35 +107,15 @@
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)))
(let ((dim (dim0 matrix)))
(multiple-value-bind (u q d)
(SV-decomposition (copy matrix))
(SV-solve u q d rhs))))
(SV-solve u q d (create-vector dim)))))
(save-test svd
(test-sv-solve-dim (create-hilbert-matrix 2))
......
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