Commit 5eb99263 authored by liam's avatar liam

Add tests to earlier special functions. Add eigensystems.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3069 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 83277ac6
;********************************************************
; file: eigensystems.lisp
; description: Eigenvectors and eigenvalues
; date: Sun May 21 2006 - 19:52
; author: Liam M. Healy
; modified: Sun May 21 2006 - 23:22
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; Hermitian matrix functions not defined, await complex matrices.
;;;;****************************************************************************
;;;; Workspace macro
;;;;****************************************************************************
(export 'with-eigensystem-workspace)
(defmacro with-eigensystem-workspace
((ws n &optional vectors hermitian) &body body)
"Create a workspace to use in finding eigenvalues and eigenvectors.
If eigenvectors are desired, set 'vectors to T."
`(let ((,ws ,(if vectors
(if hermitian
`(eigen-hermv-alloc ,n)
`(eigen-symmv-alloc ,n))
(if hermitian
`(eigen-herm-alloc ,n)
`(eigen-symm-alloc ,n)))))
(check-null-pointer ,ws :ENOMEM "workspace allocation failed")
(unwind-protect
(progn ,@body)
,(if vectors
(if hermitian
`(eigen-hermv-free ,n)
`(eigen-symmv-free ,ws))
(if hermitian
`(eigen-herm-free ,n)
`(eigen-symm-free ,ws))))))
;;;;****************************************************************************
;;;; Real Symmetric Matrices
;;;;****************************************************************************
(defun-gsl eigen-symm-alloc ((n :size))
"gsl_eigen_symm_alloc"
:return (:pointer)
:c-return-value :return
:documentation "Allocate a workspace for computing eigenvalues of
@var{n}-by-@var{n} real symmetric matrices. The size of the workspace
is @math{O(2n)}.")
(defun-gsl eigen-symm-free ((w :pointer))
"gsl_eigen_symm_free"
:c-return-value :void
:documentation "Free the memory associated with the workspace @var{w}.")
(defun-gsl eigen-symmv-alloc ((n :size))
"gsl_eigen_symmv_alloc"
:return (:pointer)
:c-return-value :return
:documentation "Allocate a workspace for computing eigenvalues and
eigenvectors of @var{n}-by-@var{n} real symmetric matrices. The size of
the workspace is @math{O(4n)}.")
(defun-gsl eigen-symmv-free ((w :pointer))
"gsl_eigen_symmv_free"
:c-return-value :void
:documentation "Free the memory associated with the workspace @var{w}.")
(defun-gsl eigenvalues-symmetric
((A gsl-matrix-c) (eval gsl-vector-c) (ws :pointer))
"gsl_eigen_symm"
:documentation "Eigenvalues of the real symmetric matrix
@var{A}. Additional workspace of the appropriate size must be provided
in @var{w}. The diagonal and lower triangular part of @var{A} are
destroyed during the computation, but the strict upper triangular part
is not referenced. The eigenvalues are stored in the vector @var{eval}
and are unordered."
:invalidate (A eval)
:return-input (eval))
(defun-gsl eigenvalues-eigenvectors-symmetric
((A gsl-matrix-c) (eval gsl-vector-c) (evec gsl-matrix-c) (ws :pointer))
"gsl_eigen_symmv"
:documentation "The eigenvalues and eigenvectors of the real
symmetric matrix @var{A}. Additional workspace of the appropriate size
must be provided in @var{w}. The diagonal and lower triangular part of
@var{A} are destroyed during the computation, but the strict upper
triangular part is not referenced. The eigenvalues are stored in the
vector @var{eval} and are unordered. The corresponding eigenvectors are
stored in the columns of the matrix @var{evec}. For example, the
eigenvector in the first column corresponds to the first eigenvalue.
The eigenvectors are guaranteed to be mutually orthogonal and normalised
to unit magnitude."
:invalidate (A eval evec)
:return-input (eval evec))
;;;;****************************************************************************
;;;; Complex Hermitian Matrices
;;;;****************************************************************************
(defun-gsl eigen-herm-alloc ((n :size))
"gsl_eigen_herm_alloc"
:return (:pointer)
:c-return-value :return
:documentation "Allocate a workspace for computing eigenvalues of
@var{n}-by-@var{n} complex hermitian matrices. The size of the workspace
is @math{O(3n)}.")
(defun-gsl eigen-herm-free ((w :pointer))
"gsl_eigen_herm_free"
:c-return-value :void
:documentation "Free the memory associated with the workspace @var{w}.")
(defun-gsl eigen-hermv-alloc ((n :size))
"gsl_eigen_hermv_alloc"
:return (:pointer)
:c-return-value :return
:documentation "Allocate a workspace for computing eigenvalues and
eigenvectors of @var{n}-by-@var{n} complex hermitian matrices. The size of
the workspace is @math{O(5n)}.")
(defun-gsl eigen-hermv-free ((w :pointer))
"gsl_eigen_hermv_free"
:c-return-value :void
:documentation "Free the memory associated with the workspace @var{w}.")
;;; The rest awaits complex matrices
;;;;****************************************************************************
;;;; Sorting Eigenvalues and Eigenvectors
;;;;****************************************************************************
(cffi:defcenum eigen-sort-type
"gsl_eigen_sort_t from /usr/include/gsl/gsl_eigen.h."
:value-ascending :value-descending :absolute-ascending :absolute-descending)
(defun-gsl sort-eigenvalues-eigenvectors
((eval gsl-vector-c) (evec gsl-matrix-c) (sort-type eigen-sort-type))
"gsl_eigen_symmv_sort"
:documentation "Simultaneously sort the eigenvalues stored in the vector
@var{eval} and the corresponding real eigenvectors stored in the columns
of the matrix @var{evec} into ascending or descending order according to
the value of the parameter @var{sort_type}, :value-ascending,
:value-descending, :absolute-ascending, :absolute-descending."
:invalidate (eval evec)
:return-input (eval evec))
;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************
#|
(with-data (mat matrix (3 3))
(setf (data mat)
#2A((20.0d0 -10.0d0 0.0d0)
(-10.0d0 30.0d0 0.0d0)
(0.0d0 0.0d0 40.0d0)))
(with-data (evecs matrix (3 3))
(with-data (evals vector 3)
(with-eigensystem-workspace (ws 3 t)
(eigenvalues-eigenvectors-symmetric mat evals evecs ws)
(values (data evals) (data evecs))))))
#(13.819660112501051d0 36.180339887498945d0 40.0d0)
#2A((0.8506508083520399d0 -0.5257311121191336d0 0.0d0)
(0.5257311121191336d0 0.8506508083520399d0 0.0d0)
(0.0d0 0.0d0 1.0d0))
|#
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Sun May 14 2006 - 00:04
; modified: Sun May 21 2006 - 23:22
;********************************************************
;;; $Id: $
......@@ -83,4 +83,5 @@
(:file "svd")
(:file "cholesky")
(:file "diagonal")
(:file "householder")))))
(:file "householder")))
(:file "eigensystems" :depends-on (init data))))
......@@ -3,11 +3,13 @@
; description: Airy functions
; date: Fri Mar 17 2006 - 18:41
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 14:09
;********************************************************
(in-package :gsl)
;;; airy-zero-Ai, airy-zero-Bi seem to only get memory fault
;;;;****************************************************************************
;;;; Airy functions
;;;;****************************************************************************
......@@ -31,7 +33,7 @@
(defun-gsl airy-Ai-scaled ((x :double))
"gsl_sf_airy_Ai_scaled_e"
:return (sf-result)
:documentation "The scaled Airy function @math{S_A(x) Ai(x)}. For @math{x>0} the scaling factor @math{S_A(x)} is @c{$\exp(+(2/3) x^{3/2})$} @math{\exp(+(2/3) x^(3/2))}, and is 1 for @math{x<0}."
:documentation "The scaled Airy function @math{S_A(x) Ai(x)}. For @math{x>0} the scaling factor @math{S_A(x)} is @math{\exp(+(2/3) x^(3/2))}, and is 1 for @math{x<0}."
:mode t)
(defun-gsl airy-Bi-scaled ((x :double))
......@@ -64,12 +66,14 @@
:documentation "The scaled Airy function derivative S_B(x) Bi'(x). For @math{x>0} the scaling factor @math{S_B(x)} is @c{$\exp(-(2/3) x^{3/2})$} @math{exp(-(2/3) x^(3/2))}, and is 1 for @math{x<0}."
:mode t)
(defun-gsl airy-zero-Ai ((s :unsigned-int))
;;; Memory fault
(defun-gsl airy-zero-Ai ((s :size))
"gsl_sf_airy_zero_Ai_e"
:return (sf-result)
:documentation "The location of the @var{s}-th zero of the Airy function @math{Ai(x)}."
:mode t)
;;; Memory fault
(defun-gsl airy-zero-Bi ((s :unsigned-int))
"gsl_sf_airy_zero_Bi_e"
:return (sf-result)
......@@ -88,3 +92,19 @@ function derivative @math{Ai'(x)}."
:return (sf-result)
:documentation "The location of the @var{s}-th zero of the Airy function derivative @math{Bi'(x)}."
:mode t)
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(lisp-unit:define-test airy
(lisp-unit:assert-first-fp-equal "0.157259233805d-01" (airy-ai 2.5d0))
(lisp-unit:assert-first-fp-equal "0.648166073846d+01" (airy-bi 2.5d0))
(lisp-unit:assert-first-fp-equal "0.219322205129d+00" (airy-ai-scaled 2.5d0))
(lisp-unit:assert-first-fp-equal "0.464750480196d+00" (airy-bi-scaled 2.5d0))
(lisp-unit:assert-first-fp-equal "-0.262508810359d-01" (airy-ai-deriv 2.5d0))
(lisp-unit:assert-first-fp-equal "0.942142331733d+01" (airy-bi-deriv 2.5d0))
(lisp-unit:assert-first-fp-equal "-0.366108938475d+00"
(airy-ai-deriv-scaled 2.5d0))
(lisp-unit:assert-first-fp-equal "0.675538444164d+00"
(airy-bi-deriv-scaled 2.5d0)))
......@@ -3,7 +3,7 @@
; description: Clausen function
; date: Sat Mar 18 2006 - 23:18
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 14:43
;********************************************************
;;; $Id:$
......@@ -14,3 +14,6 @@
:documentation
"The Clausen integral @math{Cl_2(x)}."
:return (sf-result))
(lisp-unit:define-test clausen
(lisp-unit:assert-first-fp-equal "0.433598203236d+00" (clausen 2.5d0)))
......@@ -3,7 +3,7 @@
; description: Coulumb functions
; date: Sat Mar 18 2006 - 23:23
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 14:50
;********************************************************
;;; $Id:$
......@@ -115,3 +115,17 @@ functions in the limit @math{\eta \to 0}."
:documentation
"The Coulomb wave function normalization constant @math{C_L(\eta)} for @math{L = Lmin \dots Lmin + kmax}, @math{Lmin > -1}."
:return ((:double (1+ kmax))))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(lisp-unit:define-test coulomb
(lisp-unit:assert-first-fp-equal "0.164169997248d+00"
(hydrogenicr-1 1.0d0 2.5d0))
(lisp-unit:assert-first-fp-equal "0.766646808415d-01"
(hydrogenicr 3 1 1.0d0 2.5d0))
(lisp-unit:assert-first-fp-equal "0.716177996747d-01"
(coulomb-wave-fg 1.0d0 2.0d0 2.5d0 1))
(lisp-unit:assert-first-fp-equal "0.138091464419d-02"
(coulomb-cl 1.0d0 2.5d0)))
......@@ -3,7 +3,7 @@
; description: Coupling coefficients
; date: Sun Mar 19 2006 - 13:30
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 14:53
;********************************************************
;;; $Id: $
......@@ -24,16 +24,10 @@ section are declared in the header file @file{gsl_sf_coupling.h}.
"gsl_sf_coupling_3j_e"
:documentation
"The Wigner 3-j coefficient,
@tex
\beforedisplay
$$
\pmatrix{ja & jb & jc\cr
\pmatrix{ja & jb & jc\cr
ma & mb & mc\cr}
$$
\afterdisplay
@end tex
where the arguments are given in half-integer units, @math{ja} =
@var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc."
where the arguments are given in half-integer units, @math{ja} =
@var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc."
:return (sf-result))
(defun-gsl coupling-6j ((two-ja :int) (two-jb :int) (two-jc :int)
......@@ -41,16 +35,10 @@ where the arguments are given in half-integer units, @math{ja} =
"gsl_sf_coupling_6j_e"
:documentation
"The Wigner 6-j coefficient,
@tex
\beforedisplay
$$
\left\{\matrix{ja & jb & jc\cr
\left\{\matrix{ja & jb & jc\cr
jd & je & jf\cr}\right\}
$$
\afterdisplay
@end tex
where the arguments are given in half-integer units, @math{ja} =
@var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc."
where the arguments are given in half-integer units, @math{ja} =
@var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc."
:return (sf-result))
(defun-gsl coupling-9j ((two-ja :int) (two-jb :int) (two-jc :int)
......@@ -58,18 +46,10 @@ where the arguments are given in half-integer units, @math{ja} =
(two-jg :int) (two-jh :int) (two-ji :int))
"gsl_sf_coupling_9j_e"
:documentation
"The Wigner 9-j coefficient,
@tex
\beforedisplay
$$
\left\{\matrix{ja & jb & jc\cr
"The Wigner 9-j coefficient,
\left\{\matrix{ja & jb & jc\cr
jd & je & jf\cr
jg & jh & ji\cr}\right\}
$$
\afterdisplay
@end tex
where the arguments are given in half-integer units, @math{ja} =
@var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc."
where the arguments are given in half-integer units, @math{ja} =
@var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc."
:return (sf-result))
......@@ -3,7 +3,7 @@
; description: Dawson function
; date: Sun Mar 19 2006 - 14:31
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 14:58
;********************************************************
;;; $Id: $
......@@ -21,3 +21,6 @@ Stegun, Table 7.5. The Dawson functions are declared in the header file
:documentation
"Dawson's integral for @var{x}."
:return (sf-result))
(lisp-unit:define-test dawson
(lisp-unit:assert-first-fp-equal "0.538079506913d+00" (dawson 1.0d0)))
......@@ -3,7 +3,7 @@
; description: Deybe functions
; date: Sun Mar 19 2006 - 14:34
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 14:59
;********************************************************
;;; $Id: $
......@@ -52,3 +52,9 @@ Stegun, Section 27.1.
:documentation
"The fourth-order Debye function @math{D_4(x) = (4/x^4) \int_0^x dt (t^4/(e^t - 1))}."
:return (sf-result))
(lisp-unit:define-test debye
(lisp-unit:assert-first-fp-equal "0.777504634112d+00" (debye-1 1.0d0))
(lisp-unit:assert-first-fp-equal "0.707878475628d+00" (debye-2 1.0d0))
(lisp-unit:assert-first-fp-equal "0.674415564078d+00" (debye-3 1.0d0))
(lisp-unit:assert-first-fp-equal "0.654874068887d+00" (debye-4 1.0d0)))
......@@ -3,7 +3,7 @@
; description: Dilogarithm
; date: Fri Mar 17 2006 - 18:44
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 15:00
;********************************************************
(in-package :gsl)
......@@ -31,3 +31,6 @@
(values
(complex re im)
(complex re-err im-err))))))
(lisp-unit:define-test dilogarithm
(lisp-unit:assert-first-fp-equal "0.164493406685d+01" (dilogarithm 1.0d0)))
......@@ -3,7 +3,7 @@
; description: Elliptic integrals
; date: Mon Mar 20 2006 - 21:50
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 19:06
;********************************************************
;;; $Id: $
......@@ -55,16 +55,10 @@
(defun-gsl elliptic-integral-D ((phi :double) (k :double) (n :double))
"gsl_sf_ellint_D_e"
:documentation
"The incomplete elliptic integral
@math{D(\phi,k,n)} which is defined through the Carlson form @math{RD(x,y,z)}
by the following relation,
@tex
\beforedisplay
$$
D(\phi,k,n) = RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1).
$$
\afterdisplay
@end tex"
"The incomplete elliptic integral @math{D(\phi,k,n)} which is
defined through the Carlson form @math{RD(x,y,z)}
by the following relation,
D(\phi,k,n) = RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1)."
:mode t
:return (sf-result))
......@@ -99,3 +93,29 @@ $$
"The incomplete elliptic integral @math{RJ(x,y,z,p)}."
:mode t
:return (sf-result))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(lisp-unit:define-test elliptic-integrals
(lisp-unit:assert-first-fp-equal "0.157079632679d+01"
(elliptic-integral-K-complete 0.0d0))
(lisp-unit:assert-first-fp-equal "0.157079632679d+01"
(elliptic-integral-E-complete 0.0d0))
(lisp-unit:assert-first-fp-equal "-0.677417538204d+00"
(elliptic-integral-F -0.5d0 2.0d0))
(lisp-unit:assert-first-fp-equal "-0.401819480553d+00"
(elliptic-integral-E -0.5d0 2.0d0))
(lisp-unit:assert-first-fp-equal "-0.617791316339d+00"
(elliptic-integral-P -0.5d0 2.0d0 1.0d0))
(lisp-unit:assert-first-fp-equal "-0.688995144126d-01"
(elliptic-integral-D -0.5d0 2.0d0 1.0d0))
(lisp-unit:assert-first-fp-equal "0.881373587020d+00"
(elliptic-integral-RC 2.0d0 1.0d0))
(lisp-unit:assert-first-fp-equal "0.799259963030d+00"
(elliptic-integral-RD 2.0d0 1.0d0 1.0d0))
(lisp-unit:assert-first-fp-equal "0.881373587020d+00"
(elliptic-integral-RF 2.0d0 1.0d0 1.0d0))
(lisp-unit:assert-first-fp-equal "0.522800417499d+00"
(elliptic-integral-RJ 2.0d0 1.0d0 1.0d0 2.0d0)))
......@@ -3,7 +3,7 @@
; description: Error functions
; date: Mon Mar 20 2006 - 22:31
; author: Liam M. Healy
; modified: Sat Mar 25 2006 - 22:11
; modified: Sun May 21 2006 - 19:08
;********************************************************
;;; $Id: $
......@@ -54,3 +54,11 @@ function
:documentation
"The hazard function for the normal distribution."
:return (sf-result))
(lisp-unit:define-test error-functions
(LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.842700792950d+00" (ERF 1.0d0))
(LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.157299207050d+00" (ERFC 1.0d0))
(LISP-UNIT:ASSERT-FIRST-FP-EQUAL "-0.184960550993d+01" (LOG-ERFC 1.0d0))
(LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.241970724519d+00" (ERF-Z 1.0d0))
(LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.158655253931d+00" (ERF-Q 1.0d0))
(LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.152513527616d+01" (HAZARD 1.0d0)))
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