Skip to content
Snippets Groups Projects
Commit 6f51bb7f authored by Liam Healy's avatar Liam Healy
Browse files

Eigensystems port to ffa

Eigenvalues, eigenvectors ported to ffa, using generic functions to
choose between real symmetrix and complex hermitian matrices.  The
example has been tried and works.
parent 383cfb06
No related branches found
No related tags found
No related merge requests found
;; Eigenvectors and eigenvalues ;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52 ;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2008-03-10 21:20:59EDT eigensystems.lisp> ;; Time-stamp: <2008-08-14 22:48:02EDT eigensystems.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Real Symmetric Matrices ;;;; Workspace
;;;;**************************************************************************** ;;;;****************************************************************************
(defgo-s (eigen-symm n) eigen-symm-alloc eigen-symm-free) (defgo-s (eigen-symm n) eigen-symm-alloc eigen-symm-free)
(defmfun eigen-symm-alloc (n) (defmfun eigen-symm-alloc (n)
"gsl_eigen_symm_alloc" ((n size)) "gsl_eigen_symm_alloc" ((n sizet))
:c-return :pointer :c-return :pointer
:export nil :export nil
:index '(letm eigen-symm) :index '(letm eigen-symm)
...@@ -32,7 +32,7 @@ ...@@ -32,7 +32,7 @@
(defgo-s (eigen-symmv n) eigen-symmv-alloc eigen-symmv-free) (defgo-s (eigen-symmv n) eigen-symmv-alloc eigen-symmv-free)
(defmfun eigen-symmv-alloc (n) (defmfun eigen-symmv-alloc (n)
"gsl_eigen_symmv_alloc" ((n size)) "gsl_eigen_symmv_alloc" ((n sizet))
:c-return :pointer :c-return :pointer
:export nil :export nil
:index '(letm eigen-symmv) :index '(letm eigen-symmv)
...@@ -48,44 +48,10 @@ ...@@ -48,44 +48,10 @@
:documentation ; FDL :documentation ; FDL
"Free the memory associated with the workspace w.") "Free the memory associated with the workspace w.")
(defmfun eigenvalues-symmetric (A eval ws)
"gsl_eigen_symm" ((A gsl-matrix-c) (eval gsl-vector-c) (ws :pointer))
:documentation ; FDL
"Eigenvalues of the real symmetric matrix
A. Additional workspace of the appropriate size must be provided
in w. The diagonal and lower triangular part of A are
destroyed during the computation, but the strict upper triangular part
is not referenced. The eigenvalues are stored in the vector eval
and are unordered."
:invalidate (A eval)
:return (eval))
(defmfun eigenvalues-eigenvectors-symmetric (A eval evec ws)
"gsl_eigen_symmv"
(((pointer A) :pointer) ((pointer eval) :pointer)
((pointer evec) :pointer) (ws :pointer))
:documentation ; FDL
"The eigenvalues and eigenvectors of the real
symmetric matrix A. Additional workspace of the appropriate size
must be provided in w. The diagonal and lower triangular part of
A are destroyed during the computation, but the strict upper
triangular part is not referenced. The eigenvalues are stored in the
vector eval and are unordered. The corresponding eigenvectors are
stored in the columns of the matrix 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 (eval evec))
;;;;****************************************************************************
;;;; Complex Hermitian Matrices
;;;;****************************************************************************
(defgo-s (eigen-herm n) eigen-herm-alloc eigen-herm-free) (defgo-s (eigen-herm n) eigen-herm-alloc eigen-herm-free)
(defmfun eigen-herm-alloc (n) (defmfun eigen-herm-alloc (n)
"gsl_eigen_herm_alloc" ((n size)) "gsl_eigen_herm_alloc" ((n sizet))
:documentation ; FDL :documentation ; FDL
"Allocate a workspace for computing eigenvalues of "Allocate a workspace for computing eigenvalues of
n-by-n complex hermitian matrices. The size of the workspace n-by-n complex hermitian matrices. The size of the workspace
...@@ -98,24 +64,10 @@ ...@@ -98,24 +64,10 @@
:documentation ; FDL :documentation ; FDL
"Free the memory associated with the workspace w.") "Free the memory associated with the workspace w.")
(defmfun eigenvalues-hermitian (A eval w)
"gsl_eigen_herm"
(((pointer A) gsl-matrix-c) ((pointer eval) gsl-vector-c) (w :pointer))
:documentation ; FDL
"Compute the eigenvalues of the complex hermitian matrix
A. Additional workspace of the appropriate size must be provided
in w. The diagonal and lower triangular part of A are
destroyed during the computation, but the strict upper triangular part
is not referenced. The imaginary parts of the diagonal are assumed to be
zero and are not referenced. The eigenvalues are stored in the vector
eval and are unordered."
:invalidate (eval A)
:return (eval))
(defgo-s (eigen-hermv n) eigen-hermv-alloc eigen-hermv-free) (defgo-s (eigen-hermv n) eigen-hermv-alloc eigen-hermv-free)
(defmfun eigen-hermv-alloc (n) (defmfun eigen-hermv-alloc (n)
"gsl_eigen_hermv_alloc" ((n size)) "gsl_eigen_hermv_alloc" ((n sizet))
:documentation ; FDL :documentation ; FDL
"Allocate a workspace for computing eigenvalues and "Allocate a workspace for computing eigenvalues and
eigenvectors of n-by-n complex hermitian matrices. The size of eigenvectors of n-by-n complex hermitian matrices. The size of
...@@ -128,69 +80,93 @@ ...@@ -128,69 +80,93 @@
:documentation ; FDL :documentation ; FDL
"Free the memory associated with the workspace w.") "Free the memory associated with the workspace w.")
(defmfun eigenvalues-eigenvectors-hermitian (A eval evec w) ;;;;****************************************************************************
"gsl_eigen_hermv" ;;;; Eigenvalues and eigenvectors
(((pointer A) gsl-matrix-c) ((pointer eval) gsl-vector-c) ;;;;****************************************************************************
((pointer evec) gsl-matrix-c) (w :pointer))
(defmfun eigenvalues ((A matrix) eigenvalues ws)
(double-float "gsl_eigen_symm"
complex-double-float "gsl_eigen_herm")
(((mpointer A) :pointer)
((mpointer eigenvalues) :pointer) (ws :pointer))
:definition :generic
:element-types :doubles
:inputs (A)
:outputs (A eigenvalues)
:return (eigenvalues)
:documentation ; FDL :documentation ; FDL
"Compute the eigenvalues and eigenvectors of the complex "Eigenvalues of the real symmetric or complex hermitian matrix A.
Additional workspace of the appropriate size and type must be
provided in w. The diagonal and lower triangular part of A are
destroyed during the computation, but the strict upper triangular
part is not referenced. For the complex hermitian case, The
imaginary parts of the diagonal are assumed to be zero and are not
referenced. The eigenvalues are stored in the vector eigenvalues and
are unordered.")
(defmfun eigenvalues-eigenvectors ((A matrix) eigenvalues eigenvectors ws)
(double-float "gsl_eigen_symmv"
complex-double-float "gsl_eigen_hermv")
(((mpointer A) :pointer) ((mpointer eigenvalues) :pointer)
((mpointer eigenvectors) :pointer) (ws :pointer))
:definition :generic
:element-types :doubles
:inputs (A)
:outputs (A eigenvalues eigenvectors)
:return (eigenvalues eigenvectors)
:documentation ; FDL
"The eigenvalues and eigenvectors of the real symmetric or complex
hermitian matrix A. Additional workspace of the appropriate size hermitian matrix A. Additional workspace of the appropriate size
must be provided in w. The diagonal and lower triangular part of must be provided in w. The diagonal and lower triangular part of A
A are destroyed during the computation, but the strict upper are destroyed during the computation, but the strict upper
triangular part is not referenced. The imaginary parts of the diagonal triangular part is not referenced. For complex hermitian matrices,
are assumed to be zero and are not referenced. The eigenvalues are the imaginary parts of the diagonal are assumed to be zero and are
stored in the vector eval and are unordered. The corresponding not referenced. The eigenvalues are stored in the vector
complex eigenvectors are stored in the columns of the matrix evec. eigenvalues and are unordered. The corresponding eigenvectors are
For example, the eigenvector in the first column corresponds to the stored in the columns of the matrix eigenvectors. For example, the
first eigenvalue. The eigenvectors are guaranteed to be mutually eigenvector in the first column corresponds to the first eigenvalue.
orthogonal and normalised to unit magnitude." The eigenvectors are guaranteed to be mutually orthogonal and
:invalidate (A eval evec) normalised to unit magnitude.")
:return (eval evec))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Sorting Eigenvalues and Eigenvectors ;;;; Sorting Eigenvalues and Eigenvectors
;;;;**************************************************************************** ;;;;****************************************************************************
(cffi:defcenum eigen-sort-type (cffi:defcenum eigen-sort-type
"gsl_eigen_sort_t from /usr/include/gsl/gsl_eigen.h." ;; gsl_eigen_sort_t from /usr/include/gsl/gsl_eigen.h.
:value-ascending :value-descending :absolute-ascending :absolute-descending) :value-ascending :value-descending :absolute-ascending :absolute-descending)
(export 'sort-eigenvalues-eigenvectors) (defmfun sort-eigenvalues-eigenvectors
(defgeneric sort-eigenvalues-eigenvectors (eval evec sort-type) ((eigenvalues vector) (eigenvectors matrix) sort-type)
(:documentation ; FDL (double-float "gsl_eigen_symmv_sort"
complex-double-float "gsl_eigen_hermv_sort")
(((mpointer eigenvalues) :pointer) ((mpointer eigenvectors) :pointer)
(sort-type eigen-sort-type))
:definition :generic
:element-types :doubles
:inputs (eigenvalues eigenvectors)
:outputs (eigenvalues eigenvectors)
:documentation ; FDL
"Simultaneously sort the eigenvalues stored in the vector "Simultaneously sort the eigenvalues stored in the vector
eval and the corresponding real eigenvectors stored in the columns eigenvalues and the corresponding real eigenvectors stored in the columns
of the matrix evec into ascending or descending order according to of the matrix eigenvectors into ascending or descending order according to
the value of the parameter sort-type: :value-ascending, the value of the parameter sort-type: :value-ascending,
:value-descending, :absolute-ascending, :absolute-descending.")) :value-descending, :absolute-ascending, :absolute-descending.")
(defmfun sort-eigenvalues-eigenvectors
((eval vector-double-float) evec sort-type)
"gsl_eigen_symmv_sort"
((eval gsl-vector-c) (evec gsl-matrix-c) (sort-type eigen-sort-type))
:type :method
:invalidate (eval evec))
(defmfun sort-eigenvalues-eigenvectors ((eval vector-complex) evec sort-type)
"gsl_eigen_hermv_sort"
((eval gsl-vector-c) (evec gsl-matrix-c) (sort-type eigen-sort-type))
:type :method
:invalidate (eval evec))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Example ;;;; Example
;;;;**************************************************************************** ;;;;****************************************************************************
(defun eigenvalue-eigenvectors-example () (defun eigenvalue-eigenvectors-example ()
(letm ((evecs (matrix-double-float 3 3)) (letm ((evecs (matrix-double-float '(3 3)))
(evals (vector-double-float 3)) (evals (vector-double-float 3))
(ws (eigen-symmv 3)) (ws (eigen-symmv 3))
(mat (matrix-double-float (mat (matrix-double-float
#2A((20.0d0 -10.0d0 0.0d0) (a (20.0d0 -10.0d0 0.0d0)
(-10.0d0 30.0d0 0.0d0) (-10.0d0 30.0d0 0.0d0)
(0.0d0 0.0d0 40.0d0))))) (0.0d0 0.0d0 40.0d0)))))
(eigenvalues-eigenvectors-symmetric mat evals evecs ws) (eigenvalues-eigenvectors mat evals evecs ws)
(values (data evals) (data evecs)))) (values (cl-array evals) (cl-array evecs))))
#| #|
(make-tests eigensystems (make-tests eigensystems
......
;; Definition of GSLL system ;; Definition of GSLL system
;; Liam Healy ;; Liam Healy
;; Time-stamp: <2008-08-12 22:02:21EDT gsll.asd> ;; Time-stamp: <2008-08-14 22:48:37EDT gsll.asd>
;; $Id$ ;; $Id$
(asdf:defsystem "gsll" (asdf:defsystem "gsll"
...@@ -90,7 +90,6 @@ ...@@ -90,7 +90,6 @@
(:file "cholesky") (:file "cholesky")
(:file "diagonal") (:file "diagonal")
(:file "householder"))) (:file "householder")))
#+no
(:file "eigensystems" :depends-on (init data)) (:file "eigensystems" :depends-on (init data))
;; Skip fft for now, I'm not sure how it works in C ;; Skip fft for now, I'm not sure how it works in C
#+no #+no
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment