eigensystems.lisp 11.4 KB
 liam committed Jan 21, 2008 1 2 ``````;; Eigenvectors and eigenvalues ;; Liam Healy, Sun May 21 2006 - 19:52 `````` Liam Healy committed Jan 12, 2009 3 ``````;; Time-stamp: <2009-01-12 10:22:00EST eigensystems.lisp> `````` lhealy committed Jul 24, 2008 4 ``````;; \$Id\$ `````` liam committed May 22, 2006 5 6 7 `````` (in-package :gsl) `````` 8 9 ``````;;; /usr/include/gsl/gsl_eigen.h `````` Liam Healy committed Jan 11, 2009 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? `````` liam committed May 22, 2006 14 ``````;;;;**************************************************************************** `````` Liam Healy committed Aug 14, 2008 15 ``````;;;; Workspace `````` liam committed May 22, 2006 16 17 ``````;;;;**************************************************************************** `````` Liam Healy committed Dec 25, 2008 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).") `````` Liam Healy committed Jan 11, 2009 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)) `````` Liam Healy committed Dec 25, 2008 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).") `````` Liam Healy committed Jan 11, 2009 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).") `````` Liam Healy committed Dec 25, 2008 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).") `````` liam committed May 22, 2006 63 `````` `````` Liam Healy committed Dec 25, 2008 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).") `````` liam committed May 22, 2006 70 `````` `````` Liam Healy committed Aug 14, 2008 71 72 73 74 ``````;;;;**************************************************************************** ;;;; Eigenvalues and eigenvectors ;;;;**************************************************************************** `````` Liam Healy committed Jan 04, 2009 75 76 77 ``````(defmfun eigenvalues ((A matrix) &optional `````` Liam Healy committed Jan 08, 2009 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 committed Aug 14, 2008 81 `````` (double-float "gsl_eigen_symm" `````` Liam Healy committed Jan 04, 2009 82 `````` complex-double-float "gsl_eigen_herm") `````` Liam Healy committed Aug 14, 2008 83 `````` (((mpointer A) :pointer) `````` Liam Healy committed Dec 26, 2008 84 `````` ((mpointer eigenvalues) :pointer) ((mpointer ws) :pointer)) `````` Liam Healy committed Aug 14, 2008 85 86 87 88 89 `````` :definition :generic :element-types :doubles :inputs (A) :outputs (A eigenvalues) :return (eigenvalues) `````` liam committed Jan 21, 2008 90 `````` :documentation ; FDL `````` Liam Healy committed Aug 14, 2008 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.") `````` Liam Healy committed Jan 11, 2009 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.") `````` Liam Healy committed Jan 04, 2009 175 176 177 ``````(defmfun eigenvalues-eigenvectors ((A matrix) &optional `````` Liam Healy committed Jan 08, 2009 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 committed Aug 14, 2008 182 `````` (double-float "gsl_eigen_symmv" `````` Liam Healy committed Jan 04, 2009 183 `````` complex-double-float "gsl_eigen_hermv") `````` Liam Healy committed Aug 14, 2008 184 `````` (((mpointer A) :pointer) ((mpointer eigenvalues) :pointer) `````` Liam Healy committed Dec 26, 2008 185 `````` ((mpointer eigenvectors) :pointer) ((mpointer ws) :pointer)) `````` Liam Healy committed Aug 14, 2008 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 `````` liam committed Jan 21, 2008 193 `````` hermitian matrix A. Additional workspace of the appropriate size `````` Liam Healy committed Aug 14, 2008 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.") `````` liam committed May 22, 2006 204 `````` `````` Liam Healy committed Jan 11, 2009 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 committed Jan 12, 2009 213 `````` shur-vectors `````` Liam Healy committed Jan 11, 2009 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 committed Jan 12, 2009 230 `````` :return (eigenvalues eigenvectors) `````` Liam Healy committed Jan 11, 2009 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 committed Jan 12, 2009 241 `````` eigenvectors are computed, and an error code is returned.") `````` Liam Healy committed Jan 11, 2009 242 `````` `````` liam committed May 22, 2006 243 244 245 246 247 ``````;;;;**************************************************************************** ;;;; Sorting Eigenvalues and Eigenvectors ;;;;**************************************************************************** (cffi:defcenum eigen-sort-type `````` Liam Healy committed Aug 14, 2008 248 `````` ;; gsl_eigen_sort_t from /usr/include/gsl/gsl_eigen.h. `````` liam committed May 22, 2006 249 250 `````` :value-ascending :value-descending :absolute-ascending :absolute-descending) `````` Liam Healy committed Aug 14, 2008 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 `````` liam committed Jan 21, 2008 262 `````` "Simultaneously sort the eigenvalues stored in the vector `````` Liam Healy committed Aug 14, 2008 263 264 `````` eigenvalues and the corresponding real eigenvectors stored in the columns of the matrix eigenvectors into ascending or descending order according to `````` liam committed Jan 21, 2008 265 `````` the value of the parameter sort-type: :value-ascending, `````` Liam Healy committed Aug 14, 2008 266 `````` :value-descending, :absolute-ascending, :absolute-descending.") `````` liam committed May 22, 2006 267 268 269 270 271 `````` ;;;;**************************************************************************** ;;;; Example ;;;;**************************************************************************** `````` Liam Healy committed Jan 05, 2009 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))))) |# `````` liam committed Jan 28, 2008 284 ``````(defun eigenvalue-eigenvectors-example () `````` Liam Healy committed Dec 26, 2008 285 `````` (let ((evecs (make-marray 'double-float :dimensions '(3 3))) `````` Liam Healy committed Jan 04, 2009 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)))) `````` Liam Healy committed Jan 05, 2009 290 `````` (eigenvalues-eigenvectors mat evals evecs) `````` Liam Healy committed Aug 14, 2008 291 `````` (values (cl-array evals) (cl-array evecs)))) `````` liam committed Jan 28, 2008 292 `````` `````` Liam Healy committed Oct 25, 2008 293 294 295 ``````(save-test eigensystems (eigenvalue-eigenvectors-example)) ``````