nonsymmv.c

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

/* eigen/nonsymmv.c
 * 
 * Copyright (C) 2006 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 <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_matrix.h>

/*
 * This module computes the eigenvalues and eigenvectors of a real
 * nonsymmetric matrix.
 * 
 * This file contains routines based on original code from LAPACK
 * which is distributed under the modified BSD license. The LAPACK
 * routines used are DTREVC and DLALN2.
 */

#define GSL_NONSYMMV_SMLNUM (2.0 * GSL_DBL_MIN)
#define GSL_NONSYMMV_BIGNUM ((1.0 - GSL_DBL_EPSILON) / GSL_NONSYMMV_SMLNUM)

static void nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
                                            gsl_vector_complex *eval,
                                            gsl_matrix_complex *evec,
                                            gsl_eigen_nonsymmv_workspace *w);
static void nonsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
                                            gsl_matrix_complex *evec);

/*
gsl_eigen_nonsymmv_alloc()

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

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_nonsymmv_workspace *
gsl_eigen_nonsymmv_alloc(const size_t n)
{
  gsl_eigen_nonsymmv_workspace *w;

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

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

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

  w->size = n;
  w->Z = NULL;
  w->nonsymm_workspace_p = gsl_eigen_nonsymm_alloc(n);

  if (w->nonsymm_workspace_p == 0)
    {
      gsl_eigen_nonsymmv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for nonsymm workspace", GSL_ENOMEM);
    }

  /*
   * set parameters to compute the full Schur form T and balance
   * the matrices
   */
  gsl_eigen_nonsymm_params(1, 1, w->nonsymm_workspace_p);

  w->work = gsl_vector_alloc(n);
  w->work2 = gsl_vector_alloc(n);
  w->work3 = gsl_vector_alloc(n);
  if (w->work == 0 || w->work2 == 0 || w->work3 == 0)
    {
      gsl_eigen_nonsymmv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for nonsymmv additional workspace", GSL_ENOMEM);
    }

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

/*
gsl_eigen_nonsymmv_free()
  Free workspace w
*/

void
gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w)
{
  if (w->nonsymm_workspace_p)
    gsl_eigen_nonsymm_free(w->nonsymm_workspace_p);

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

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

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

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

/*
gsl_eigen_nonsymmv()

Solve the nonsymmetric eigensystem problem

A x = \lambda x

for the eigenvalues \lambda and right eigenvectors x

Inputs: A    - general real matrix
        eval - where to store eigenvalues
        evec - where to store eigenvectors
        w    - workspace

Return: success or error
*/

int
gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
                    gsl_matrix_complex * evec,
                    gsl_eigen_nonsymmv_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 (eval->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  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 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;

      /* compute eigenvalues, Schur form, and Schur vectors */
      s = gsl_eigen_nonsymm_Z(A, eval, &Z, w->nonsymm_workspace_p);

      if (w->Z)
        {
          /*
           * save the Schur vectors in user supplied matrix, since
           * they will be destroyed when computing eigenvectors
           */
          gsl_matrix_memcpy(w->Z, &Z);
        }

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

          /* normalize so that Euclidean norm is 1 */
          nonsymmv_normalize_eigenvectors(eval, evec);
        }

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

/*
gsl_eigen_nonsymmv_Z()
  Compute eigenvalues and eigenvectors of a real nonsymmetric matrix
and also save the Schur vectors. See comments in gsl_eigen_nonsymm_Z
for more information.

Inputs: A    - real nonsymmetric matrix
        eval - where to store eigenvalues
        evec - where to store eigenvectors
        Z    - where to store Schur vectors
        w    - nonsymmv workspace

Return: success or error
*/

int
gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
                      gsl_matrix_complex * evec, gsl_matrix * Z,
                      gsl_eigen_nonsymmv_workspace * w)
{
  /* check matrix and vector sizes */

  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues/eigenvectors", GSL_ENOTSQR);
    }
  else if (eval->size != A->size1)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (evec->size1 != A->size1)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
    {
      GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      int s;

      w->Z = Z;

      s = gsl_eigen_nonsymmv(A, eval, evec, w);

      w->Z = NULL;

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

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

/*
nonsymmv_get_right_eigenvectors()
  Compute the right eigenvectors of the Schur form T and then
backtransform them using the Schur vectors to get right eigenvectors of
the original matrix.

Inputs: T    - Schur form
        Z    - Schur vectors
        eval - where to store eigenvalues (to ensure that the
               correct eigenvalue is stored in the same position
               as the eigenvectors)
        evec - where to store eigenvectors
        w    - nonsymmv workspace

Return: none

Notes: 1) based on LAPACK routine DTREVC - the algorithm used is
          backsubstitution on the upper quasi triangular system T
          followed by backtransformation by Z to get vectors of the
          original matrix.

       2) The Schur vectors in Z are destroyed and replaced with
          eigenvectors stored with the same storage scheme as DTREVC.
          The eigenvectors are also stored in 'evec'

       3) The matrix T is unchanged on output

       4) Each eigenvector is normalized so that the element of
          largest magnitude has magnitude 1; here the magnitude of
          a complex number (x,y) is taken to be |x| + |y|
*/

static void
nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
                                gsl_vector_complex *eval,
                                gsl_matrix_complex *evec,
                                gsl_eigen_nonsymmv_workspace *w)
{
  const size_t N = T->size1;
  const double smlnum = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
  const double bignum = (1.0 - GSL_DBL_EPSILON) / smlnum;
  int i;              /* looping */
  size_t iu,          /* looping */
         ju,