htswanalysis/MACS/lib/gsl/gsl-1.11/doc/eigen.texi

   1 @cindex eigenvalues and eigenvectors
   2 This chapter describes functions for computing eigenvalues and
   3 eigenvectors of matrices.  There are routines for real symmetric,
   4 real nonsymmetric, complex hermitian, real generalized symmetric-definite,
   5 complex generalized hermitian-definite, and real generalized nonsymmetric
   6 eigensystems. Eigenvalues can be computed with or without eigenvectors.
   7 The hermitian and real symmetric matrix algorithms are symmetric bidiagonalization
   8 followed by QR reduction. The nonsymmetric algorithm is the Francis QR
   9 double-shift.  The generalized nonsymmetric algorithm is the QZ method due
  10 to Moler and Stewart.
  11
  12 The functions described in this chapter are declared in the header file
  13 @file{gsl_eigen.h}.
  14
  15 @menu
  16 * Real Symmetric Matrices::
  17 * Complex Hermitian Matrices::
  18 * Real Nonsymmetric Matrices::
  19 * Real Generalized Symmetric-Definite Eigensystems::
  20 * Complex Generalized Hermitian-Definite Eigensystems::
  21 * Real Generalized Nonsymmetric Eigensystems::
  22 * Sorting Eigenvalues and Eigenvectors::
  23 * Eigenvalue and Eigenvector Examples::
  24 * Eigenvalue and Eigenvector References::
  25 @end menu
  26
  27 @node Real Symmetric Matrices
  28 @section Real Symmetric Matrices
  29 @cindex symmetric matrix, real, eigensystem
  30 @cindex real symmetric matrix, eigensystem
  31
  32 For real symmetric matrices, the library uses the symmetric
  33 bidiagonalization and QR reduction method.  This is described in Golub
  34 & van Loan, section 8.3.  The computed eigenvalues are accurate to an
  35 absolute accuracy of @math{\epsilon ||A||_2}, where @math{\epsilon} is
  36 the machine precision.
  37
  38 @deftypefun {gsl_eigen_symm_workspace *} gsl_eigen_symm_alloc (const size_t @var{n})
  39 This function allocates a workspace for computing eigenvalues of
  40 @var{n}-by-@var{n} real symmetric matrices.  The size of the workspace
  41 is @math{O(2n)}.
  42 @end deftypefun
  43
  44 @deftypefun void gsl_eigen_symm_free (gsl_eigen_symm_workspace * @var{w})
  45 This function frees the memory associated with the workspace @var{w}.
  46 @end deftypefun
  47
  48 @deftypefun int gsl_eigen_symm (gsl_matrix * @var{A}, gsl_vector * @var{eval}, gsl_eigen_symm_workspace * @var{w})
  49 This function computes the eigenvalues of the real symmetric matrix
  50 @var{A}.  Additional workspace of the appropriate size must be provided
  51 in @var{w}.  The diagonal and lower triangular part of @var{A} are
  52 destroyed during the computation, but the strict upper triangular part
  53 is not referenced.  The eigenvalues are stored in the vector @var{eval}
  54 and are unordered.
  55 @end deftypefun
  56
  57 @deftypefun {gsl_eigen_symmv_workspace *} gsl_eigen_symmv_alloc (const size_t @var{n})
  58 This function allocates a workspace for computing eigenvalues and
  59 eigenvectors of @var{n}-by-@var{n} real symmetric matrices.  The size of
  60 the workspace is @math{O(4n)}.
  61 @end deftypefun
  62
  63 @deftypefun void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * @var{w})
  64 This function frees the memory associated with the workspace @var{w}.
  65 @end deftypefun
  66
  67 @deftypefun int gsl_eigen_symmv (gsl_matrix * @var{A}, gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_symmv_workspace * @var{w})
  68 This function computes the eigenvalues and eigenvectors of the real
  69 symmetric matrix @var{A}.  Additional workspace of the appropriate size
  70 must be provided in @var{w}.  The diagonal and lower triangular part of
  71 @var{A} are destroyed during the computation, but the strict upper
  72 triangular part is not referenced.  The eigenvalues are stored in the
  73 vector @var{eval} and are unordered.  The corresponding eigenvectors are
  74 stored in the columns of the matrix @var{evec}.  For example, the
  75 eigenvector in the first column corresponds to the first eigenvalue.
  76 The eigenvectors are guaranteed to be mutually orthogonal and normalised
  77 to unit magnitude.
  78 @end deftypefun
  79
  80 @node Complex Hermitian Matrices
  81 @section Complex Hermitian Matrices
  82
  83 @cindex hermitian matrix, complex, eigensystem
  84 @cindex complex hermitian matrix, eigensystem
  85
  86 @deftypefun {gsl_eigen_herm_workspace *} gsl_eigen_herm_alloc (const size_t @var{n})
  87 This function allocates a workspace for computing eigenvalues of
  88 @var{n}-by-@var{n} complex hermitian matrices.  The size of the workspace
  89 is @math{O(3n)}.
  90 @end deftypefun
  91
  92 @deftypefun void gsl_eigen_herm_free (gsl_eigen_herm_workspace * @var{w})
  93 This function frees the memory associated with the workspace @var{w}.
  94 @end deftypefun
  95
  96 @deftypefun int gsl_eigen_herm (gsl_matrix_complex * @var{A}, gsl_vector * @var{eval}, gsl_eigen_herm_workspace * @var{w})
  97 This function computes the eigenvalues of the complex hermitian matrix
  98 @var{A}.  Additional workspace of the appropriate size must be provided
  99 in @var{w}.  The diagonal and lower triangular part of @var{A} are
 100 destroyed during the computation, but the strict upper triangular part
 101 is not referenced.  The imaginary parts of the diagonal are assumed to be
 102 zero and are not referenced. The eigenvalues are stored in the vector
 103 @var{eval} and are unordered.
 104 @end deftypefun
 105
 106 @deftypefun {gsl_eigen_hermv_workspace *} gsl_eigen_hermv_alloc (const size_t @var{n})
 107 This function allocates a workspace for computing eigenvalues and
 108 eigenvectors of @var{n}-by-@var{n} complex hermitian matrices.  The size of
 109 the workspace is @math{O(5n)}.
 110 @end deftypefun
 111
 112 @deftypefun void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * @var{w})
 113 This function frees the memory associated with the workspace @var{w}.
 114 @end deftypefun
 115
 116 @deftypefun int gsl_eigen_hermv (gsl_matrix_complex * @var{A}, gsl_vector * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_hermv_workspace * @var{w})
 117 This function computes the eigenvalues and eigenvectors of the complex
 118 hermitian matrix @var{A}.  Additional workspace of the appropriate size
 119 must be provided in @var{w}.  The diagonal and lower triangular part of
 120 @var{A} are destroyed during the computation, but the strict upper
 121 triangular part is not referenced. The imaginary parts of the diagonal
 122 are assumed to be zero and are not referenced.  The eigenvalues are
 123 stored in the vector @var{eval} and are unordered.  The corresponding
 124 complex eigenvectors are stored in the columns of the matrix @var{evec}.
 125 For example, the eigenvector in the first column corresponds to the
 126 first eigenvalue.  The eigenvectors are guaranteed to be mutually
 127 orthogonal and normalised to unit magnitude.
 128 @end deftypefun
 129
 130 @node Real Nonsymmetric Matrices
 131 @section Real Nonsymmetric Matrices
 132 @cindex nonsymmetric matrix, real, eigensystem
 133 @cindex real nonsymmetric matrix, eigensystem
 134
 135 The solution of the real nonsymmetric eigensystem problem for a
 136 matrix @math{A} involves computing the Schur decomposition
 137 @tex
 138 \beforedisplay
 139 $$
 140 A = Z T Z^T
 141 $$
 142 \afterdisplay
 143 @end tex
 144 @ifinfo
 145
 146 @example
 147 A = Z T Z^T
 148 @end example
 149
 150 @end ifinfo
 151 where @math{Z} is an orthogonal matrix of Schur vectors and @math{T},
 152 the Schur form, is quasi upper triangular with diagonal
 153 @math{1}-by-@math{1} blocks which are real eigenvalues of @math{A}, and
 154 diagonal @math{2}-by-@math{2} blocks whose eigenvalues are complex
 155 conjugate eigenvalues of @math{A}. The algorithm used is the double
 156 shift Francis method.
 157
 158 @deftypefun {gsl_eigen_nonsymm_workspace *} gsl_eigen_nonsymm_alloc (const size_t @var{n})
 159 This function allocates a workspace for computing eigenvalues of
 160 @var{n}-by-@var{n} real nonsymmetric matrices. The size of the workspace
 161 is @math{O(2n)}.
 162 @end deftypefun
 163
 164 @deftypefun void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * @var{w})
 165 This function frees the memory associated with the workspace @var{w}.
 166 @end deftypefun
 167
 168 @deftypefun void gsl_eigen_nonsymm_params (const int @var{compute_t}, const int @var{balance}, gsl_eigen_nonsymm_workspace * @var{w})
 169 This function sets some parameters which determine how the eigenvalue
 170 problem is solved in subsequent calls to @code{gsl_eigen_nonsymm}.
 171
 172 If @var{compute_t} is set to 1, the full Schur form @math{T} will be
 173 computed by @code{gsl_eigen_nonsymm}. If it is set to 0,
 174 @math{T} will not be computed (this is the default setting). Computing
 175 the full Schur form @math{T} requires approximately 1.5-2 times the
 176 number of flops.
 177
 178 If @var{balance} is set to 1, a balancing transformation is applied
 179 to the matrix prior to computing eigenvalues. This transformation is
 180 designed to make the rows and columns of the matrix have comparable
 181 norms, and can result in more accurate eigenvalues for matrices
 182 whose entries vary widely in magnitude. See @ref{Balancing} for more
 183 information. Note that the balancing transformation does not preserve
 184 the orthogonality of the Schur vectors, so if you wish to compute the
 185 Schur vectors with @code{gsl_eigen_nonsymm_Z} you will obtain the Schur
 186 vectors of the balanced matrix instead of the original matrix. The
 187 relationship will be
 188 @tex
 189 \beforedisplay
 190 $$
 191 T = Q^t D^{-1} A D Q
 192 $$
 193 \afterdisplay
 194 @end tex
 195 @ifinfo
 196
 197 @example
 198 T = Q^t D^(-1) A D Q
 199 @end example
 200
 201 @end ifinfo
 202 @noindent
 203 where @var{Q} is the matrix of Schur vectors for the balanced matrix, and
 204 @var{D} is the balancing transformation. Then @code{gsl_eigen_nonsymm_Z}
 205 will compute a matrix @var{Z} which satisfies
 206 @tex
 207 \beforedisplay
 208 $$
 209 T = Z^{-1} A Z
 210 $$
 211 \afterdisplay
 212 @end tex
 213 @ifinfo
 214
 215 @example
 216 T = Z^(-1) A Z
 217 @end example
 218
 219 @end ifinfo
 220 @noindent
 221 with @math{Z = D Q}. Note that @var{Z} will not be orthogonal. For
 222 this reason, balancing is not performed by default.
 223 @end deftypefun
 224
 225 @deftypefun int gsl_eigen_nonsymm (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_eigen_nonsymm_workspace * @var{w})
 226 This function computes the eigenvalues of the real nonsymmetric matrix
 227 @var{A} and stores them in the vector @var{eval}. If @math{T} is
 228 desired, it is stored in the upper portion of @var{A} on output.
 229 Otherwise, on output, the diagonal of @var{A} will contain the
 230 @math{1}-by-@math{1} real eigenvalues and @math{2}-by-@math{2}
 231 complex conjugate eigenvalue systems, and the rest of @var{A} is
 232 destroyed. In rare cases, this function may fail to find all
 233 eigenvalues. If this happens, an error code is returned
 234 and the number of converged eigenvalues is stored in @code{w->n_evals}.
 235 The converged eigenvalues are stored in the beginning of @var{eval}.
 236 @end deftypefun
 237
 238 @deftypefun int gsl_eigen_nonsymm_Z (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix * @var{Z}, gsl_eigen_nonsymm_workspace * @var{w})
 239 This function is identical to @code{gsl_eigen_nonsymm} except it also
 240 computes the Schur vectors and stores them into @var{Z}.
 241 @end deftypefun
 242
 243 @deftypefun {gsl_eigen_nonsymmv_workspace *} gsl_eigen_nonsymmv_alloc (const size_t @var{n})
 244 This function allocates a workspace for computing eigenvalues and
 245 eigenvectors of @var{n}-by-@var{n} real nonsymmetric matrices. The
 246 size of the workspace is @math{O(5n)}.
 247 @end deftypefun
 248
 249 @deftypefun void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * @var{w})
 250 This function frees the memory associated with the workspace @var{w}.
 251 @end deftypefun
 252
 253 @deftypefun int gsl_eigen_nonsymmv (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_nonsymmv_workspace * @var{w})
 254 This function computes eigenvalues and right eigenvectors of the
 255 @var{n}-by-@var{n} real nonsymmetric matrix @var{A}. It first calls
 256 @code{gsl_eigen_nonsymm} to compute the eigenvalues, Schur form @math{T}, and
 257 Schur vectors. Then it finds eigenvectors of @math{T} and backtransforms
 258 them using the Schur vectors. The Schur vectors are destroyed in the
 259 process, but can be saved by using @code{gsl_eigen_nonsymmv_Z}. The
 260 computed eigenvectors are normalized to have unit magnitude. On
 261 output, the upper portion of @var{A} contains the Schur form
 262 @math{T}. If @code{gsl_eigen_nonsymm} fails, no eigenvectors are
 263 computed, and an error code is returned.
 264 @end deftypefun
 265
 266 @deftypefun int gsl_eigen_nonsymmv_Z (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_matrix * @var{Z}, gsl_eigen_nonsymmv_workspace * @var{w})
 267 This function is identical to @code{gsl_eigen_nonsymmv} except it also saves
 268 the Schur vectors into @var{Z}.
 269 @end deftypefun
 270
 271 @node Real Generalized Symmetric-Definite Eigensystems
 272 @section Real Generalized Symmetric-Definite Eigensystems
 273 @cindex generalized symmetric eigensystems
 274
 275 The real generalized symmetric-definite eigenvalue problem is to find
 276 eigenvalues @math{\lambda} and eigenvectors @math{x} such that
 277 @tex
 278 \beforedisplay
 279 $$
 280 A x = \lambda B x
 281 $$
 282 \afterdisplay
 283 @end tex
 284 @ifinfo
 285 @example
 286 A x = \lambda B x
 287 @end example
 288 @end ifinfo
 289 where @math{A} and @math{B} are symmetric matrices, and @math{B} is
 290 positive-definite. This problem reduces to the standard symmetric
 291 eigenvalue problem by applying the Cholesky decomposition to @math{B}:
 292 @tex
 293 \beforedisplay
 294 $$
 295 \eqalign{
 296 A x & = \lambda B x \cr
 297 A x & = \lambda L L^t x \cr
 298 \left( L^{-1} A L^{-t} \right) L^t x & = \lambda L^t x
 299 }
 300 $$
 301 \afterdisplay
 302 @end tex
 303 @ifinfo
 304 @example
 305                       A x = \lambda B x
 306                       A x = \lambda L L^t x
 307 ( L^@{-1@} A L^@{-t@} ) L^t x = \lambda L^t x
 308 @end example
 309 @end ifinfo
 310 Therefore, the problem becomes @math{C y = \lambda y} where
 311 @c{$C = L^{-1} A L^{-t}$}
 312 @math{C = L^@{-1@} A L^@{-t@}}
 313 is symmetric, and @math{y = L^t x}. The standard
 314 symmetric eigensolver can be applied to the matrix @math{C}.
 315 The resulting eigenvectors are backtransformed to find the
 316 vectors of the original problem. The eigenvalues and eigenvectors
 317 of the generalized symmetric-definite eigenproblem are always real.
 318
 319 @deftypefun {gsl_eigen_gensymm_workspace *} gsl_eigen_gensymm_alloc (const size_t @var{n})
 320 This function allocates a workspace for computing eigenvalues of
 321 @var{n}-by-@var{n} real generalized symmetric-definite eigensystems. The
 322 size of the workspace is @math{O(2n)}.
 323 @end deftypefun
 324
 325 @deftypefun void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * @var{w})
 326 This function frees the memory associated with the workspace @var{w}.
 327 @end deftypefun
 328
 329 @deftypefun int gsl_eigen_gensymm (gsl_matrix * @var{A}, gsl_matrix * @var{B}, gsl_vector * @var{eval}, gsl_eigen_gensymm_workspace * @var{w})
 330 This function computes the eigenvalues of the real generalized
 331 symmetric-definite matrix pair (@var{A}, @var{B}), and stores them
 332 in @var{eval}, using the method outlined above. On output, @var{B}
 333 contains its Cholesky decomposition and @var{A} is destroyed.
 334 @end deftypefun
 335
 336 @deftypefun {gsl_eigen_gensymmv_workspace *} gsl_eigen_gensymmv_alloc (const size_t @var{n})
 337 This function allocates a workspace for computing eigenvalues and
 338 eigenvectors of @var{n}-by-@var{n} real generalized symmetric-definite
 339 eigensystems. The size of the workspace is @math{O(4n)}.
 340 @end deftypefun
 341
 342 @deftypefun void gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * @var{w})
 343 This function frees the memory associated with the workspace @var{w}.
 344 @end deftypefun
 345
 346 @deftypefun int gsl_eigen_gensymmv (gsl_matrix * @var{A}, gsl_matrix * @var{B}, gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_gensymmv_workspace * @var{w})
 347 This function computes the eigenvalues and eigenvectors of the real
 348 generalized symmetric-definite matrix pair (@var{A}, @var{B}), and
 349 stores them in @var{eval} and @var{evec} respectively. The computed
 350 eigenvectors are normalized to have unit magnitude. On output,
 351 @var{B} contains its Cholesky decomposition and @var{A} is destroyed.
 352 @end deftypefun
 353
 354 @node Complex Generalized Hermitian-Definite Eigensystems
 355 @section Complex Generalized Hermitian-Definite Eigensystems
 356 @cindex generalized hermitian definite eigensystems
 357
 358 The complex generalized hermitian-definite eigenvalue problem is to find
 359 eigenvalues @math{\lambda} and eigenvectors @math{x} such that
 360 @tex
 361 \beforedisplay
 362 $$
 363 A x = \lambda B x
 364 $$
 365 \afterdisplay
 366 @end tex
 367 @ifinfo
 368 @example
 369 A x = \lambda B x
 370 @end example
 371 @end ifinfo
 372 where @math{A} and @math{B} are hermitian matrices, and @math{B} is
 373 positive-definite. Similarly to the real case, this can be reduced
 374 to @math{C y = \lambda y} where
 375 @c{$C = L^{-1} A L^{-\dagger}$}
 376 @math{C = L^@{-1@} A L^@{-H@}}
 377 is hermitian, and
 378 @c{$y = L^{\dagger} x$}
 379 @math{y = L^H x}. The standard
 380 hermitian eigensolver can be applied to the matrix @math{C}.
 381 The resulting eigenvectors are backtransformed to find the
 382 vectors of the original problem. The eigenvalues
 383 of the generalized hermitian-definite eigenproblem are always real.
 384
 385 @deftypefun {gsl_eigen_genherm_workspace *} gsl_eigen_genherm_alloc (const size_t @var{n})
 386 This function allocates a workspace for computing eigenvalues of
 387 @var{n}-by-@var{n} complex generalized hermitian-definite eigensystems. The
 388 size of the workspace is @math{O(3n)}.
 389 @end deftypefun
 390
 391 @deftypefun void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * @var{w})
 392 This function frees the memory associated with the workspace @var{w}.
 393 @end deftypefun
 394
 395 @deftypefun int gsl_eigen_genherm (gsl_matrix_complex * @var{A}, gsl_matrix_complex * @var{B}, gsl_vector * @var{eval}, gsl_eigen_genherm_workspace * @var{w})
 396 This function computes the eigenvalues of the complex generalized
 397 hermitian-definite matrix pair (@var{A}, @var{B}), and stores them
 398 in @var{eval}, using the method outlined above. On output, @var{B}
 399 contains its Cholesky decomposition and @var{A} is destroyed.
 400 @end deftypefun
 401
 402 @deftypefun {gsl_eigen_genhermv_workspace *} gsl_eigen_genhermv_alloc (const size_t @var{n})
 403 This function allocates a workspace for computing eigenvalues and
 404 eigenvectors of @var{n}-by-@var{n} complex generalized hermitian-definite
 405 eigensystems. The size of the workspace is @math{O(5n)}.
 406 @end deftypefun
 407
 408 @deftypefun void gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * @var{w})
 409 This function frees the memory associated with the workspace @var{w}.
 410 @end deftypefun
 411
 412 @deftypefun int gsl_eigen_genhermv (gsl_matrix_complex * @var{A}, gsl_matrix_complex * @var{B}, gsl_vector * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_genhermv_workspace * @var{w})
 413 This function computes the eigenvalues and eigenvectors of the complex
 414 generalized hermitian-definite matrix pair (@var{A}, @var{B}), and
 415 stores them in @var{eval} and @var{evec} respectively. The computed
 416 eigenvectors are normalized to have unit magnitude. On output,
 417 @var{B} contains its Cholesky decomposition and @var{A} is destroyed.
 418 @end deftypefun
 419
 420 @node Real Generalized Nonsymmetric Eigensystems
 421 @section Real Generalized Nonsymmetric Eigensystems
 422 @cindex generalized eigensystems
 423
 424 Given two square matrices (@math{A}, @math{B}), the generalized
 425 nonsymmetric eigenvalue problem is to find eigenvalues @math{\lambda} and
 426 eigenvectors @math{x} such that
 427 @tex
 428 \beforedisplay
 429 $$
 430 A x = \lambda B x
 431 $$
 432 \afterdisplay
 433 @end tex
 434 @ifinfo
 435 @example
 436 A x = \lambda B x
 437 @end example
 438 @end ifinfo
 439 We may also define the problem as finding eigenvalues @math{\mu} and
 440 eigenvectors @math{y} such that
 441 @tex
 442 \beforedisplay
 443 $$
 444 \mu A y = B y
 445 $$
 446 \afterdisplay
 447 @end tex
 448 @ifinfo
 449 @example
 450 \mu A y = B y
 451 @end example
 452 @end ifinfo
 453 Note that these two problems are equivalent (with @math{\lambda = 1/\mu})
 454 if neither @math{\lambda} nor @math{\mu} is zero. If say, @math{\lambda}
 455 is zero, then it is still a well defined eigenproblem, but its alternate
 456 problem involving @math{\mu} is not. Therefore, to allow for zero
 457 (and infinite) eigenvalues, the problem which is actually solved is
 458 @tex
 459 \beforedisplay
 460 $$
 461 \beta A x = \alpha B x
 462 $$
 463 \afterdisplay
 464 @end tex
 465 @ifinfo
 466 @example
 467 \beta A x = \alpha B x
 468 @end example
 469 @end ifinfo
 470 The eigensolver routines below will return two values @math{\alpha}
 471 and @math{\beta} and leave it to the user to perform the divisions
 472 @math{\lambda = \alpha / \beta} and @math{\mu = \beta / \alpha}.
 473
 474 If the determinant of the matrix pencil @math{A - \lambda B} is zero
 475 for all @math{\lambda}, the problem is said to be singular; otherwise
 476 it is called regular.  Singularity normally leads to some
 477 @math{\alpha = \beta = 0} which means the eigenproblem is ill-conditioned
 478 and generally does not have well defined eigenvalue solutions. The
 479 routines below are intended for regular matrix pencils and could yield
 480 unpredictable results when applied to singular pencils.
 481
 482 The solution of the real generalized nonsymmetric eigensystem problem for a
 483 matrix pair @math{(A, B)} involves computing the generalized Schur
 484 decomposition
 485 @tex
 486 \beforedisplay
 487 $$
 488 A = Q S Z^T
 489 $$
 490 $$
 491 B = Q T Z^T
 492 $$
 493 \afterdisplay
 494 @end tex
 495 @ifinfo
 496 @example
 497 A = Q S Z^T
 498 B = Q T Z^T
 499 @end example
 500 @end ifinfo
 501 where @math{Q} and @math{Z} are orthogonal matrices of left and right
 502 Schur vectors respectively, and @math{(S, T)} is the generalized Schur
 503 form whose diagonal elements give the @math{\alpha} and @math{\beta}
 504 values. The algorithm used is the QZ method due to Moler and Stewart
 505 (see references).
 506
 507 @deftypefun {gsl_eigen_gen_workspace *} gsl_eigen_gen_alloc (const size_t @var{n})
 508 This function allocates a workspace for computing eigenvalues of
 509 @var{n}-by-@var{n} real generalized nonsymmetric eigensystems. The
 510 size of the workspace is @math{O(n)}.
 511 @end deftypefun
 512
 513 @deftypefun void gsl_eigen_gen_free (gsl_eigen_gen_workspace * @var{w})
 514 This function frees the memory associated with the workspace @var{w}.
 515 @end deftypefun
 516
 517 @deftypefun void gsl_eigen_gen_params (const int @var{compute_s}, const int @var{compute_t}, const int @var{balance}, gsl_eigen_gen_workspace * @var{w})
 518 This function sets some parameters which determine how the eigenvalue
 519 problem is solved in subsequent calls to @code{gsl_eigen_gen}.
 520
 521 If @var{compute_s} is set to 1, the full Schur form @math{S} will be
 522 computed by @code{gsl_eigen_gen}. If it is set to 0,
 523 @math{S} will not be computed (this is the default setting). @math{S}
 524 is a quasi upper triangular matrix with 1-by-1 and 2-by-2 blocks
 525 on its diagonal. 1-by-1 blocks correspond to real eigenvalues, and
 526 2-by-2 blocks correspond to complex eigenvalues.
 527
 528 If @var{compute_t} is set to 1, the full Schur form @math{T} will be
 529 computed by @code{gsl_eigen_gen}. If it is set to 0,
 530 @math{T} will not be computed (this is the default setting). @math{T}
 531 is an upper triangular matrix with non-negative elements on its diagonal.
 532 Any 2-by-2 blocks in @math{S} will correspond to a 2-by-2 diagonal
 533 block in @math{T}.
 534
 535 The @var{balance} parameter is currently ignored, since generalized
 536 balancing is not yet implemented.
 537 @end deftypefun
 538
 539 @deftypefun int gsl_eigen_gen (gsl_matrix * @var{A}, gsl_matrix * @var{B}, gsl_vector_complex * @var{alpha}, gsl_vector * @var{beta}, gsl_eigen_gen_workspace * @var{w})
 540 This function computes the eigenvalues of the real generalized nonsymmetric
 541 matrix pair (@var{A}, @var{B}), and stores them as pairs in
 542 (@var{alpha}, @var{beta}), where @var{alpha} is complex and @var{beta} is
 543 real. If @math{\beta_i} is non-zero, then
 544 @math{\lambda = \alpha_i / \beta_i} is an eigenvalue. Likewise,
 545 if @math{\alpha_i} is non-zero, then
 546 @math{\mu = \beta_i / \alpha_i} is an eigenvalue of the alternate
 547 problem @math{\mu A y = B y}. The elements of @var{beta} are normalized
 548 to be non-negative.
 549
 550 If @math{S} is desired, it is stored in @var{A} on output. If @math{T}
 551 is desired, it is stored in @var{B} on output. The ordering of
 552 eigenvalues in (@var{alpha}, @var{beta}) follows the ordering
 553 of the diagonal blocks in the Schur forms @math{S} and @math{T}. In rare
 554 cases, this function may fail to find all eigenvalues. If this occurs, an
 555 error code is returned.
 556 @end deftypefun
 557
 558 @deftypefun int gsl_eigen_gen_QZ (gsl_matrix * @var{A}, gsl_matrix * @var{B}, gsl_vector_complex * @var{alpha}, gsl_vector * @var{beta}, gsl_matrix * @var{Q}, gsl_matrix * @var{Z}, gsl_eigen_gen_workspace * @var{w})
 559 This function is identical to @code{gsl_eigen_gen} except it also
 560 computes the left and right Schur vectors and stores them into @var{Q}
 561 and @var{Z} respectively.
 562 @end deftypefun
 563
 564 @deftypefun {gsl_eigen_genv_workspace *} gsl_eigen_genv_alloc (const size_t @var{n})
 565 This function allocates a workspace for computing eigenvalues and
 566 eigenvectors of @var{n}-by-@var{n} real generalized nonsymmetric
 567 eigensystems. The size of the workspace is @math{O(7n)}.
 568 @end deftypefun
 569
 570 @deftypefun void gsl_eigen_genv_free (gsl_eigen_genv_workspace * @var{w})
 571 This function frees the memory associated with the workspace @var{w}.
 572 @end deftypefun
 573
 574 @deftypefun int gsl_eigen_genv (gsl_matrix * @var{A}, gsl_matrix * @var{B}, gsl_vector_complex * @var{alpha}, gsl_vector * @var{beta}, gsl_matrix_complex * @var{evec}, gsl_eigen_genv_workspace * @var{w})
 575 This function computes eigenvalues and right eigenvectors of the
 576 @var{n}-by-@var{n} real generalized nonsymmetric matrix pair
 577 (@var{A}, @var{B}). The eigenvalues are stored in (@var{alpha}, @var{beta})
 578 and the eigenvectors are stored in @var{evec}. It first calls
 579 @code{gsl_eigen_gen} to compute the eigenvalues, Schur forms, and
 580 Schur vectors. Then it finds eigenvectors of the Schur forms and
 581 backtransforms them using the Schur vectors. The Schur vectors are
 582 destroyed in the process, but can be saved by using
 583 @code{gsl_eigen_genv_QZ}. The computed eigenvectors are normalized
 584 to have unit magnitude. On output, (@var{A}, @var{B}) contains
 585 the generalized Schur form (@math{S}, @math{T}). If @code{gsl_eigen_gen}
 586 fails, no eigenvectors are computed, and an error code is returned.
 587 @end deftypefun
 588
 589 @deftypefun int gsl_eigen_genv_QZ (gsl_matrix * @var{A}, gsl_matrix * @var{B}, gsl_vector_complex * @var{alpha}, gsl_vector * @var{beta}, gsl_matrix_complex * @var{evec}, gsl_matrix * @var{Q}, gsl_matrix * @var{Z}, gsl_eigen_genv_workspace * @var{w})
 590 This function is identical to @code{gsl_eigen_genv} except it also
 591 computes the left and right Schur vectors and stores them into @var{Q}
 592 and @var{Z} respectively.
 593 @end deftypefun
 594
 595 @node Sorting Eigenvalues and Eigenvectors
 596 @section Sorting Eigenvalues and Eigenvectors
 597 @cindex sorting eigenvalues and eigenvectors
 598
 599 @deftypefun int gsl_eigen_symmv_sort (gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_sort_t @var{sort_type})
 600 This function simultaneously sorts the eigenvalues stored in the vector
 601 @var{eval} and the corresponding real eigenvectors stored in the columns
 602 of the matrix @var{evec} into ascending or descending order according to
 603 the value of the parameter @var{sort_type},
 604
 605 @table @code
 606 @item GSL_EIGEN_SORT_VAL_ASC
 607 ascending order in numerical value
 608 @item GSL_EIGEN_SORT_VAL_DESC
 609 descending order in numerical value
 610 @item GSL_EIGEN_SORT_ABS_ASC
 611 ascending order in magnitude
 612 @item GSL_EIGEN_SORT_ABS_DESC
 613 descending order in magnitude
 614 @end table
 615
 616 @end deftypefun
 617
 618 @deftypefun int gsl_eigen_hermv_sort (gsl_vector * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_sort_t @var{sort_type})
 619 This function simultaneously sorts the eigenvalues stored in the vector
 620 @var{eval} and the corresponding complex eigenvectors stored in the
 621 columns of the matrix @var{evec} into ascending or descending order
 622 according to the value of the parameter @var{sort_type} as shown above.
 623 @end deftypefun
 624
 625 @deftypefun int gsl_eigen_nonsymmv_sort (gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_sort_t @var{sort_type})
 626 This function simultaneously sorts the eigenvalues stored in the vector
 627 @var{eval} and the corresponding complex eigenvectors stored in the
 628 columns of the matrix @var{evec} into ascending or descending order
 629 according to the value of the parameter @var{sort_type} as shown above.
 630 Only @code{GSL_EIGEN_SORT_ABS_ASC} and @code{GSL_EIGEN_SORT_ABS_DESC} are
 631 supported due to the eigenvalues being complex.
 632 @end deftypefun
 633
 634 @deftypefun int gsl_eigen_gensymmv_sort (gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_sort_t @var{sort_type})
 635 This function simultaneously sorts the eigenvalues stored in the vector
 636 @var{eval} and the corresponding real eigenvectors stored in the columns
 637 of the matrix @var{evec} into ascending or descending order according to
 638 the value of the parameter @var{sort_type} as shown above.
 639 @end deftypefun
 640
 641 @deftypefun int gsl_eigen_genhermv_sort (gsl_vector * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_sort_t @var{sort_type})
 642 This function simultaneously sorts the eigenvalues stored in the vector
 643 @var{eval} and the corresponding complex eigenvectors stored in the columns
 644 of the matrix @var{evec} into ascending or descending order according to
 645 the value of the parameter @var{sort_type} as shown above.
 646 @end deftypefun
 647
 648 @deftypefun int gsl_eigen_genv_sort (gsl_vector_complex * @var{alpha}, gsl_vector * @var{beta}, gsl_matrix_complex * @var{evec}, gsl_eigen_sort_t @var{sort_type})
 649 This function simultaneously sorts the eigenvalues stored in the vectors
 650 (@var{alpha}, @var{beta}) and the corresponding complex eigenvectors
 651 stored in the columns of the matrix @var{evec} into ascending or
 652 descending order according to the value of the parameter @var{sort_type}
 653 as shown above. Only @code{GSL_EIGEN_SORT_ABS_ASC} and
 654 @code{GSL_EIGEN_SORT_ABS_DESC} are supported due to the eigenvalues being
 655 complex.
 656 @end deftypefun
 657
 658 @comment @deftypefun int gsl_eigen_jacobi (gsl_matrix * @var{matrix}, gsl_vector * @var{eval}, gsl_matrix * @var{evec}, unsigned int @var{max_rot}, unsigned int * @var{nrot})
 659 @comment This function finds the eigenvectors and eigenvalues of a real symmetric
 660 @comment matrix by Jacobi iteration. The data in the input matrix is destroyed.
 661 @comment @end deftypefun
 662
 663 @comment @deftypefun int gsl_la_invert_jacobi (const gsl_matrix * @var{matrix}, gsl_matrix * @var{ainv}, unsigned int @var{max_rot})
 664 @comment Invert a matrix by Jacobi iteration.
 665 @comment @end deftypefun
 666
 667 @comment @deftypefun int gsl_eigen_sort (gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_sort_t @var{sort_type})
 668 @comment This functions sorts the eigensystem results based on eigenvalues.
 669 @comment Sorts in order of increasing value or increasing
 670 @comment absolute value, depending on the value of
 671 @comment @var{sort_type}, which can be @code{GSL_EIGEN_SORT_VALUE}
 672 @comment or @code{GSL_EIGEN_SORT_ABSVALUE}.
 673 @comment @end deftypefun
 674
 675 @node Eigenvalue and Eigenvector Examples
 676 @section Examples
 677
 678 The following program computes the eigenvalues and eigenvectors of the 4-th order Hilbert matrix, @math{H(i,j) = 1/(i + j + 1)}.
 679
 680 @example
 681 @verbatiminclude examples/eigen.c
 682 @end example
 683
 684 @noindent
 685 Here is the beginning of the output from the program,
 686
 687 @example
 688 $ ./a.out
 689 eigenvalue = 9.67023e-05
 690 eigenvector =
 691 -0.0291933
 692 0.328712
 693 -0.791411
 694 0.514553
 695 ...
 696 @end example
 697
 698 @noindent
 699 This can be compared with the corresponding output from @sc{gnu octave},
 700
 701 @example
 702 octave> [v,d] = eig(hilb(4));
 703 octave> diag(d)
 704 ans =
 705
 706    9.6702e-05
 707    6.7383e-03
 708    1.6914e-01
 709    1.5002e+00
 710
 711 octave> v
 712 v =
 713
 714    0.029193   0.179186  -0.582076   0.792608
 715   -0.328712  -0.741918   0.370502   0.451923
 716    0.791411   0.100228   0.509579   0.322416
 717   -0.514553   0.638283   0.514048   0.252161
 718 @end example
 719
 720 @noindent
 721 Note that the eigenvectors can differ by a change of sign, since the
 722 sign of an eigenvector is arbitrary.
 723
 724 The following program illustrates the use of the nonsymmetric
 725 eigensolver, by computing the eigenvalues and eigenvectors of
 726 the Vandermonde matrix
 727 @c{$V(x;i,j) = x_i^{n - j}$}
 728 @math{V(x;i,j) = x_i^@{n - j@}}
 729 with @math{x = (-1,-2,3,4)}.
 730
 731 @example
 732 @verbatiminclude examples/eigen_nonsymm.c
 733 @end example
 734
 735 @noindent
 736 Here is the beginning of the output from the program,
 737
 738 @example
 739 $ ./a.out
 740 eigenvalue = -6.41391 + 0i
 741 eigenvector =
 742 -0.0998822 + 0i
 743 -0.111251 + 0i
 744 0.292501 + 0i
 745 0.944505 + 0i
 746 eigenvalue = 5.54555 + 3.08545i
 747 eigenvector =
 748 -0.043487 + -0.0076308i
 749 0.0642377 + -0.142127i
 750 -0.515253 + 0.0405118i
 751 -0.840592 + -0.00148565i
 752 ...
 753 @end example
 754
 755 @noindent
 756 This can be compared with the corresponding output from @sc{gnu octave},
 757
 758 @example
 759 octave> [v,d] = eig(vander([-1 -2 3 4]));
 760 octave> diag(d)
 761 ans =
 762
 763   -6.4139 + 0.0000i
 764    5.5456 + 3.0854i
 765    5.5456 - 3.0854i
 766    2.3228 + 0.0000i
 767
 768 octave> v
 769 v =
 770
 771  Columns 1 through 3:
 772
 773   -0.09988 + 0.00000i  -0.04350 - 0.00755i  -0.04350 + 0.00755i
 774   -0.11125 + 0.00000i   0.06399 - 0.14224i   0.06399 + 0.14224i
 775    0.29250 + 0.00000i  -0.51518 + 0.04142i  -0.51518 - 0.04142i
 776    0.94451 + 0.00000i  -0.84059 + 0.00000i  -0.84059 - 0.00000i
 777
 778  Column 4:
 779
 780   -0.14493 + 0.00000i
 781    0.35660 + 0.00000i
 782    0.91937 + 0.00000i
 783    0.08118 + 0.00000i
 784
 785 @end example
 786 Note that the eigenvectors corresponding to the eigenvalue
 787 @math{5.54555 + 3.08545i} are slightly different. This is because
 788 they differ by the multiplicative constant
 789 @math{0.9999984 + 0.0017674i} which has magnitude 1.
 790
 791 @node Eigenvalue and Eigenvector References
 792 @section References and Further Reading
 793
 794 Further information on the algorithms described in this section can be
 795 found in the following book,
 796
 797 @itemize @asis
 798 @item
 799 G. H. Golub, C. F. Van Loan, @cite{Matrix Computations} (3rd Ed, 1996),
 800 Johns Hopkins University Press, ISBN 0-8018-5414-8.
 801 @end itemize
 802
 803 @noindent
 804 Further information on the generalized eigensystems QZ algorithm
 805 can be found in this paper,
 806
 807 @itemize @asis
 808 @item
 809 C. Moler, G. Stewart, "An Algorithm for Generalized Matrix Eigenvalue
 810 Problems," SIAM J. Numer. Anal., Vol 10, No 2, 1973.
 811 @end itemize
 812
 813 @noindent
 814 @cindex LAPACK
 815 Eigensystem routines for very large matrices can be found in the
 816 Fortran library @sc{lapack}. The @sc{lapack} library is described in,
 817
 818 @itemize @asis
 819 @item
 820 @cite{LAPACK Users' Guide} (Third Edition, 1999), Published by SIAM,
 821 ISBN 0-89871-447-8.
 822
 823 @uref{http://www.netlib.org/lapack}
 824 @end itemize
 825
 826 @noindent
 827 The @sc{lapack} source code can be found at the website above along with
 828 an online copy of the users guide.