Commit bcb8e0b5 authored by Liam Healy's avatar Liam Healy
Browse files

Cholesky: add tests; return argument for matrix-product

Added eight tests from GSL's test for Cholesky decomposition and
solving.  The four decomposition tests incorrectly produce zero in
CCL.  The first optional argument in matrix-product and
matrix-product-symmetric can now be nil, to indicate a matrix of the
right dimensions should be constructed.  This is the default; the
default previously was to make that matrix directly, which meant if
any other optional arguments were supplied, the complicated default
construction for the return matrix would have to be re-entered.
parent c2a88130
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-09-24 22:22:50EDT gsll-tests.asd>
;; Time-stamp: <2009-09-26 12:46:47EDT gsll-tests.asd>
(asdf:defsystem "gsll-tests"
:name "gsll-tests"
......@@ -30,7 +30,7 @@
(:file "cauchy")
(:file "cdot")
(:file "chebyshev")
(:file "chi-squared")
(:file "cholesky")
(:file "clausen")
(:file "column")
(:file "combination")
......
;; BLAS level 2, Matrix-vector operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2009-05-25 14:23:16EDT blas2.lisp>
;; Time-stamp: <2009-09-26 12:52:44EDT blas2.lisp>
;; $Id$
(in-package :gsl)
......@@ -33,17 +33,20 @@
(defmfun matrix-product
((A matrix) (x vector)
&optional
(y (make-marray
element-type :dimensions (matrix-product-dimensions A x)
:initial-element 0))
(alpha 1) (beta 1) (TransA :notrans) TransB)
y
(alpha 1) (beta 1) (TransA :notrans) TransB
&aux
(yarr (or y
(make-marray
element-type :dimensions (matrix-product-dimensions A x)
:initial-element 0))))
("gsl_blas_" :type "gemv")
((transa cblas-transpose) (alpha :element-c-type) ((mpointer A) :pointer)
((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer))
((mpointer x) :pointer) (beta :element-c-type) ((mpointer yarr) :pointer))
:definition :generic
:element-types #+fsbv :float-complex #-fsbv :float
:inputs (A x y)
:outputs (y)
:inputs (A x yarr)
:outputs (yarr)
:documentation ; FDL
"If the second and third arguments are vectors, compute
the matrix-vector product and sum
......@@ -120,17 +123,20 @@
(defmfun matrix-product-symmetric
((A matrix) (x vector)
&optional
(y (make-marray element-type :dimensions (matrix-product-dimensions A x)
:initial-element 0))
(alpha 1) (beta 1) (uplo :upper) (side :left))
&optional y
(alpha 1) (beta 1) (uplo :upper) (side :left)
&aux
(yarr
(or y
(make-marray element-type :dimensions (matrix-product-dimensions A x)
:initial-element 0))))
("gsl_blas_" :type "symv")
((uplo cblas-uplo) (alpha :element-c-type) ((mpointer A) :pointer)
((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer))
((mpointer x) :pointer) (beta :element-c-type) ((mpointer yarr) :pointer))
:definition :generic
:element-types :float
:inputs (A x y)
:outputs (y)
:inputs (A x yarr)
:outputs (yarr)
:documentation ; FDL
"If the second and third arguments are vectors, compute
the matrix-vector product and sum y = alpha A
......
;; BLAS level 3, Matrix-matrix operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2009-05-25 14:26:09EDT blas3.lisp>
;; Time-stamp: <2009-09-26 12:52:25EDT blas3.lisp>
;; $Id$
(in-package :gsl)
......@@ -21,33 +21,38 @@
(defmfun matrix-product
((A matrix) (B matrix)
&optional
(C (make-marray
element-type
:dimensions (matrix-product-dimensions A B) :initial-element 0))
(alpha 1) (beta 1) (TransA :notrans) (TransB :notrans))
C
(alpha 1) (beta 1) (TransA :notrans) (TransB :notrans)
&aux
(Carr
(or C
(make-marray element-type :dimensions (matrix-product-dimensions A B)
:initial-element 0))))
("gsl_blas_" :type "gemm")
((TransA cblas-transpose) (TransB cblas-transpose)
(alpha :element-c-type) ((mpointer A) :pointer)
((mpointer B) :pointer) (beta :element-c-type) ((mpointer C) :pointer))
((mpointer B) :pointer) (beta :element-c-type) ((mpointer Carr) :pointer))
:definition :methods
:element-types #+fsbv :float-complex #-fsbv :float
:inputs (A B C)
:outputs (C))
:inputs (A B Carr)
:outputs (Carr))
(defmfun matrix-product-symmetric
((A matrix) (B matrix)
&optional
(C (make-marray element-type :dimensions (matrix-product-dimensions A B)
:initial-element 0))
(alpha 1) (beta 1) (uplo :upper) (side :left))
&optional C (alpha 1) (beta 1) (uplo :upper) (side :left)
&aux
(Carr
(or C
(make-marray element-type :dimensions (matrix-product-dimensions A B)
:initial-element 0))))
("gsl_blas_" :type "symm")
((side cblas-side) (uplo cblas-uplo) (alpha :element-c-type)
((mpointer A) :pointer) ((mpointer B) :pointer)
(beta :element-c-type) ((mpointer C) :pointer))
(beta :element-c-type) ((mpointer Carr) :pointer))
:definition :methods
:element-types #+fsbv :float-complex #-fsbv :float
:inputs (A B C)
:outputs (C))
:inputs (A B Carr)
:outputs (Carr))
#+fsbv
(defmfun matrix-product-hermitian
......
;; Cholesky Decomposition
;; Liam Healy, Wed May 3 2006 - 16:38
;; Time-stamp: <2009-02-17 21:08:51EST cholesky.lisp>
;; Time-stamp: <2009-09-26 12:36:23EDT cholesky.lisp>
;; $Id$
(in-package :gsl)
......@@ -58,3 +58,31 @@
will replace b. If x-spec is T, then an array will be created and the
solution returned in it. If x-spec is a marray, the solution will
be returned in it.")
;;; Examples and unit test, from linalg/test.c
(defun test-cholesky-solve-dim (matrix)
"Solve the linear equation using Cholesky with the supplied matrix and
a right-hand side vector which is the reciprocal of one more than
the index."
(cholesky-solve
(cholesky-decomposition (copy matrix))
(create-rhs-vector (dim0 matrix))))
(defun test-cholesky-decomp-dim (matrix)
"Decompose using Cholesky and then multiply."
(let ((decomp (cholesky-decomposition (copy matrix))))
(dotimes (row (dim0 matrix) decomp)
(loop for col from (1+ row) below (dim1 matrix) do
(setf (maref decomp row col) 0.0d0)))
(matrix-product decomp decomp nil 1.0d0 0.0d0 :notrans :trans)))
(save-test cholesky
(test-cholesky-solve-dim *hilb2*)
(test-cholesky-solve-dim *hilb3*)
(test-cholesky-solve-dim *hilb4*)
(test-cholesky-solve-dim *hilb12*)
(test-cholesky-decomp-dim *hilb2*)
(test-cholesky-decomp-dim *hilb3*)
(test-cholesky-decomp-dim *hilb4*)
(test-cholesky-decomp-dim *hilb12*))
;; Regression test CHOLESKY for GSLL, automatically generated
(in-package :gsl)
(LISP-UNIT:DEFINE-TEST CHOLESKY
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *hilb2-soln*)
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM *HILB2*))))
(let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *hilb3-soln*)
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM *HILB3*))))
(let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *hilb4-soln*)
(MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *hilb4*))))
(let ((lisp-unit:*epsilon* 0.5d0))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *hilb12-soln*)
(MULTIPLE-VALUE-LIST (TEST-SV-SOLVE-DIM *hilb12*))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *hilb2*)
(MULTIPLE-VALUE-LIST (TEST-CHOLESKY-DECOMP-DIM *HILB2*))))
(let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *HILB3*)
(MULTIPLE-VALUE-LIST (TEST-CHOLESKY-DECOMP-DIM *HILB3*))))
(let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *HILB4*)
(MULTIPLE-VALUE-LIST (TEST-CHOLESKY-DECOMP-DIM *HILB4*))))
(let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST *HILB12*)
(MULTIPLE-VALUE-LIST (TEST-CHOLESKY-DECOMP-DIM *HILB12*)))))
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