genv.c

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

/* eigen/genv.c
 * 
 * Copyright (C) 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>
#include <gsl/gsl_errno.h>

/*
 * This module computes the eigenvalues and eigenvectors of a
 * real generalized eigensystem A x = \lambda B x. Left and right
 * Schur vectors are optionally computed as well.
 *
 * This file contains routines based on original code from LAPACK
 * which is distributed under the modified BSD license.
 */

static int genv_get_right_eigenvectors(const gsl_matrix *S,
                                       const gsl_matrix *T,
                                       gsl_matrix *Z,
                                       gsl_matrix_complex *evec,
                                       gsl_eigen_genv_workspace *w);
static void genv_normalize_eigenvectors(gsl_vector_complex *alpha,
                                        gsl_matrix_complex *evec);

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

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_genv_workspace *
gsl_eigen_genv_alloc(const size_t n)
{
  gsl_eigen_genv_workspace *w;

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

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

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

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

  w->gen_workspace_p = gsl_eigen_gen_alloc(n);

  if (w->gen_workspace_p == 0)
    {
      gsl_eigen_genv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for gen workspace", GSL_ENOMEM);
    }

  /* compute the full Schur forms */
  gsl_eigen_gen_params(1, 1, 1, w->gen_workspace_p);

  w->work1 = gsl_vector_alloc(n);
  w->work2 = gsl_vector_alloc(n);
  w->work3 = gsl_vector_alloc(n);
  w->work4 = gsl_vector_alloc(n);
  w->work5 = gsl_vector_alloc(n);
  w->work6 = gsl_vector_alloc(n);

  if (w->work1 == 0 || w->work2 == 0 || w->work3 == 0 ||
      w->work4 == 0 || w->work5 == 0 || w->work6 == 0)
    {
      gsl_eigen_genv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for additional workspace", GSL_ENOMEM);
    }

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

/*
gsl_eigen_genv_free()
  Free workspace w
*/

void
gsl_eigen_genv_free(gsl_eigen_genv_workspace *w)
{
  if (w->gen_workspace_p)
    gsl_eigen_gen_free(w->gen_workspace_p);

  if (w->work1)
    gsl_vector_free(w->work1);

  if (w->work2)
    gsl_vector_free(w->work2);

  if (w->work3)
    gsl_vector_free(w->work3);

  if (w->work4)
    gsl_vector_free(w->work4);

  if (w->work5)
    gsl_vector_free(w->work5);

  if (w->work6)
    gsl_vector_free(w->work6);

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

/*
gsl_eigen_genv()

Solve the generalized eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda and right eigenvectors x.

Inputs: A     - general real matrix
        B     - general real matrix
        alpha - (output) where to store eigenvalue numerators
        beta  - (output) where to store eigenvalue denominators
        evec  - (output) where to store eigenvectors
        w     - workspace

Return: success or error
*/

int
gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B, gsl_vector_complex * alpha,
                gsl_vector * beta, gsl_matrix_complex *evec,
                gsl_eigen_genv_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 if (evec->size1 != N)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else
    {
      int s;
      gsl_matrix Z;

      /*
       * We need a place to store the right Schur vectors, so we will
       * treat evec as a real matrix and store them in the left
       * half - the factor of 2 in the tda corresponds to the
       * complex multiplicity
       */
      Z.size1 = N;
      Z.size2 = N;
      Z.tda = 2 * N;
      Z.data = evec->data;
      Z.block = 0;
      Z.owner = 0;

      s = gsl_eigen_gen_QZ(A, B, alpha, beta, w->Q, &Z, w->gen_workspace_p);

      if (w->Z)
        {
          /* save right Schur vectors */
          gsl_matrix_memcpy(w->Z, &Z);
        }

      /* only compute eigenvectors if we found all eigenvalues */
      if (s == GSL_SUCCESS)
        {
          /* compute eigenvectors */
          s = genv_get_right_eigenvectors(A, B, &Z, evec, w);

          if (s == GSL_SUCCESS)
            genv_normalize_eigenvectors(alpha, evec);
        }

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

/*
gsl_eigen_genv_QZ()

Solve the generalized eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda and right eigenvectors x. 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 - (output) where to store eigenvalue numerators
        beta  - (output) where to store eigenvalue denominators
        evec  - (output) where to store eigenvectors
        Q     - (output) if non-null, where to store left Schur vectors
        Z     - (output) if non-null, where to store right Schur vectors
        w     - workspace

Return: success or error
*/

int
gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
                   gsl_vector_complex * alpha, gsl_vector * beta,
                   gsl_matrix_complex * evec,
                   gsl_matrix * Q, gsl_matrix * Z,
                   gsl_eigen_genv_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_genv(A, B, alpha, beta, evec, w);

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

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

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

/*
genv_get_right_eigenvectors()
  Compute right eigenvectors of the Schur form (S, T) and then
backtransform them using the right Schur vectors to get right
eigenvectors of the original system.

Inputs: S     - upper quasi-triangular Schur form of A
        T     - upper triangular Schur form of B
        Z     - right Schur vectors
        evec  - (output) where to store eigenvectors
        w     - workspace

Return: success or error

Notes: 1) based on LAPACK routine DTGEVC
       2) eigenvectors are stored in the order that their
          eigenvalues appear in the Schur form
*/

static int
genv_get_right_eigenvectors(const gsl_matrix *S, const gsl_matrix *T,
                            gsl_matrix *Z,
                            gsl_matrix_complex *evec,
                            gsl_eigen_genv_workspace *w)
{
  const size_t N = w->size;
  const double small = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
  const double big = 1.0 / small;
  const double bignum = 1.0 / (GSL_DBL_MIN * N);
  size_t i, j, k, end;
  int is;
  double anorm, bnorm;
  double temp, temp2, temp2r, temp2i;
  double ascale, bscale;
  double salfar, sbeta;
  double acoef, bcoefr, bcoefi, acoefa, bcoefa;
  double creala, cimaga, crealb, cimagb, cre2a, cim2a, cre2b, cim2b;
  double dmin, xmax;
  double scale;
  size_t nw, na;
  int lsa, lsb;
  int complex_pair;
  gsl_complex z_zero, z_one;
  double bdiag[2] = { 0.0, 0.0 };
  double sum[4];
  int il2by2;
  size_t jr, jc, ja;
  double xscale;
  gsl_vector_complex_view ecol;
  gsl_vector_view re, im, re2, im2;

  GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);
  GSL_SET_COMPLEX(&z_one, 1.0, 0.0);

  /*
   * Compute the 1-norm of each column of (S, T) excluding elements
   * belonging to the diagonal blocks to check for possible overflow
   * in the triangular solver
   */

  anorm = fabs(gsl_matrix_get(S, 0, 0));
  if (N > 1)
    anorm += fabs(gsl_matrix_get(S, 1, 0));
  bnorm = fabs(gsl_matrix_get(T, 0, 0));

  gsl_vector_set(w->work1, 0, 0.0);
  gsl_vector_set(w->work2, 0, 0.0);

  for (j = 1; j < N; ++j)
    {
      temp = temp2 = 0.0;
      if (gsl_matrix_get(S, j, j - 1) == 0.0)
        end = j;
      else
        end = j - 1;

      for (i = 0; i < end; ++i)
        {
          temp += fabs(gsl_matrix_get(S, i, j));
          temp2 += fabs(gsl_matrix_get(T, i, j));
        }

      gsl_vector_set(w->work1, j, temp);
      gsl_vector_set(w->work2, j, temp2);

      for (i = end; i < GSL_MIN(j + 2, N); ++i)
        {
          temp += fabs(gsl_matrix_get(S, i, j));
          temp2 += fabs(gsl_matrix_get(T, i, j));
        }

      anorm = GSL_MAX(anorm, temp);
      bnorm = GSL_MAX(bnorm, temp2);
    }

  ascale = 1.0 / GSL_MAX(anorm, GSL_DBL_MIN);
  bscale = 1.0 / GSL_MAX(bnorm, GSL_DBL_MIN);

  complex_pair = 0;
  for (k = 0; k < N; ++k)
    {
      size_t je = N - 1 - k;

      if (complex_pair)
        {
          complex_pair = 0;
          continue;
        }

      nw = 1;
      if (je > 0)
        {
          if (gsl_matrix_get(S, je, je - 1) != 0.0)
            {
              complex_pair = 1;
              nw = 2;
            }
        }

      if (!complex_pair)
        {
          if (fabs(gsl_matrix_get(S, je, je)) <= GSL_DBL_MIN &&
              fabs(gsl_matrix_get(T, je, je)) <= GSL_DBL_MIN)
            {
              /* singular matrix pencil - unit eigenvector */
              for (i = 0; i < N; ++i)
                gsl_matrix_complex_set(evec, i, je, z_zero);

              gsl_matrix_complex_set(evec, je, je, z_one);

              continue;
            }

          /* clear vector */
          for (i = 0; i < N; ++i)
            gsl_vector_set(w->work3, i, 0.0);
        }
      else
        {
          /* clear vectors */
          for (i = 0; i < N; ++i)
            {
              gsl_vector_set(w->work3, i, 0.0);
              gsl_vector_set(w->work4, i, 0.0);
            }
        }

      if (!complex_pair)
        {
          /* real eigenvalue */

          temp = 1.0 / GSL_MAX(GSL_DBL_MIN,
                               GSL_MAX(fabs(gsl_matrix_get(S, je, je)) * ascale,
                                       fabs(gsl_matrix_get(T, je, je)) * bscale));
          salfar = (temp * gsl_matrix_get(S, je, je)) * ascale;
          sbeta = (temp * gsl_matrix_get(T, je, je)) * bscale;
          acoef = sbeta * ascale;
          bcoefr = salfar * bscale;
          bcoefi = 0.0;

          /* scale to avoid underflow */
          scale = 1.0;
          lsa = fabs(sbeta) >= GSL_DBL_MIN && fabs(acoef) < small;
          lsb = fabs(salfar) >= GSL_DBL_MIN && fabs(bcoefr) < small;
          if (lsa)
            scale = (small / fabs(sbeta)) * GSL_MIN(anorm, big);
          if (lsb)
            scale = GSL_MAX(scale, (small / fabs(salfar)) * GSL_MIN(bnorm, big));

          if (lsa || lsb)
            {