eigensystems.lisp 7.49 KB
Newer Older
1 2
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
3
;; Time-stamp: <2008-02-03 14:14:10EST eigensystems.lisp>
4
;; $Id: $
5 6 7 8 9 10 11

(in-package :gsl)

;;;;****************************************************************************
;;;; Real Symmetric Matrices
;;;;****************************************************************************

12
(defgo-s (eigen-symm n) eigen-symm-alloc eigen-symm-free)
13

liam's avatar
liam committed
14 15
(defun-gsl eigen-symm-alloc (n)
  "gsl_eigen_symm_alloc" ((n :size))
16 17
  :c-return :pointer
  :export nil
18
  :index '(letm eigen-symm)
19 20 21
  :documentation			; FDL
  "Allocate a workspace for computing eigenvalues of
  n-by-n real symmetric matrices.  The size of the workspace
22
  is O(2n).")
23

liam's avatar
liam committed
24 25 26
(defun-gsl eigen-symm-free (w)
  "gsl_eigen_symm_free" ((w :pointer))
  :c-return :void
27
  :export nil
28
  :index '(letm eigen-symm)
29 30
  :documentation			; FDL
  "Free the memory associated with the workspace w.")
31

32
(defgo-s (eigen-symmv n) eigen-symmv-alloc eigen-symmv-free)
33

liam's avatar
liam committed
34 35
(defun-gsl eigen-symmv-alloc (n)
  "gsl_eigen_symmv_alloc" ((n :size))
36 37
  :c-return :pointer
  :export nil
38
  :index '(letm eigen-symmv)
39 40 41
  :documentation			; FDL
  "Allocate a workspace for computing eigenvalues and
  eigenvectors of n-by-n real symmetric matrices.  The size of
42
  the workspace is O(4n).")
43

liam's avatar
liam committed
44 45
(defun-gsl eigen-symmv-free (w)
  "gsl_eigen_symmv_free" ((w :pointer))
46
  :index '(letm eigen-symmv)
liam's avatar
liam committed
47
  :c-return :void
48 49
  :documentation			; FDL
  "Free the memory associated with the workspace w.")
50

liam's avatar
liam committed
51 52
(defun-gsl eigenvalues-symmetric (A eval ws)
  "gsl_eigen_symm" ((A gsl-matrix-c) (eval gsl-vector-c) (ws :pointer))
53 54 55 56
  :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
57 58 59 60
  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)
liam's avatar
liam committed
61
  :return (eval))
62

liam's avatar
liam committed
63
(defun-gsl eigenvalues-eigenvectors-symmetric (A eval evec ws)
64
  "gsl_eigen_symmv"
65 66
  (((pointer A) :pointer) ((pointer eval) :pointer)
   ((pointer evec) :pointer) (ws :pointer))
67 68 69 70 71
  :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
72
  triangular part is not referenced.  The eigenvalues are stored in the
73 74
  vector eval and are unordered.  The corresponding eigenvectors are
  stored in the columns of the matrix evec.  For example, the
75 76 77 78
  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)
liam's avatar
liam committed
79
  :return (eval evec))
80 81 82 83 84

;;;;****************************************************************************
;;;; Complex Hermitian Matrices
;;;;****************************************************************************

85
(defgo-s (eigen-herm n) eigen-herm-alloc eigen-herm-free)
86

liam's avatar
liam committed
87 88
(defun-gsl eigen-herm-alloc (n)
  "gsl_eigen_herm_alloc" ((n :size))
89 90 91 92
  :documentation			; FDL
  "Allocate a workspace for computing eigenvalues of
  n-by-n complex hermitian matrices.  The size of the workspace
  is O(3n)."
liam's avatar
liam committed
93
  :c-return :pointer)
94

liam's avatar
liam committed
95 96 97
(defun-gsl eigen-herm-free (w)
  "gsl_eigen_herm_free" ((w :pointer))
  :c-return :void
98 99
  :documentation			; FDL
  "Free the memory associated with the workspace w.")
100

liam's avatar
liam committed
101 102 103
(defun-gsl eigenvalues-hermitian (A eval w)
  "gsl_eigen_herm"
  (((pointer A) gsl-matrix-c) ((pointer eval) gsl-vector-c) (w :pointer))
104
  :documentation			; FDL
liam's avatar
liam committed
105
  "Compute the eigenvalues of the complex hermitian matrix
106 107
   A.  Additional workspace of the appropriate size must be provided
   in w.  The diagonal and lower triangular part of A are
liam's avatar
liam committed
108 109 110 111 112 113 114
   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
   @var{eval} and are unordered."
  :invalidate (eval A)
  :return (eval))

115
(defgo-s (eigen-hermv n) eigen-hermv-alloc eigen-hermv-free)
116

liam's avatar
liam committed
117 118
(defun-gsl eigen-hermv-alloc (n)
  "gsl_eigen_hermv_alloc" ((n :size))
119 120 121 122
  :documentation			; FDL
  "Allocate a workspace for computing eigenvalues and
  eigenvectors of n-by-n complex hermitian matrices.  The size of
  the workspace is O(5n)."
liam's avatar
liam committed
123
  :c-return :pointer)
124

liam's avatar
liam committed
125 126 127
(defun-gsl eigen-hermv-free (w)
  "gsl_eigen_hermv_free" ((w :pointer))
  :c-return :void
128 129
  :documentation			; FDL
  "Free the memory associated with the workspace w.")
130

liam's avatar
liam committed
131 132 133 134
(defun-gsl 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))
135
  :documentation			; FDL
liam's avatar
liam committed
136
  "Compute the eigenvalues and eigenvectors of the complex
137 138 139
  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
liam's avatar
liam committed
140 141
  triangular part is not referenced. The imaginary parts of the diagonal
  are assumed to be zero and are not referenced.  The eigenvalues are
142 143
  stored in the vector eval and are unordered.  The corresponding
  complex eigenvectors are stored in the columns of the matrix evec.
liam's avatar
liam committed
144 145 146 147 148
  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))
149 150 151 152 153 154 155 156 157

;;;;****************************************************************************
;;;; 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)

158 159 160 161 162 163 164 165 166
(export 'sort-eigenvalues-eigenvectors)
(defgeneric sort-eigenvalues-eigenvectors (eval evec sort-type)
  (: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
  the value of the parameter sort-type: :value-ascending,
  :value-descending, :absolute-ascending, :absolute-descending."))

167
(defun-gsl sort-eigenvalues-eigenvectors
liam's avatar
liam committed
168
    ((eval gsl-vector-double) evec sort-type)
169
  "gsl_eigen_symmv_sort"
liam's avatar
liam committed
170 171 172 173 174 175 176 177 178 179
  ((eval gsl-vector-c) (evec gsl-matrix-c) (sort-type eigen-sort-type))
  :type :method
  :invalidate (eval evec))

(defun-gsl sort-eigenvalues-eigenvectors
    ((eval gsl-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))
180 181 182 183 184

;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************
 
185
(defun eigenvalue-eigenvectors-example ()
186 187 188 189 190 191 192 193 194
  (letm ((evecs (matrix-double 3 3))
	 (evals (vector-double 3))
	 (ws (eigen-symmv 3))
	 (mat (matrix-double
	       #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))))
195 196

#|
197 198 199 200 201 202
#(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))

|#