|
Packit |
67cb25 |
/* eigen/nonsymm.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2006 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_eigen.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector_complex.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* This module computes the eigenvalues of a real nonsymmetric
|
|
Packit |
67cb25 |
* matrix, using the double shift Francis method.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* See the references in francis.c.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This module gets the matrix ready by balancing it and
|
|
Packit |
67cb25 |
* reducing it to Hessenberg form before passing it to the
|
|
Packit |
67cb25 |
* francis module.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_alloc()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Allocate a workspace for solving the nonsymmetric eigenvalue problem.
|
|
Packit |
67cb25 |
The size of this workspace is O(2n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: n - size of matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: pointer to workspace
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_workspace *
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_alloc(const size_t n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_workspace *w;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("matrix dimension must be positive integer",
|
|
Packit |
67cb25 |
GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w = (gsl_eigen_nonsymm_workspace *)
|
|
Packit |
67cb25 |
calloc (1, sizeof (gsl_eigen_nonsymm_workspace));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->size = n;
|
|
Packit |
67cb25 |
w->Z = NULL;
|
|
Packit |
67cb25 |
w->do_balance = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->diag = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->diag == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_free(w);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for balancing vector", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->tau = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->tau == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_free(w);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for hessenberg coefficients", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->francis_workspace_p = gsl_eigen_francis_alloc();
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->francis_workspace_p == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_free(w);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for francis workspace", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (w);
|
|
Packit |
67cb25 |
} /* gsl_eigen_nonsymm_alloc() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_free()
|
|
Packit |
67cb25 |
Free workspace w
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
RETURN_IF_NULL (w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->tau)
|
|
Packit |
67cb25 |
gsl_vector_free(w->tau);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->diag)
|
|
Packit |
67cb25 |
gsl_vector_free(w->diag);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->francis_workspace_p)
|
|
Packit |
67cb25 |
gsl_eigen_francis_free(w->francis_workspace_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
free(w);
|
|
Packit |
67cb25 |
} /* gsl_eigen_nonsymm_free() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_params()
|
|
Packit |
67cb25 |
Set some parameters which define how we solve the eigenvalue
|
|
Packit |
67cb25 |
problem.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: compute_t - 1 if we want to compute T, 0 if not
|
|
Packit |
67cb25 |
balance - 1 if we want to balance the matrix, 0 if not
|
|
Packit |
67cb25 |
w - nonsymm workspace
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_params (const int compute_t, const int balance,
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_francis_T(compute_t, w->francis_workspace_p);
|
|
Packit |
67cb25 |
w->do_balance = balance;
|
|
Packit |
67cb25 |
} /* gsl_eigen_nonsymm_params() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solve the nonsymmetric eigenvalue problem
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A x = \lambda x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for the eigenvalues \lambda using the Francis method.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here we compute the real Schur form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
T = Z^t A Z
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with the diagonal blocks of T giving us the eigenvalues.
|
|
Packit |
67cb25 |
Z is a matrix of Schur vectors which is not computed by
|
|
Packit |
67cb25 |
this algorithm. See gsl_eigen_nonsymm_Z().
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - general real matrix
|
|
Packit |
67cb25 |
eval - where to store eigenvalues
|
|
Packit |
67cb25 |
w - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: If T is computed, it is stored in A on output. Otherwise
|
|
Packit |
67cb25 |
the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
|
|
Packit |
67cb25 |
blocks.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check matrix and vector sizes */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (eval->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->do_balance)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* balance the matrix */
|
|
Packit |
67cb25 |
gsl_linalg_balance_matrix(A, w->diag);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute the Hessenberg reduction of A */
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_decomp(A, w->tau);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->Z)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* initialize the matrix Z to U, which is the matrix used
|
|
Packit |
67cb25 |
* to construct the Hessenberg reduction.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute U and store it in Z */
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_unpack(A, w->tau, w->Z);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* find the eigenvalues and Schur vectors */
|
|
Packit |
67cb25 |
s = gsl_eigen_francis_Z(A, eval, w->Z, w->francis_workspace_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->do_balance)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* The Schur vectors in Z are the vectors for the balanced
|
|
Packit |
67cb25 |
* matrix. We now must undo the balancing to get the
|
|
Packit |
67cb25 |
* vectors for the original matrix A.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_linalg_balance_accum(w->Z, w->diag);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* find the eigenvalues only */
|
|
Packit |
67cb25 |
s = gsl_eigen_francis(A, eval, w->francis_workspace_p);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->n_evals = w->francis_workspace_p->n_evals;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_eigen_nonsymm() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_Z()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solve the nonsymmetric eigenvalue problem
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A x = \lambda x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for the eigenvalues \lambda.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here we compute the real Schur form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
T = Z^t A Z
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with the diagonal blocks of T giving us the eigenvalues.
|
|
Packit |
67cb25 |
Z is the matrix of Schur vectors.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - general real matrix
|
|
Packit |
67cb25 |
eval - where to store eigenvalues
|
|
Packit |
67cb25 |
Z - where to store Schur vectors
|
|
Packit |
67cb25 |
w - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: If T is computed, it is stored in A on output. Otherwise
|
|
Packit |
67cb25 |
the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
|
|
Packit |
67cb25 |
blocks.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* check matrix and vector sizes */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (A->size1 != A->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (eval->size != A->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->Z = Z;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s = gsl_eigen_nonsymm(A, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->Z = NULL;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_eigen_nonsymm_Z() */
|