Commit 12c2be7e authored by Liam Healy's avatar Liam Healy

Non-symmetric eigensystems, array outputs for special functions

Added some functions introduced in GSL 1.9 to solve eigenvectors and
eigenvalues for non-symmetric matrices.  These are untested as yet.
Changed array output for special functions to be optional arguments
which can take a size (integer) and create a marray of the right
size.
parent 0b69f21e
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2009-01-08 21:49:15EST eigensystems.lisp>
;; Time-stamp: <2009-01-11 10:35:06EST eigensystems.lisp>
;; $Id$
(in-package :gsl)
;;; /usr/include/gsl/gsl_eigen.h
;;; Should symmetric matrices form a subclass of matrices, so that
;;; both eigenvalues and eigenvalues-nonsymm could be methods of the
;;; same function?
;;;;****************************************************************************
;;;; Workspace
;;;;****************************************************************************
......@@ -18,6 +22,23 @@
n-by-n real symmetric matrices. The size of the workspace
is O(2n).")
;; V 1.9
(defmobject eigen-nonsymm
"gsl_eigen_nonsymm" ((n sizet))
"non-symmetric eigenvalue workspace" ; FDL
"Make a workspace for computing eigenvalues of
n-by-n real non-symmetric matrices. The size of the workspace
is O(2n).")
(cffi:defcstruct gsl-nonsymm-ws
(size sizet) ; size of matrices
(diag :pointer) ; diagonal matrix elements from balancing
(tau :pointer) ; Householder coefficients
(Z :pointer) ; pointer to Z matrix
(balancep :int) ; perform balancing transformation?
(n-evals sizet) ; number of eigenvalues found
(francis-ws :pointer))
(defmobject eigen-symmv
"gsl_eigen_symmv" ((n sizet))
"symmetric eigensystem workspace" ; FDL
......@@ -25,6 +46,14 @@
eigenvectors of n-by-n real symmetric matrices. The size of
the workspace is O(4n).")
;; V 1.9
(defmobject eigen-nonsymmv
"gsl_eigen_nonsymmv" ((n sizet))
"non-symmetric eigenvalue workspace" ; FDL
"Make a workspace for computing for computing eigenvalues and
eigenvectors of n-by-n real nonsymmetric matrices. The size of the
workspace is O(5n).")
(defmobject eigen-herm
"gsl_eigen_herm" ((n sizet))
"Hermitian eigenvalue workspace" ; FDL
......@@ -68,6 +97,81 @@
referenced. The eigenvalues are stored in the vector eigenvalues and
are unordered.")
(defmfun set-parameters-nonsymmetric
(ws &optional compute-shur-form balance)
"gsl_eigen_nonsymm_params"
(((if compute-shur-form 1 0) :int)
((if balance 1 0) :int) ((mpointer ws) :pointer))
:gsl-version (1 9)
:c-return :void
:export nil
:index eigenvalues-nonsymm)
(defmfun eigenvalues-nonsymm
(A
&optional
(eigenvalues
(make-marray '(complex double-float) :dimensions (dim0 A)))
(ws (make-eigen-nonsymm (dim0 A)))
compute-shur-form balance shur-vectors
&aux
(sv
(if (eql shur-vectors t)
(make-marray 'double-float :dimensions (dimensions A))
shur-vectors)))
("gsl_eigen_nonsymm" "gsl_eigen_nonsymm_Z")
((((mpointer A) :pointer)
((mpointer eigenvalues) :pointer) ((mpointer ws) :pointer))
(((mpointer A) :pointer)
((mpointer eigenvalues) :pointer) ((mpointer ws) :pointer)
((mpointer sv) :pointer)))
:before
((set-parameters-nonsymmetric ws compute-shur-form balance))
:gsl-version (1 9)
:switch (shur-vectors)
:inputs (A)
:outputs (A eigenvalues)
:return
(eigenvalues
(cffi:foreign-slot-value (mpointer ws) 'gsl-nonsymm-ws 'n-evals))
:documentation ; FDL
"Compute the eigenvalues of the real nonsymmetric matrix A and
stores them in the vector 'eigenvalues. If T is desired, it is
stored in the upper portion of A on output. Otherwise, on output,
the diagonal of A will contain the 1-by-1 real eigenvalues and
2-by-2 complex conjugate eigenvalue systems, and the rest of A is
destroyed. In rare cases, this function may fail to find all
eigenvalues. If this happens, a warning is signalled and the number
of converged eigenvalues is returned as a second value. The
converged eigenvalues are stored in the beginning of eval.
If compute-shur-form is true, the full Schur form T will be computed.
If it is set to nil, T will not be computed (this is
the default setting). Computing the full Schur form requires
approximately 1.5-2 times the number of flops.
If balance is true, a balancing transformation is applied to the
matrix prior to computing eigenvalues. This transformation is
designed to make the rows and columns of the matrix have comparable
norms, and can result in more accurate eigenvalues for matrices
whose entries vary widely in magnitude. See Balancing for more
information. Note that the balancing transformation does not
preserve the orthogonality of the Schur vectors, so if you wish to
compute the Schur vectors with you will obtain the Schur vectors of
the balanced matrix instead of the original matrix. The relationship
will be
T = Q^t D^(-1) A D Q
where Q is the matrix of Schur vectors for the balanced matrix, and D
is the balancing transformation. Then this function will compute
a matrix Z which satisfies
T = Z^(-1) A Z
with Z = D Q. Note that Z will not be orthogonal. For this reason,
balancing is not performed by default.")
(defmfun eigenvalues-eigenvectors
((A matrix)
&optional
......@@ -98,6 +202,75 @@
The eigenvectors are guaranteed to be mutually orthogonal and
normalised to unit magnitude.")
(defmfun eigenvalues-eigenvectors-nonsymm
(A
&optional
(eigenvalues
(make-marray '(complex double-float) :dimensions (dim0 A)))
(eigenvectors
(make-marray '(complex double-float) :dimensions (dimensions A)))
(ws (make-eigen-nonsymmv (dim0 A)))
compute-shur-form balance shur-vectors
&aux
(sv
(if (eql shur-vectors t)
(make-marray 'double-float :dimensions (dimensions A))
shur-vectors)))
("gsl_eigen_nonsymmv" "gsl_eigen_nonsymmv_Z")
((((mpointer A) :pointer)
((mpointer eigenvalues) :pointer) ((mpointer eigenvectors) :pointer)
((mpointer ws) :pointer))
(((mpointer A) :pointer)
((mpointer eigenvalues) :pointer) ((mpointer eigenvectors) :pointer)
((mpointer ws) :pointer) ((mpointer sv) :pointer)))
:before ; this applies for evec-eval too, right?
((set-parameters-nonsymmetric ws compute-shur-form balance))
:gsl-version (1 9)
:switch (shur-vectors)
:inputs (A)
:outputs (A eigenvalues eigenvectors)
:return
(eigenvalues
(cffi:foreign-slot-value (mpointer ws) 'gsl-nonsymm-ws 'n-evals))
:documentation ; FDL
"Compute eigenvalues and right eigenvectors of the n-by-n real
nonsymmetric matrix A. It first calls #'eigenvalues-nonsymm to
compute the eigenvalues, Schur form T, and Schur vectors. Then it
finds eigenvectors of T and backtransforms them using the Schur
vectors. The Schur vectors are destroyed in the process, but can be
saved by specifying binding shur-vectors to a vector of length n, or
t to have it automatically made. The computed eigenvectors are
normalized to have unit magnitude. On output, the upper portion of A
contains the Schur form T. If #'eigenvalues-nonsymm fails, no
eigenvectors are computed, and an error code is returned.
If compute-shur-form is true, the full Schur form T will be computed.
If it is set to nil, T will not be computed (this is
the default setting). Computing the full Schur form requires
approximately 1.5-2 times the number of flops.
If balance is true, a balancing transformation is applied to the
matrix prior to computing eigenvalues. This transformation is
designed to make the rows and columns of the matrix have comparable
norms, and can result in more accurate eigenvalues for matrices
whose entries vary widely in magnitude. See Balancing for more
information. Note that the balancing transformation does not
preserve the orthogonality of the Schur vectors, so if you wish to
compute the Schur vectors with you will obtain the Schur vectors of
the balanced matrix instead of the original matrix. The relationship
will be
T = Q^t D^(-1) A D Q
where Q is the matrix of Schur vectors for the balanced matrix, and D
is the balancing transformation. Then this function will compute
a matrix Z which satisfies
T = Z^(-1) A Z
with Z = D Q. Note that Z will not be orthogonal. For this reason,
balancing is not performed by default.")
;;;;****************************************************************************
;;;; Sorting Eigenvalues and Eigenvectors
;;;;****************************************************************************
......
;; Helpers that define a single GSL function interface
;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
;; Time-stamp: <2009-01-10 20:15:42EST defmfun-single.lisp>
;; Time-stamp: <2009-01-11 10:07:06EST defmfun-single.lisp>
;; $Id: $
(in-package :gsl)
......@@ -91,7 +91,9 @@
(let ((auxstart (position '&aux arglist)))
;; &aux bindings are checked
(when auxstart
(mapcan 'rest (subseq arglist (1+ auxstart)))))))))))
(apply
'append
(mapcar 'rest (subseq arglist (1+ auxstart))))))))))))
,@(when documentation (list documentation))
#-native
,(funcall body-maker name arglist gsl-name c-arguments key-args)
......
;; Macro for defining GSL functions.
;; Liam Healy 2008-04-16 20:49:50EDT defmfun.lisp
;; Time-stamp: <2009-01-10 20:29:54EST defmfun.lisp>
;; Time-stamp: <2009-01-10 21:47:36EST defmfun.lisp>
;; $Id$
(in-package :gsl)
......@@ -182,18 +182,24 @@
selects one of two GSL functions."
(let ((optpos (position-if (lambda (s) (member s *defmfun-optk*)) arglist)))
(if optpos
(let ((mandatory-arglist (subseq arglist 0 optpos))
(optional-arglist (subseq arglist (1+ optpos))))
`(if ,(first optional-arglist)
,(body-no-optional-arg
name
(append mandatory-arglist optional-arglist)
(second gsl-name)
(second c-arguments)
key-args)
,(body-no-optional-arg
name
mandatory-arglist
(first gsl-name)
(first c-arguments)
key-args))))))
(with-defmfun-key-args key-args
(let ((mandatory-arglist (subseq arglist 0 optpos))
(optional-arglist (subseq arglist (1+ optpos))))
`(if ,(or (first switch) (first optional-arglist))
,(body-no-optional-arg
name
(append mandatory-arglist optional-arglist)
(second gsl-name)
(second c-arguments)
key-args)
,(body-no-optional-arg
name
(append mandatory-arglist
(when switch
(remove-if
(lambda (arg)
(member (if (listp arg) (first arg) arg) switch))
optional-arglist)))
(first gsl-name)
(first c-arguments)
key-args)))))))
;; Coulumb functions
;; Liam Healy, Sat Mar 18 2006 - 23:23
;; Time-stamp: <2008-12-29 20:49:45EST coulomb.lisp>
;; Time-stamp: <2009-01-11 10:14:54EST coulomb.lisp>
;; $Id$
(in-package :gsl)
......@@ -48,7 +48,9 @@
overflow occurs, the condition 'overflow is signalled and scaling
exponents are stored in the modifiable parameters exp-F, exp-G.")
(defmfun coulomb-wave-F-array (L-min eta x fc-array)
(defmfun coulomb-wave-F-array
(L-min eta x &optional (size-or-array *default-sf-array-size*)
&aux (fc-array (vdf size-or-array)))
"gsl_sf_coulomb_wave_F_array"
((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
((c-pointer fc-array) :pointer) (F-exponent :double))
......@@ -59,7 +61,12 @@
L = Lmin ... Lmin + kmax, storing the results in fc-array.
In the case of overflow the exponent is stored in the second value returned.")
(defmfun coulomb-wave-FG-array (L-min eta x fc-array gc-array)
(defmfun coulomb-wave-FG-array
(L-min eta x
&optional (fc-size-or-array *default-sf-array-size*)
gc-size-or-array
&aux (fc-array (vdf fc-size-or-array))
(gc-array (vdf (or gc-size-or-array (dim0 fc-array)))))
"gsl_sf_coulomb_wave_FG_array"
((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
((c-pointer fc-array) :pointer) ((c-pointer gc-array) :pointer)
......@@ -72,7 +79,16 @@
results in fc_array and gc_array. In the case of overflow the
exponents are stored in F_exponent and G_exponent.")
(defmfun coulomb-wave-FGp-array (L-min eta x fc-array fcp-array gc-array gcp-array)
(defmfun coulomb-wave-FGp-array
(L-min eta x
&optional (fc-size-or-array *default-sf-array-size*)
fcp-size-or-array
gc-size-or-array
gcp-size-or-array
&aux (fc-array (vdf fc-size-or-array))
(fcp-array (vdf (or fcp-size-or-array (dim0 fc-array))))
(gc-array (vdf (or gc-size-or-array (dim0 fc-array))))
(gcp-array (vdf (or gcp-size-or-array (dim0 fc-array)))))
"gsl_sf_coulomb_wave_FGp_array"
((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
((c-pointer fc-array) :pointer) ((c-pointer fcp-array) :pointer)
......@@ -90,7 +106,9 @@
In the case of overflow the exponents are stored in F_exponent
and G_exponent.")
(defmfun coulomb-wave-sphF-array (L-min eta x fc-array)
(defmfun coulomb-wave-sphF-array
(L-min eta x &optional (size-or-array *default-sf-array-size*)
&aux (fc-array (vdf size-or-array)))
"gsl_sf_coulomb_wave_sphF_array"
((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
((c-pointer fc-array) :pointer) (F-exponent :double))
......@@ -113,14 +131,16 @@
"The Coulomb wave function normalization constant C_L(\eta)
for L > -1.")
(defmfun coulomb-CL-array (L-min eta cl)
(defmfun coulomb-CL-array
(L-min eta &optional (size-or-array *default-sf-array-size*)
&aux (array (vdf size-or-array)))
"gsl_sf_coulomb_CL_array"
((L-min :double) ((1- (dim0 cl)) :int) (eta :double)
((c-pointer cl) :pointer))
((L-min :double) ((1- (dim0 array)) :int) (eta :double)
((c-pointer array) :pointer))
:outputs (array)
:documentation ; FDL
"The Coulomb wave function normalization constant C_L(\eta)
for L = Lmin ... Lmin + kmax, Lmin > -1."
:outputs (cl))
for L = Lmin ... Lmin + kmax, Lmin > -1.")
;;;;****************************************************************************
;;;; Examples and unit test
......
;; Gegenbauer polynomials
;; Liam Healy, Fri Apr 28 2006 - 20:40
;; Time-stamp: <2008-12-26 12:25:31EST gegenbauer.lisp>
;; Time-stamp: <2009-01-11 09:45:00EST gegenbauer.lisp>
;; $Id$
(in-package :gsl)
......@@ -27,17 +27,16 @@
"The Gegenbauer polynomial C^{(\lambda)}_n(x)} for a specific value of n,
lambda, x subject to \lambda > -1/2, n >= 0.")
(defmfun gegenbauer-array (lambda x result)
(defmfun gegenbauer-array
(lambda x &optional (size-or-array *default-sf-array-size*)
&aux (array (vdf size-or-array)))
"gsl_sf_gegenpoly_array"
(((1- (dim0 result)) :int)
(lambda :double) (x :double) ((c-pointer result) :pointer))
:outputs (result)
(((1- (dim0 array)) :int)
(lambda :double) (x :double) ((c-pointer array) :pointer))
:outputs (array)
:documentation ; FDL
"Compute an array of Gegenbauer polynomials C^{(\lambda)}_n(X)}
for n = 0, 1, 2, ..., length(result)-1}, subject to \lambda > -1/2.")
;;; (defparameter vec (make-data 'vector nil 3))
;;; (gegenbauer-array 1.0d0 3.0d0 vec)
for n = 0, 1, 2, ..., length(array)-1}, subject to \lambda > -1/2.")
;;;;****************************************************************************
;;;; Examples and unit test
......@@ -48,5 +47,4 @@
(gegenbauer-2 1.0d0 3.0d0)
(gegenbauer-3 1.0d0 3.0d0)
(gegenbauer 4 1.0d0 3.0d0)
(let ((arr (make-marray 'double-float :dimensions 4)))
(gegenbauer-array 1.0d0 3.0d0 arr) (cl-array arr)))
(cl-array (gegenbauer-array 1.0d0 3.0d0 4)))
;; Legendre functions
;; Liam Healy, Sat Apr 29 2006 - 19:16
;; Time-stamp: <2008-12-26 11:52:33EST legendre.lisp>
;; Time-stamp: <2009-01-11 10:21:24EST legendre.lisp>
;; $Id$
(in-package :gsl)
......@@ -35,21 +35,25 @@
"The Legendre polynomial P_l(x) for a specific value of l,
x subject to l >= 0, |x| <= 1.")
(defmfun legendre-Pl-array (x array)
(defmfun legendre-Pl-array
(x &optional (size-or-array *default-sf-array-size*)
&aux (array (vdf size-or-array)))
"gsl_sf_legendre_Pl_array"
(((1- (dim0 array)) :int) (x :double) ((c-pointer array) :pointer))
:outputs (array)
:documentation ; FDL
"Compute an array of Legendre polynomials
P_l(x) for l = 0, ..., length(array), |x| <= 1."
:outputs (array))
P_l(x) for l = 0, ..., length(array), |x| <= 1.")
(defmfun legendre-Pl-deriv-array (x array)
(defmfun legendre-Pl-deriv-array
(x &optional (size-or-array *default-sf-array-size*)
&aux (array (vdf size-or-array)))
"gsl_sf_legendre_Pl_deriv_array"
(((1- (dim0 array)) :int) (x :double) ((c-pointer array) :pointer))
:outputs (array)
:documentation ; FDL
"Compute an array of Legendre polynomials derivatives
dP_l(x)/dx, for l = 0, ..., length(array), |x| <= 1."
:outputs (array))
dP_l(x)/dx, for l = 0, ..., length(array), |x| <= 1.")
(defmfun legendre-Q0 (x)
"gsl_sf_legendre_Q0_e" ((x :double) (ret sf-result))
......@@ -91,25 +95,31 @@
"The associated Legendre polynomial
P_l^m(x) for m >= 0, l >= m, |x| <= 1.")
(defmfun legendre-Plm-array (m x array)
(defmfun legendre-Plm-array
(m x &optional (size-or-array *default-sf-array-size*)
&aux (array (vdf size-or-array)))
"gsl_sf_legendre_Plm_array"
(((+ (dim0 array) m -1) :int) (m :int) (x :double)
((c-pointer array) :pointer))
:outputs (array)
:documentation ; FDL
"An array of Legendre polynomials
P_l^m(x), for m >= 0,
l = |m|, ..., |m|+length(array)-1} and |x| <= 1."
:outputs (array))
l = |m|, ..., |m|+length(array)-1} and |x| <= 1.")
(defmfun legendre-Plm-deriv-array (m x values derivatives)
(defmfun legendre-Plm-deriv-array
(m x &optional (values-size-or-array *default-sf-array-size*)
(derivatives-size-or-array *default-sf-array-size*)
&aux (values (vdf values-size-or-array))
(derivatives (vdf derivatives-size-or-array)))
"gsl_sf_legendre_Plm_deriv_array"
(((+ (dim0 values) m -1) :int) (m :int) (x :double)
((c-pointer values) :pointer) ((c-pointer derivatives) :pointer))
:outputs (values derivatives)
:documentation ; FDL
"An array of Legendre polynomials
values and derivatives dP_l^m(x)/dx for m >= 0,
l = |m|, ..., length(values) and |x| <= 1."
:outputs (values derivatives))
l = |m|, ..., length(values) and |x| <= 1.")
(defmfun legendre-sphPlm (l m x)
"gsl_sf_legendre_sphPlm_e" ((l :int) (m :int) (x :double) (ret sf-result))
......@@ -120,25 +130,31 @@
m >= 0, l >= m, |x| <= 1. These routines avoid the overflows
that occur for the standard normalization of P_l^m(x).")
(defmfun legendre-sphPlm-array (m x array)
(defmfun legendre-sphPlm-array
(m x &optional (size-or-array *default-sf-array-size*)
&aux (array (vdf size-or-array)))
"gsl_sf_legendre_sphPlm_array"
(((+ (dim0 array) m -1) :int) (m :int) (x :double)
((c-pointer array) :pointer))
:outputs (array)
:documentation ; FDL
"An array of normalized associated Legendre functions
\sqrt(2l+1)/(4\pi) \sqrt(l-m)!/(l+m)! P_l^m(x),
for m >= 0, l = |m|, ..., length(array)}, |x| <= 1.0."
:outputs (array))
for m >= 0, l = |m|, ..., length(array)}, |x| <= 1.0.")
(defmfun legendre-sphPlm-deriv-array (m x values derivatives)
(defmfun legendre-sphPlm-deriv-array
(m x &optional (values-size-or-array *default-sf-array-size*)
(derivatives-size-or-array *default-sf-array-size*)
&aux (values (vdf values-size-or-array))
(derivatives (vdf derivatives-size-or-array)))
"gsl_sf_legendre_sphPlm_deriv_array"
(((+ (dim0 values) m -1) :int) (m :int) (x :double)
((c-pointer values) :pointer) ((c-pointer derivatives) :pointer))
:outputs (values derivatives)
:documentation ; FDL
"An array of normalized associated Legendre functions
values and derivatives for m >= 0,
l = |m|, ..., length(array)}, |x| <= 1.0."
:outputs (values derivatives))
l = |m|, ..., length(array)}, |x| <= 1.0.")
(defmfun legendre-array-size (lmax m)
"gsl_sf_legendre_array_size" ((lmax :int) (m :int))
......@@ -234,7 +250,9 @@
\eta >= 0, l >= 0. In the flat limit this takes the form
L^{H3d}_l(\lambda,\eta) = j_l(\lambda\eta).")
(defmfun legendre-H3d-array (lambda eta array)
(defmfun legendre-H3d-array
(lambda eta &optional (size-or-array *default-sf-array-size*)
&aux (array (vdf size-or-array)))
"gsl_sf_legendre_H3d_array"
(((1- (dim0 array)) :int) (lambda :double) (eta :double)
((c-pointer array) :pointer))
......@@ -243,10 +261,6 @@
"An array of radial eigenfunctions
L^{H3d}_l(\lambda, \eta) for 0 <= l <= length(array).")
;;; (defparameter hleg (make-data 'vector nil 3))
;;; (legendre-H3d-array 1.0d0 0.5d0 hleg)
;;; #<GSL-VECTOR #(0.9200342692589383d0 0.21694026450392123d0 0.047950660488307775d0) {C07CB51}>
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
......
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