eigensystems.lisp 11.6 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-25 10:37:09EST 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
(defmobject eigen-symm
    "gsl_eigen_symm" ((n sizet))
Liam Healy's avatar
Liam Healy committed
20 21
    "symmetric eigenvalue workspace"
    :documentation			; FDL
22 23 24 25
    "Make a workspace for computing eigenvalues of
  n-by-n real symmetric matrices.  The size of the workspace
  is O(2n).")

26 27
(defmobject eigen-nonsymm
    "gsl_eigen_nonsymm" ((n sizet))
Liam Healy's avatar
Liam Healy committed
28 29 30
    "non-symmetric eigenvalue workspace"
    :gsl-version (1 9)
    :documentation			; FDL
31 32 33 34 35 36 37 38 39 40 41 42 43
    "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))

44 45
(defmobject eigen-symmv
    "gsl_eigen_symmv" ((n sizet))
Liam Healy's avatar
Liam Healy committed
46 47
    "symmetric eigensystem workspace"
    :documentation			; FDL
48 49 50 51
    "Make a workspace for computing eigenvalues and
  eigenvectors of n-by-n real symmetric matrices.  The size of
  the workspace is O(4n).")

52 53
(defmobject eigen-nonsymmv
    "gsl_eigen_nonsymmv" ((n sizet))
Liam Healy's avatar
Liam Healy committed
54 55 56
    "non-symmetric eigenvalue workspace"
    :gsl-version (1 9)
    :documentation			; FDL
57 58 59 60
    "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).")

61 62 63
(defmobject eigen-herm
    "gsl_eigen_herm" ((n sizet))
    "Hermitian eigenvalue workspace"	; FDL
Liam Healy's avatar
Liam Healy committed
64
    :documentation			; FDL
65 66 67
    "Make a workspace for computing eigenvalues of
  n-by-n complex Hermitian matrices.  The size of the workspace
  is O(3n).")
68

69 70
(defmobject eigen-hermv
    "gsl_eigen_hermv" ((n sizet))
Liam Healy's avatar
Liam Healy committed
71 72
    "Hermitian eigensystem workspace"
    :documentation			; FDL
73 74 75
    "Make a workspace for computing eigenvalues and
  eigenvectors of n-by-n complex hermitian matrices.  The size of
  the workspace is O(5n).")
76

Liam Healy's avatar
Liam Healy committed
77 78 79 80
;;;;****************************************************************************
;;;; Eigenvalues and eigenvectors
;;;;****************************************************************************

81 82 83
(defmfun eigenvalues
    ((A matrix)
     &optional
84 85 86
     (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
87
  (double-float "gsl_eigen_symm"
88
		complex-double-float "gsl_eigen_herm")
Liam Healy's avatar
Liam Healy committed
89
  (((mpointer A) :pointer)
90
   ((mpointer eigenvalues) :pointer) ((mpointer ws) :pointer))
Liam Healy's avatar
Liam Healy committed
91 92 93 94 95
  :definition :generic
  :element-types :doubles
  :inputs (A)
  :outputs (A eigenvalues)
  :return (eigenvalues)
96
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
97 98 99 100 101 102 103 104 105
  "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.")

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 175 176 177 178 179 180
(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.")

181 182 183
(defmfun eigenvalues-eigenvectors
    ((A matrix)
     &optional
184 185 186 187
     (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
188
  (double-float "gsl_eigen_symmv"
189
		complex-double-float "gsl_eigen_hermv")  
Liam Healy's avatar
Liam Healy committed
190
  (((mpointer A) :pointer) ((mpointer eigenvalues) :pointer)
191
   ((mpointer eigenvectors) :pointer) ((mpointer ws) :pointer))
Liam Healy's avatar
Liam Healy committed
192 193 194 195 196 197 198
  :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
199
  hermitian matrix A.  Additional workspace of the appropriate size
Liam Healy's avatar
Liam Healy committed
200 201 202 203 204 205 206 207 208 209
  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.")
210

211 212 213 214 215 216 217 218
(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
219
     shur-vectors
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235
     &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
236
  :return (eigenvalues eigenvectors)
237 238 239 240 241 242 243 244 245 246
  :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
247
  eigenvectors are computed, and an error code is returned.")
248

249 250 251 252 253
;;;;****************************************************************************
;;;; Sorting Eigenvalues and Eigenvectors
;;;;****************************************************************************

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

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

;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************
 
278 279 280 281 282 283 284 285 286 287 288 289
#|
;; 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)))))
|#

290
(defun eigenvalue-eigenvectors-example ()
291
  (let ((evecs (make-marray 'double-float :dimensions '(3 3)))
292 293 294 295
	(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))))
296
    (eigenvalues-eigenvectors mat evals evecs)
Liam Healy's avatar
Liam Healy committed
297
    (values (cl-array evals) (cl-array evecs))))
298

Liam Healy's avatar
Liam Healy committed
299 300 301
(save-test eigensystems
	   (eigenvalue-eigenvectors-example))