eigensystems.lisp 11.4 KB
Newer Older
1 2
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2009-01-12 10:22:00EST eigensystems.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_eigen.h

10 11 12 13
;;; Should symmetric matrices form a subclass of matrices, so that
;;; both eigenvalues and eigenvalues-nonsymm could be methods of the
;;; same function?

14
;;;;****************************************************************************
Liam Healy's avatar
Liam Healy committed
15
;;;; Workspace 
16 17
;;;;****************************************************************************

18 19 20 21 22 23 24
(defmobject eigen-symm
    "gsl_eigen_symm" ((n sizet))
    "symmetric eigenvalue workspace"	; FDL
    "Make a workspace for computing eigenvalues of
  n-by-n real symmetric matrices.  The size of the workspace
  is O(2n).")

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
;; 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))

42 43 44 45 46 47 48
(defmobject eigen-symmv
    "gsl_eigen_symmv" ((n sizet))
    "symmetric eigensystem workspace"	; FDL
    "Make a workspace for computing eigenvalues and
  eigenvectors of n-by-n real symmetric matrices.  The size of
  the workspace is O(4n).")

49 50 51 52 53 54 55 56
;; 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).")

57 58 59 60 61 62
(defmobject eigen-herm
    "gsl_eigen_herm" ((n sizet))
    "Hermitian eigenvalue workspace"	; FDL
    "Make a workspace for computing eigenvalues of
  n-by-n complex Hermitian matrices.  The size of the workspace
  is O(3n).")
63

64 65 66 67 68 69
(defmobject eigen-hermv
    "gsl_eigen_hermv" ((n sizet))
    "Hermitian eigensystem workspace"	; FDL
    "Make a workspace for computing eigenvalues and
  eigenvectors of n-by-n complex hermitian matrices.  The size of
  the workspace is O(5n).")
70

Liam Healy's avatar
Liam Healy committed
71 72 73 74
;;;;****************************************************************************
;;;; Eigenvalues and eigenvectors
;;;;****************************************************************************

75 76 77
(defmfun eigenvalues
    ((A matrix)
     &optional
78 79 80
     (eigenvalues (make-marray element-type :dimensions (dim0 A)))
     (ws (eltcase complex (make-eigen-herm (dim0 A))
		  t (make-eigen-symm (dim0 A)))))
Liam Healy's avatar
Liam Healy committed
81
  (double-float "gsl_eigen_symm"
82
		complex-double-float "gsl_eigen_herm")
Liam Healy's avatar
Liam Healy committed
83
  (((mpointer A) :pointer)
84
   ((mpointer eigenvalues) :pointer) ((mpointer ws) :pointer))
Liam Healy's avatar
Liam Healy committed
85 86 87 88 89
  :definition :generic
  :element-types :doubles
  :inputs (A)
  :outputs (A eigenvalues)
  :return (eigenvalues)
90
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
91 92 93 94 95 96 97 98 99
  "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.")

100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
(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.")

175 176 177
(defmfun eigenvalues-eigenvectors
    ((A matrix)
     &optional
178 179 180 181
     (eigenvalues (make-marray element-type :dimensions (dim0 A)))
     (eigenvectors (make-marray element-type :dimensions (dimensions A)))
     (ws (eltcase complex (make-eigen-hermv (dim0 A))
		  t (make-eigen-symmv (dim0 A)))))
Liam Healy's avatar
Liam Healy committed
182
  (double-float "gsl_eigen_symmv"
183
		complex-double-float "gsl_eigen_hermv")  
Liam Healy's avatar
Liam Healy committed
184
  (((mpointer A) :pointer) ((mpointer eigenvalues) :pointer)
185
   ((mpointer eigenvectors) :pointer) ((mpointer ws) :pointer))
Liam Healy's avatar
Liam Healy committed
186 187 188 189 190 191 192
  :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
193
  hermitian matrix A.  Additional workspace of the appropriate size
Liam Healy's avatar
Liam Healy committed
194 195 196 197 198 199 200 201 202 203
  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.")
204

205 206 207 208 209 210 211 212
(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)))
Liam Healy's avatar
Liam Healy committed
213
     shur-vectors
214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
     &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)))
  :gsl-version (1 9)
  :switch (shur-vectors)
  :inputs (A)
  :outputs (A eigenvalues eigenvectors)
Liam Healy's avatar
Liam Healy committed
230
  :return (eigenvalues eigenvectors)
231 232 233 234 235 236 237 238 239 240
  :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
Liam Healy's avatar
Liam Healy committed
241
  eigenvectors are computed, and an error code is returned.")
242

243 244 245 246 247
;;;;****************************************************************************
;;;; Sorting Eigenvalues and Eigenvectors
;;;;****************************************************************************

(cffi:defcenum eigen-sort-type
Liam Healy's avatar
Liam Healy committed
248
  ;; gsl_eigen_sort_t from /usr/include/gsl/gsl_eigen.h.
249 250
  :value-ascending :value-descending :absolute-ascending :absolute-descending)

Liam Healy's avatar
Liam Healy committed
251 252 253 254 255 256 257 258 259 260 261
(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
262
   "Simultaneously sort the eigenvalues stored in the vector
Liam Healy's avatar
Liam Healy committed
263 264
  eigenvalues and the corresponding real eigenvectors stored in the columns
  of the matrix eigenvectors into ascending or descending order according to
265
  the value of the parameter sort-type: :value-ascending,
Liam Healy's avatar
Liam Healy committed
266
  :value-descending, :absolute-ascending, :absolute-descending.")
267 268 269 270 271

;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************
 
272 273 274 275 276 277 278 279 280 281 282 283
#|
;; When GSL 1.12 is in place and elt+ is supported for complex, this
;; can be expanded to :double-types, if we can create a hermitian
;; matrix.

(generate-all-array-tests eigensystems5 (double-float)
 (let ((m1 (array-default '(5 5))))
   (multiple-value-bind (eval evec)
       (eigenvalues-eigenvectors (elt+ m1 (matrix-transpose-copy m1)))
     (list (cl-array eval) (cl-array evec)))))
|#

284
(defun eigenvalue-eigenvectors-example ()
285
  (let ((evecs (make-marray 'double-float :dimensions '(3 3)))
286 287 288 289
	(evals (make-marray 'double-float :dimensions 3))
	(mat #m((20.0d0 -10.0d0 0.0d0)
		(-10.0d0 30.0d0 0.0d0)
		(0.0d0 0.0d0 40.0d0))))
290
    (eigenvalues-eigenvectors mat evals evecs)
Liam Healy's avatar
Liam Healy committed
291
    (values (cl-array evals) (cl-array evecs))))
292

Liam Healy's avatar
Liam Healy committed
293 294 295
(save-test eigensystems
	   (eigenvalue-eigenvectors-example))