nonsymmetric-generalized.lisp 7.16 KB
Newer Older
1 2
;; Generalized eigensystems for nonsymmetric real matrices
;; Liam Healy 2009-02-16 14:27:20EST nonsymmetric-generalized.lisp
3
;; Time-stamp: <2012-01-13 12:01:33EST nonsymmetric-generalized.lisp>
4
;;
5
;; Copyright 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21 22 23

(in-package :gsl)

(defmobject eigen-gen
24
    "gsl_eigen_gen" ((n :sizet))
25 26 27 28 29 30 31 32
    "generalized nonsymmetric eigenvalue workspace"
    :gsl-version (1 10)
    :documentation			; FDL
    "Make a workspace for computing eigenvalues of n-by-n real
    generalized nonsymmetric eigensystems. The size of the workspace
    is O(n).")

(defmobject eigen-genv
33
    "gsl_eigen_genv" ((n :sizet))
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
    "generalized nonsymmetric eigenvector and eigenvalue workspace"
    :gsl-version (1 10)
    :documentation			; FDL
    "Make a workspace for computin geigenvalues and eigenvectors of
    n-by-n real generalized nonsymmetric eigensystems. The size of the
    workspace is O(7n).")

(defmfun set-parameters-gen
    (ws &optional compute-shur-form-s compute-shur-form-t balance)
  "gsl_eigen_gen_params"
  (((if compute-shur-form-s 1 0) :int)
   ((if compute-shur-form-t 1 0) :int)
   ((if balance 1 0) :int)		; currently ignored
   ((mpointer ws) :pointer))
  :gsl-version (1 10)
  :c-return :void
  :export nil
  :index eigenvalues-gen)

(defmfun eigenvalues-gen
    (A B
       &optional
       (alpha
57
	(grid:make-foreign-array '(complex double-float) :dimensions (dim0 A)))
58
       (beta
59
	(grid:make-foreign-array 'double-float :dimensions (dim0 A)))
60 61 62 63 64 65
       (ws (make-eigen-gen (dim0 A)))
       compute-shur-form-s compute-shur-form-t shur-vectors
       &aux
       (balance nil)			; unused by GSL
       (Q
	(if (eql shur-vectors t)
66
	    (grid:make-foreign-array 'double-float :dimensions (grid:dimensions A))
67 68 69
	    shur-vectors))
       (Z
	(if (eql shur-vectors t)
70
	    (grid:make-foreign-array 'double-float :dimensions (grid:dimensions A))
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
	    shur-vectors)))
  ("gsl_eigen_gen" "gsl_eigen_gen_QZ")
  ((((mpointer A) :pointer) ((mpointer B) :pointer)
    ((mpointer alpha) :pointer) ((mpointer beta) :pointer) ((mpointer ws) :pointer))
   (((mpointer A) :pointer) ((mpointer B) :pointer)
    ((mpointer alpha) :pointer) ((mpointer beta) :pointer)
    ((mpointer Q) :pointer) ((mpointer Z) :pointer)
    ((mpointer ws) :pointer)))
  :before
  ((set-parameters-gen ws compute-shur-form-s compute-shur-form-t balance))
  :gsl-version (1 10)
  :switch (shur-vectors)
  :inputs (A B)
  :outputs (A B alpha beta)
  :return (alpha beta)
  :documentation			; FDL
  "Compute the eigenvalues of the real generalized nonsymmetric matrix
   pair (A, B), and store them as pairs in (alpha, beta), where alpha
   is complex and beta is real. If \beta_i is non-zero, then \lambda =
   \alpha_i / \beta_i is an eigenvalue. Likewise, if \alpha_i is
   non-zero, then \mu = \beta_i / \alpha_i is an eigenvalue of the
   alternate problem \mu A y = B y. The elements of beta are
   normalized to be non-negative.

   If S is desired, it is stored in A on output. If T is desired, it
   is stored in B on output. The ordering of eigenvalues in (alpha,
   beta) follows the ordering of the diagonal blocks in the Schur
   forms S and T. In rare cases, this function may fail to find all
   eigenvalues. If this occurs, an error code is returned.

   If compute-shur-form-s is true, the full Schur form S will be
   computed. If it is NIL, S will not be computed (this is the default
   setting). S is a quasi upper triangular matrix with 1-by-1 and
   2-by-2 blocks on its diagonal. 1-by-1 blocks correspond to real
   eigenvalues, and 2-by-2 blocks correspond to complex eigenvalues.

   If compute-shur-form-t true, the full Schur form T will be
   computed.  If it is NIL, T will not be
   computed (this is the default setting). T is an upper triangular
   matrix with non-negative elements on its diagonal. Any 2-by-2
   blocks in S will correspond to a 2-by-2 diagonal block in T.")

(defmfun eigenvalues-eigenvectors-gen
    (A B
       &optional
       (alpha
117
	(grid:make-foreign-array '(complex double-float) :dimensions (dim0 A)))
118
       (beta
119
	(grid:make-foreign-array 'double-float :dimensions (dim0 A)))
120
       (eigenvectors
121
	(grid:make-foreign-array 'double-float :dimensions (grid:dimensions A)))
122 123 124 125 126 127
       (ws (make-eigen-genv (dim0 A)))
       compute-shur-form-s compute-shur-form-t shur-vectors
       &aux
       (balance nil)			; unused by GSL
       (Q
	(if (eql shur-vectors t)
128
	    (grid:make-foreign-array 'double-float :dimensions (grid:dimensions A))
129 130 131
	    shur-vectors))
       (Z
	(if (eql shur-vectors t)
132
	    (grid:make-foreign-array 'double-float :dimensions (grid:dimensions A))
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
	    shur-vectors)))
  ("gsl_eigen_genv" "gsl_eigen_genv_QZ")
  ((((mpointer A) :pointer) ((mpointer B) :pointer)
    ((mpointer alpha) :pointer) ((mpointer beta) :pointer)
    ((mpointer eigenvectors) :pointer)
    ((mpointer ws) :pointer))
   (((mpointer A) :pointer) ((mpointer B) :pointer)
    ((mpointer alpha) :pointer) ((mpointer beta) :pointer)
    ((mpointer eigenvectors) :pointer)
    ((mpointer Q) :pointer) ((mpointer Z) :pointer)
    ((mpointer ws) :pointer)))
  :before
  ((set-parameters-gen ws compute-shur-form-s compute-shur-form-t balance))
  :gsl-version (1 10)
  :switch (shur-vectors)
  :inputs (A B)
  :outputs (A B alpha beta eigenvectors)
  :return (alpha beta eigenvectors)
  :documentation			; FDL
  "Compute eigenvalues and right eigenvectors of the n-by-n real
  generalized nonsymmetric matrix pair (A, B). The eigenvalues are
  stored in (alpha, beta) and the eigenvectors are stored in evec. It
  first calls eigenvalues-gen to compute the eigenvalues, Schur forms,
  and Schur vectors. Then it finds eigenvectors of the Schur forms and
  backtransforms them using the Schur vectors. The Schur vectors are
  destroyed in the process, but can be saved by using setting
  shur-vectors true. The computed eigenvectors are normalized to have
  unit magnitude. On output, (A, B) contains the generalized Schur
  form (S, T). 

   If compute-shur-form-s is true, the full Schur form S will be
   computed. If it is NIL, S will not be computed (this is the default
   setting). S is a quasi upper triangular matrix with 1-by-1 and
   2-by-2 blocks on its diagonal. 1-by-1 blocks correspond to real
   eigenvalues, and 2-by-2 blocks correspond to complex eigenvalues.

   If compute-shur-form-t true, the full Schur form T will be
   computed.  If it is NIL, T will not be
   computed (this is the default setting). T is an upper triangular
   matrix with non-negative elements on its diagonal. Any 2-by-2
   blocks in S will correspond to a 2-by-2 diagonal block in T.")