From 6f51bb7f180a04fe98c33c7a8b7f4fe8c6c962d3 Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Thu, 14 Aug 2008 22:51:12 -0400 Subject: [PATCH] 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. --- eigensystems.lisp | 170 ++++++++++++++++++++-------------------------- gsll.asd | 3 +- 2 files changed, 74 insertions(+), 99 deletions(-) diff --git a/eigensystems.lisp b/eigensystems.lisp index 8c1665dc..415da6af 100644 --- a/eigensystems.lisp +++ b/eigensystems.lisp @@ -1,18 +1,18 @@ ;; Eigenvectors and eigenvalues ;; 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$ (in-package :gsl) ;;;;**************************************************************************** -;;;; Real Symmetric Matrices +;;;; Workspace ;;;;**************************************************************************** (defgo-s (eigen-symm n) eigen-symm-alloc eigen-symm-free) (defmfun eigen-symm-alloc (n) - "gsl_eigen_symm_alloc" ((n size)) + "gsl_eigen_symm_alloc" ((n sizet)) :c-return :pointer :export nil :index '(letm eigen-symm) @@ -32,7 +32,7 @@ (defgo-s (eigen-symmv n) eigen-symmv-alloc eigen-symmv-free) (defmfun eigen-symmv-alloc (n) - "gsl_eigen_symmv_alloc" ((n size)) + "gsl_eigen_symmv_alloc" ((n sizet)) :c-return :pointer :export nil :index '(letm eigen-symmv) @@ -48,44 +48,10 @@ :documentation ; FDL "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) (defmfun eigen-herm-alloc (n) - "gsl_eigen_herm_alloc" ((n size)) + "gsl_eigen_herm_alloc" ((n sizet)) :documentation ; FDL "Allocate a workspace for computing eigenvalues of n-by-n complex hermitian matrices. The size of the workspace @@ -98,24 +64,10 @@ :documentation ; FDL "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) (defmfun eigen-hermv-alloc (n) - "gsl_eigen_hermv_alloc" ((n size)) + "gsl_eigen_hermv_alloc" ((n sizet)) :documentation ; FDL "Allocate a workspace for computing eigenvalues and eigenvectors of n-by-n complex hermitian matrices. The size of @@ -128,69 +80,93 @@ :documentation ; FDL "Free the memory associated with the workspace w.") -(defmfun eigenvalues-eigenvectors-hermitian (A eval evec w) - "gsl_eigen_hermv" - (((pointer A) gsl-matrix-c) ((pointer eval) gsl-vector-c) - ((pointer evec) gsl-matrix-c) (w :pointer)) +;;;;**************************************************************************** +;;;; Eigenvalues and eigenvectors +;;;;**************************************************************************** + +(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 - "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 - 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. The corresponding - complex 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)) + 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 complex hermitian matrices, + 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. The corresponding eigenvectors are + stored in the columns of the matrix eigenvectors. 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.") ;;;;**************************************************************************** ;;;; Sorting Eigenvalues and Eigenvectors ;;;;**************************************************************************** (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) -(export 'sort-eigenvalues-eigenvectors) -(defgeneric sort-eigenvalues-eigenvectors (eval evec sort-type) - (:documentation ; FDL +(defmfun sort-eigenvalues-eigenvectors + ((eigenvalues vector) (eigenvectors matrix) sort-type) + (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 - eval and the corresponding real eigenvectors stored in the columns - of the matrix evec into ascending or descending order according to + eigenvalues and the corresponding real eigenvectors stored in the columns + of the matrix eigenvectors into ascending or descending order according to the value of the parameter sort-type: :value-ascending, - :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)) + :value-descending, :absolute-ascending, :absolute-descending.") ;;;;**************************************************************************** ;;;; Example ;;;;**************************************************************************** (defun eigenvalue-eigenvectors-example () - (letm ((evecs (matrix-double-float 3 3)) + (letm ((evecs (matrix-double-float '(3 3))) (evals (vector-double-float 3)) (ws (eigen-symmv 3)) - (mat (matrix-double-float - #2A((20.0d0 -10.0d0 0.0d0) - (-10.0d0 30.0d0 0.0d0) - (0.0d0 0.0d0 40.0d0))))) - (eigenvalues-eigenvectors-symmetric mat evals evecs ws) - (values (data evals) (data evecs)))) + (mat (matrix-double-float + (a (20.0d0 -10.0d0 0.0d0) + (-10.0d0 30.0d0 0.0d0) + (0.0d0 0.0d0 40.0d0))))) + (eigenvalues-eigenvectors mat evals evecs ws) + (values (cl-array evals) (cl-array evecs)))) #| (make-tests eigensystems diff --git a/gsll.asd b/gsll.asd index c443bab7..7029ba87 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2008-08-12 22:02:21EDT gsll.asd> +;; Time-stamp: <2008-08-14 22:48:37EDT gsll.asd> ;; $Id$ (asdf:defsystem "gsll" @@ -90,7 +90,6 @@ (:file "cholesky") (:file "diagonal") (:file "householder"))) - #+no (:file "eigensystems" :depends-on (init data)) ;; Skip fft for now, I'm not sure how it works in C #+no -- GitLab