gen.c

来自「math library from gnu」· C语言 代码 · 共 2,115 行 · 第 1/4 页

/* eigen/gen.c
 * 
 * Copyright (C) 2006, 2007 Patrick Alken
 * 
 * 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <stdlib.h>
#include <math.h>

#include <config.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_matrix.h>

/*
 * This module computes the eigenvalues of a real generalized
 * eigensystem A x = \lambda B x. Left and right Schur vectors
 * are optionally computed as well.
 * 
 * Based on the algorithm from Moler and Stewart
 * [1] C. Moler, G. Stewart, "An Algorithm for Generalized Matrix
 *     Eigenvalue Problems", SIAM J. Numer. Anal., Vol 10, No 2, 1973.
 *
 * This algorithm is also described in the book
 * [2] Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm 7.7.3
 *
 * This file contains routines based on original code from LAPACK
 * which is distributed under the modified BSD license.
 */

#define GEN_ESHIFT_COEFF     (1.736)

static void gen_schur_decomp(gsl_matrix *H, gsl_matrix *R,
                             gsl_vector_complex *alpha, gsl_vector *beta,
                             gsl_eigen_gen_workspace *w);
static inline int gen_qzstep(gsl_matrix *H, gsl_matrix *R,
                             gsl_eigen_gen_workspace *w);
static inline void gen_qzstep_d(gsl_matrix *H, gsl_matrix *R,
                                gsl_eigen_gen_workspace *w);
static void gen_tri_split_top(gsl_matrix *H, gsl_matrix *R,
                              gsl_eigen_gen_workspace *w);
static inline void gen_tri_chase_zero(gsl_matrix *H, gsl_matrix *R,
                                      size_t q,
                                      gsl_eigen_gen_workspace *w);
static inline void gen_tri_zero_H(gsl_matrix *H, gsl_matrix *R,
                                  gsl_eigen_gen_workspace *w);
static inline size_t gen_search_small_elements(gsl_matrix *H,
                                               gsl_matrix *R,
                                               int *flag,
                                               gsl_eigen_gen_workspace *w);
static int gen_schur_standardize1(gsl_matrix *A, gsl_matrix *B,
                                  double *alphar, double *beta,
                                  gsl_eigen_gen_workspace *w);
static int gen_schur_standardize2(gsl_matrix *A, gsl_matrix *B,
                                  gsl_complex *alpha1,
                                  gsl_complex *alpha2,
                                  double *beta1, double *beta2,
                                  gsl_eigen_gen_workspace *w);
static int gen_compute_eigenvals(gsl_matrix *A, gsl_matrix *B,
                                 gsl_complex *alpha1,
                                 gsl_complex *alpha2, double *beta1,
                                 double *beta2);
static void gen_store_eigval1(const gsl_matrix *H, const double a,
                              const double b, gsl_vector_complex *alpha,
                              gsl_vector *beta,
                              gsl_eigen_gen_workspace *w);
static void gen_store_eigval2(const gsl_matrix *H,
                              const gsl_complex *alpha1,
                              const double beta1,
                              const gsl_complex *alpha2,
                              const double beta2,
                              gsl_vector_complex *alpha,
                              gsl_vector *beta,
                              gsl_eigen_gen_workspace *w);
static inline size_t gen_get_submatrix(const gsl_matrix *A,
                                       const gsl_matrix *B);

/*FIX**/
inline static double normF (gsl_matrix * A);
inline static void create_givens (const double a, const double b, double *c, double *s);

/*
gsl_eigen_gen_alloc()

Allocate a workspace for solving the generalized eigenvalue problem.
The size of this workspace is O(n)

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_gen_workspace *
gsl_eigen_gen_alloc(const size_t n)
{
  gsl_eigen_gen_workspace *w;

  if (n == 0)
    {
      GSL_ERROR_NULL ("matrix dimension must be positive integer",
                      GSL_EINVAL);
    }

  w = (gsl_eigen_gen_workspace *) calloc (1, sizeof (gsl_eigen_gen_workspace));

  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->size = n;
  w->max_iterations = 30 * n;
  w->n_evals = 0;
  w->n_iter = 0;
  w->needtop = 0;
  w->atol = 0.0;
  w->btol = 0.0;
  w->ascale = 0.0;
  w->bscale = 0.0;
  w->eshift = 0.0;
  w->H = NULL;
  w->R = NULL;
  w->compute_s = 0;
  w->compute_t = 0;
  w->Q = NULL;
  w->Z = NULL;

  w->work = gsl_vector_alloc(n);

  if (w->work == 0)
    {
      gsl_eigen_gen_free(w);
      GSL_ERROR_NULL ("failed to allocate space for additional workspace", GSL_ENOMEM);
    }

  return (w);
} /* gsl_eigen_gen_alloc() */

/*
gsl_eigen_gen_free()
  Free workspace w
*/

void
gsl_eigen_gen_free (gsl_eigen_gen_workspace * w)
{
  if (w->work)
    gsl_vector_free(w->work);

  free(w);
} /* gsl_eigen_gen_free() */

/*
gsl_eigen_gen_params()
  Set parameters which define how we solve the eigenvalue problem

Inputs: compute_s - 1 if we want to compute S, 0 if not
        compute_t - 1 if we want to compute T, 0 if not
        balance   - 1 if we want to balance matrices, 0 if not
        w         - gen workspace

Return: none
*/

void
gsl_eigen_gen_params (const int compute_s, const int compute_t,
                      const int balance, gsl_eigen_gen_workspace *w)
{
  w->compute_s = compute_s;
  w->compute_t = compute_t;
} /* gsl_eigen_gen_params() */

/*
gsl_eigen_gen()

Solve the generalized eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda.

Inputs: A     - general real matrix
        B     - general real matrix
        alpha - where to store eigenvalue numerators
        beta  - where to store eigenvalue denominators
        w     - workspace

Return: success or error
*/

int
gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B, gsl_vector_complex * alpha,
               gsl_vector * beta, gsl_eigen_gen_workspace * w)
{
  const size_t N = A->size1;

  /* check matrix and vector sizes */

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if ((N != B->size1) || (N != B->size2))
    {
      GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
    }
  else if (alpha->size != N || beta->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (w->size != N)
    {
      GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
    }
  else
    {
      double anorm, bnorm;

      /* compute the Hessenberg-Triangular reduction of (A, B) */
      gsl_linalg_hesstri_decomp(A, B, w->Q, w->Z, w->work);

      /* save pointers to original matrices */
      w->H = A;
      w->R = B;

      w->n_evals = 0;
      w->n_iter = 0;
      w->eshift = 0.0;

      /* determine if we need to compute top indices in QZ step */
      w->needtop = w->Q != 0 || w->Z != 0 || w->compute_t || w->compute_s;

      /* compute matrix norms */
      anorm = normF(A);
      bnorm = normF(B);

      /* compute tolerances and scaling factors */
      w->atol = GSL_MAX(GSL_DBL_MIN, GSL_DBL_EPSILON * anorm);
      w->btol = GSL_MAX(GSL_DBL_MIN, GSL_DBL_EPSILON * bnorm);
      w->ascale = 1.0 / GSL_MAX(GSL_DBL_MIN, anorm);
      w->bscale = 1.0 / GSL_MAX(GSL_DBL_MIN, bnorm);

      /* compute the generalized Schur decomposition and eigenvalues */
      gen_schur_decomp(A, B, alpha, beta, w);

      if (w->n_evals != N)
        return GSL_EMAXITER;

      return GSL_SUCCESS;
    }
} /* gsl_eigen_gen() */

/*
gsl_eigen_gen_QZ()

Solve the generalized eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda. Optionally compute left and/or right
Schur vectors Q and Z which satisfy:

A = Q S Z^t
B = Q T Z^t

where (S, T) is the generalized Schur form of (A, B)

Inputs: A     - general real matrix
        B     - general real matrix
        alpha - where to store eigenvalue numerators
        beta  - where to store eigenvalue denominators
        Q     - if non-null, where to store left Schur vectors
        Z     - if non-null, where to store right Schur vectors
        w     - workspace

Return: success or error
*/

int
gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
                  gsl_vector_complex * alpha, gsl_vector * beta,
                  gsl_matrix * Q, gsl_matrix * Z,
                  gsl_eigen_gen_workspace * w)
{
  if (Q && (A->size1 != Q->size1 || A->size1 != Q->size2))
    {
      GSL_ERROR("Q matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (Z && (A->size1 != Z->size1 || A->size1 != Z->size2))
    {
      GSL_ERROR("Z matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      int s;

      w->Q = Q;
      w->Z = Z;

      s = gsl_eigen_gen(A, B, alpha, beta, w);

      w->Q = NULL;
      w->Z = NULL;

      return s;
    }
} /* gsl_eigen_gen_QZ() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

/*
gen_schur_decomp()
  Compute the generalized Schur decomposition of the matrix pencil
(H, R) which is in Hessenberg-Triangular form

Inputs: H     - upper hessenberg matrix
        R     - upper triangular matrix
        alpha - (output) where to store eigenvalue numerators
        beta  - (output) where to store eigenvalue denominators
        w     - workspace

Return: none

Notes: 1) w->n_evals is updated to keep track of how many eigenvalues
          are found
*/

static void
gen_schur_decomp(gsl_matrix *H, gsl_matrix *R, gsl_vector_complex *alpha,
                 gsl_vector *beta, gsl_eigen_gen_workspace *w)
{
  size_t N;
  gsl_matrix_view h, r;
  gsl_matrix_view vh, vr;
  size_t q;             /* index of small subdiagonal element */
  gsl_complex z1, z2;   /* complex values */
  double a, b;
  int s;
  int flag;

  N = H->size1;

  h = gsl_matrix_submatrix(H, 0, 0, N, N);
  r = gsl_matrix_submatrix(R, 0, 0, N, N);

  while ((N > 1) && (w->n_iter)++ < w->max_iterations)
    {
      q = gen_search_small_elements(&h.matrix, &r.matrix, &flag, w);

      if (flag == 0)
        {
          /* no small elements found - do a QZ sweep */
          s = gen_qzstep(&h.matrix, &r.matrix, w);

          if (s == GSL_CONTINUE)
            {
              /*
               * (h, r) is a 2-by-2 block with complex eigenvalues -
               * standardize and read off eigenvalues
               */
              s = gen_schur_standardize2(&h.matrix,
                                         &r.matrix,
                                         &z1,
                                         &z2,
                                         &a,
                                         &b,
                                         w);
              if (s != GSL_SUCCESS)
                {
                  /*
                   * if we get here, then the standardization process
                   * perturbed the eigenvalues onto the real line -
                   * continue QZ iteration to break them into 1-by-1
                   * blocks
                   */
                  continue;
                }

              gen_store_eigval2(&h.matrix, &z1, a, &z2, b, alpha, beta, w);
              N = 0;
            } /* if (s) */

          continue;
        } /* if (flag == 0) */
      else if (flag == 2)
        {
          if (q == 0)
            {
              /*
               * the leading element of R is zero, split off a block
               * at the top
               */
              gen_tri_split_top(&h.matrix, &r.matrix, w);
            }
          else
            {
              /*
               * we found a small element on the diagonal of R - chase the
               * zero to the bottom of the active block and then zero
               * H(n, n - 1) to split off a 1-by-1 block
               */

              if (q != N - 1)
                gen_tri_chase_zero(&h.matrix, &r.matrix, q, w);

              /* now zero H(n, n - 1) */
              gen_tri_zero_H(&h.matrix, &r.matrix, w);
            }

          /* continue so the next iteration detects the zero in H */
          continue;
        }

      /*
       * a small subdiagonal element of H was found - one or two
       * eigenvalues have converged or the matrix has split into
       * two smaller matrices
       */

      if (q == (N - 1))
        {
          /*
           * the last subdiagonal element of the hessenberg matrix is 0 -
           * H_{NN} / R_{NN} is a real eigenvalue - standardize so
           * R_{NN} > 0
           */

          vh = gsl_matrix_submatrix(&h.matrix, q, q, 1, 1);
          vr = gsl_matrix_submatrix(&r.matrix, q, q, 1, 1);
          gen_schur_standardize1(&vh.matrix, &vr.matrix, &a, &b, w);

          gen_store_eigval1(&vh.matrix, a, b, alpha, beta, w);

          --N;
          h = gsl_matrix_submatrix(&h.matrix, 0, 0, N, N);
          r = gsl_matrix_submatrix(&r.matrix, 0, 0, N, N);
        }
      else if (q == (N - 2))
        {
          /* bottom right 2-by-2 block may have converged */

          vh = gsl_matrix_submatrix(&h.matrix, q, q, 2, 2);
          vr = gsl_matrix_submatrix(&r.matrix, q, q, 2, 2);
          s = gen_schur_standardize2(&vh.matrix,
                                     &vr.matrix,
                                     &z1,
                                     &z2,
                                     &a,
                                     &b,
                                     w);
          if (s != GSL_SUCCESS)
            {
              /*
               * this 2-by-2 block contains real eigenvalues that
               * have not yet separated into 1-by-1 blocks -
               * recursively call gen_schur_decomp() to finish off
               * this block
               */
              gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);
            }
          else
            {
              /* we got 2 complex eigenvalues */

              gen_store_eigval2(&vh.matrix, &z1, a, &z2, b, alpha, beta, w);
            }

          N -= 2;
          h = gsl_matrix_submatrix(&h.matrix, 0, 0, N, N);
          r = gsl_matrix_submatrix(&r.matrix, 0, 0, N, N);
        }
      else if (q == 1)
        {
          /* H_{11} / R_{11} is an eigenvalue */

          vh = gsl_matrix_submatrix(&h.matrix, 0, 0, 1, 1);
          vr = gsl_matrix_submatrix(&r.matrix, 0, 0, 1, 1);
          gen_schur_standardize1(&vh.matrix, &vr.matrix, &a, &b, w);

          gen_store_eigval1(&vh.matrix, a, b, alpha, beta, w);

          --N;
          h = gsl_matrix_submatrix(&h.matrix, 1, 1, N, N);
          r = gsl_matrix_submatrix(&r.matrix, 1, 1, N, N);
        }
      else if (q == 2)
        {
          /* upper left 2-by-2 block may have converged */

          vh = gsl_matrix_submatrix(&h.matrix, 0, 0, 2, 2);
          vr = gsl_matrix_submatrix(&r.matrix, 0, 0, 2, 2);
          s = gen_schur_standardize2(&vh.matrix,
                                     &vr.matrix,
                                     &z1,
                                     &z2,
                                     &a,
                                     &b,
                                     w);
          if (s != GSL_SUCCESS)
            {
              /*
               * this 2-by-2 block contains real eigenvalues that
               * have not yet separated into 1-by-1 blocks -
               * recursively call gen_schur_decomp() to finish off
               * this block
               */
              gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);
            }
          else
            {