|
Packit |
67cb25 |
/* eigen/francis.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2006, 2007 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 |
#include <gsl/gsl_complex.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_complex_math.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* This module computes the eigenvalues of a real upper hessenberg
|
|
Packit |
67cb25 |
* matrix, using the classical double shift Francis QR algorithm.
|
|
Packit |
67cb25 |
* It will also optionally compute the full Schur form and matrix of
|
|
Packit |
67cb25 |
* Schur vectors.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* See Golub & Van Loan, "Matrix Computations" (3rd ed),
|
|
Packit |
67cb25 |
* algorithm 7.5.2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* exceptional shift coefficients - these values are from LAPACK DLAHQR */
|
|
Packit |
67cb25 |
#define GSL_FRANCIS_COEFF1 (0.75)
|
|
Packit |
67cb25 |
#define GSL_FRANCIS_COEFF2 (-0.4375)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline void francis_schur_decomp(gsl_matrix * H,
|
|
Packit |
67cb25 |
gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace * w);
|
|
Packit |
67cb25 |
static inline size_t francis_search_subdiag_small_elements(gsl_matrix * A);
|
|
Packit |
67cb25 |
static inline int francis_qrstep(gsl_matrix * H,
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace * w);
|
|
Packit |
67cb25 |
static inline void francis_schur_standardize(gsl_matrix *A,
|
|
Packit |
67cb25 |
gsl_complex *eval1,
|
|
Packit |
67cb25 |
gsl_complex *eval2,
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace *w);
|
|
Packit |
67cb25 |
static inline size_t francis_get_submatrix(gsl_matrix *A, gsl_matrix *B);
|
|
Packit |
67cb25 |
static void francis_standard_form(gsl_matrix *A, double *cs, double *sn);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_francis_alloc()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Allocate a workspace for solving the nonsymmetric eigenvalue problem.
|
|
Packit |
67cb25 |
The size of this workspace is O(1)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: none
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: pointer to workspace
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace *
|
|
Packit |
67cb25 |
gsl_eigen_francis_alloc(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace *w;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w = (gsl_eigen_francis_workspace *)
|
|
Packit |
67cb25 |
calloc (1, sizeof (gsl_eigen_francis_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 |
/* these are filled in later */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->size = 0;
|
|
Packit |
67cb25 |
w->max_iterations = 0;
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
w->n_evals = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->compute_t = 0;
|
|
Packit |
67cb25 |
w->Z = NULL;
|
|
Packit |
67cb25 |
w->H = NULL;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (w);
|
|
Packit |
67cb25 |
} /* gsl_eigen_francis_alloc() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_francis_free()
|
|
Packit |
67cb25 |
Free francis workspace w
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_eigen_francis_free (gsl_eigen_francis_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
RETURN_IF_NULL (w);
|
|
Packit |
67cb25 |
free(w);
|
|
Packit |
67cb25 |
} /* gsl_eigen_francis_free() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_francis_T()
|
|
Packit |
67cb25 |
Called when we want to compute the Schur form T, or no longer
|
|
Packit |
67cb25 |
compute the Schur form T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: compute_t - 1 to compute T, 0 to not compute T
|
|
Packit |
67cb25 |
w - francis workspace
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_eigen_francis_T (const int compute_t, gsl_eigen_francis_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
w->compute_t = compute_t;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_francis()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solve the nonsymmetric eigenvalue problem
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H x = \lambda x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for the eigenvalues \lambda using algorithm 7.5.2 of
|
|
Packit |
67cb25 |
Golub & Van Loan, "Matrix Computations" (3rd ed)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: H - upper hessenberg matrix
|
|
Packit |
67cb25 |
eval - where to store eigenvalues
|
|
Packit |
67cb25 |
w - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error - if error code is returned,
|
|
Packit |
67cb25 |
then the QR procedure did not converge in the
|
|
Packit |
67cb25 |
allowed number of iterations. In the event of non-
|
|
Packit |
67cb25 |
convergence, the number of eigenvalues found will
|
|
Packit |
67cb25 |
still be stored in the beginning of eval,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: On output, the diagonal of H contains 1-by-1 or 2-by-2
|
|
Packit |
67cb25 |
blocks containing the eigenvalues. If T is desired,
|
|
Packit |
67cb25 |
H will contain the full Schur form on output.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* check matrix and vector sizes */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (H->size1 != H->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (eval->size != H->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = H->size1;
|
|
Packit |
67cb25 |
int j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Set internal parameters which depend on matrix size.
|
|
Packit |
67cb25 |
* The Francis solver can be called with any size matrix
|
|
Packit |
67cb25 |
* since the workspace does not depend on N.
|
|
Packit |
67cb25 |
* Furthermore, multishift solvers which call the Francis
|
|
Packit |
67cb25 |
* solver may need to call it with different sized matrices
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
w->size = N;
|
|
Packit |
67cb25 |
w->max_iterations = 30 * N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* save a pointer to original matrix since francis_schur_decomp
|
|
Packit |
67cb25 |
* is recursive
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
w->H = H;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
w->n_evals = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* zero out the first two subdiagonals (below the main subdiagonal)
|
|
Packit |
67cb25 |
* needed as scratch space by the QR sweep routine
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
for (j = 0; j < (int) N - 3; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(H, (size_t) j + 2, (size_t) j, 0.0);
|
|
Packit |
67cb25 |
gsl_matrix_set(H, (size_t) j + 3, (size_t) j, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N > 2)
|
|
Packit |
67cb25 |
gsl_matrix_set(H, N - 1, N - 3, 0.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* compute Schur decomposition of H and store eigenvalues
|
|
Packit |
67cb25 |
* into eval
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
francis_schur_decomp(H, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->n_evals != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("maximum iterations reached without finding all eigenvalues", GSL_EMAXITER);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_eigen_francis() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_eigen_francis_Z()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solve the nonsymmetric eigenvalue problem for a Hessenberg
|
|
Packit |
67cb25 |
matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H x = \lambda x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for the eigenvalues \lambda using the Francis double-shift
|
|
Packit |
67cb25 |
method.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here we compute the real Schur form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
T = Q^t H Q
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with the diagonal blocks of T giving us the eigenvalues.
|
|
Packit |
67cb25 |
Q is the matrix of Schur vectors.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Originally, H was obtained from a general nonsymmetric matrix
|
|
Packit |
67cb25 |
A via a transformation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H = U^t A U
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
so that
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
T = (UQ)^t A (UQ) = Z^t A Z
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Z is the matrix of Schur vectors computed by this algorithm
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: H - upper hessenberg matrix
|
|
Packit |
67cb25 |
eval - where to store eigenvalues
|
|
Packit |
67cb25 |
Z - where to store Schur vectors
|
|
Packit |
67cb25 |
w - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: 1) If T is computed, it is stored in H on output. Otherwise,
|
|
Packit |
67cb25 |
the diagonal of H will contain 1-by-1 and 2-by-2 blocks
|
|
Packit |
67cb25 |
containing the eigenvalues.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
2) The matrix Z must be initialized to the Hessenberg
|
|
Packit |
67cb25 |
similarity matrix U. Or if you want the eigenvalues
|
|
Packit |
67cb25 |
of H, initialize Z to the identity matrix.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
gsl_matrix * Z, gsl_eigen_francis_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* set internal Z pointer so we know to accumulate transformations */
|
|
Packit |
67cb25 |
w->Z = Z;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s = gsl_eigen_francis(H, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->Z = NULL;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
} /* gsl_eigen_francis_Z() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/********************************************
|
|
Packit |
67cb25 |
* INTERNAL ROUTINES *
|
|
Packit |
67cb25 |
********************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
francis_schur_decomp()
|
|
Packit |
67cb25 |
Compute the Schur decomposition of the matrix H
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: H - hessenberg matrix
|
|
Packit |
67cb25 |
eval - where to store eigenvalues
|
|
Packit |
67cb25 |
w - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: none
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline void
|
|
Packit |
67cb25 |
francis_schur_decomp(gsl_matrix * H, gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view m; /* active matrix we are working on */
|
|
Packit |
67cb25 |
size_t N; /* size of matrix */
|
|
Packit |
67cb25 |
size_t q; /* index of small subdiagonal element */
|
|
Packit |
67cb25 |
gsl_complex lambda1, /* eigenvalues */
|
|
Packit |
67cb25 |
lambda2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
N = H->size1;
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(H, 0, 0, N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while ((N > 2) && ((w->n_iter)++ < w->max_iterations))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
q = francis_search_subdiag_small_elements(&m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (q == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* no small subdiagonal element found - perform a QR
|
|
Packit |
67cb25 |
* sweep on the active reduced hessenberg matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
francis_qrstep(&m.matrix, w);
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* a small subdiagonal element was found - one or two eigenvalues
|
|
Packit |
67cb25 |
* have converged or the matrix has split into two smaller matrices
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (q == (N - 1))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* the last subdiagonal element of the matrix is 0 -
|
|
Packit |
67cb25 |
* m_{NN} is a real eigenvalue
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&lambda1,
|
|
Packit |
67cb25 |
gsl_matrix_get(&m.matrix, q, q), 0.0);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals, lambda1);
|
|
Packit |
67cb25 |
w->n_evals += 1;
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
--N;
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(&m.matrix, 0, 0, N, N);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (q == (N - 2))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* The bottom right 2-by-2 block of m is an eigenvalue
|
|
Packit |
67cb25 |
* system
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v = gsl_matrix_submatrix(&m.matrix, q, q, 2, 2);
|
|
Packit |
67cb25 |
francis_schur_standardize(&v.matrix, &lambda1, &lambda2, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals, lambda1);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals + 1, lambda2);
|
|
Packit |
67cb25 |
w->n_evals += 2;
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
N -= 2;
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(&m.matrix, 0, 0, N, N);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (q == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* the first matrix element is an eigenvalue */
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&lambda1,
|
|
Packit |
67cb25 |
gsl_matrix_get(&m.matrix, 0, 0), 0.0);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals, lambda1);
|
|
Packit |
67cb25 |
w->n_evals += 1;
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
--N;
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(&m.matrix, 1, 1, N, N);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (q == 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* the upper left 2-by-2 block is an eigenvalue system */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v = gsl_matrix_submatrix(&m.matrix, 0, 0, 2, 2);
|
|
Packit |
67cb25 |
francis_schur_standardize(&v.matrix, &lambda1, &lambda2, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals, lambda1);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals + 1, lambda2);
|
|
Packit |
67cb25 |
w->n_evals += 2;
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
N -= 2;
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(&m.matrix, 2, 2, N, N);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* There is a zero element on the subdiagonal somewhere
|
|
Packit |
67cb25 |
* in the middle of the matrix - we can now operate
|
|
Packit |
67cb25 |
* separately on the two submatrices split by this
|
|
Packit |
67cb25 |
* element. q is the row index of the zero element.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* operate on lower right (N - q)-by-(N - q) block first */
|
|
Packit |
67cb25 |
v = gsl_matrix_submatrix(&m.matrix, q, q, N - q, N - q);
|
|
Packit |
67cb25 |
francis_schur_decomp(&v.matrix, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* operate on upper left q-by-q block */
|
|
Packit |
67cb25 |
v = gsl_matrix_submatrix(&m.matrix, 0, 0, q, q);
|
|
Packit |
67cb25 |
francis_schur_decomp(&v.matrix, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
N = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* handle special cases of N = 1 or 2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(&m.matrix, 0, 0), 0.0);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals, lambda1);
|
|
Packit |
67cb25 |
w->n_evals += 1;
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N == 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
francis_schur_standardize(&m.matrix, &lambda1, &lambda2, w);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals, lambda1);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(eval, w->n_evals + 1, lambda2);
|
|
Packit |
67cb25 |
w->n_evals += 2;
|
|
Packit |
67cb25 |
w->n_iter = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* francis_schur_decomp() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
francis_qrstep()
|
|
Packit |
67cb25 |
Perform a Francis QR step.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
See Golub & Van Loan, "Matrix Computations" (3rd ed),
|
|
Packit |
67cb25 |
algorithm 7.5.1
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: H - upper Hessenberg matrix
|
|
Packit |
67cb25 |
w - workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: The matrix H must be "reduced", ie: have no tiny subdiagonal
|
|
Packit |
67cb25 |
elements. When computing the first householder reflection,
|
|
Packit |
67cb25 |
we divide by H_{21} so it is necessary that this element
|
|
Packit |
67cb25 |
is not zero. When a subdiagonal element becomes negligible,
|
|
Packit |
67cb25 |
the calling function should call this routine with the
|
|
Packit |
67cb25 |
submatrices split by that element, so that we don't divide
|
|
Packit |
67cb25 |
by zeros.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline int
|
|
Packit |
67cb25 |
francis_qrstep(gsl_matrix * H, gsl_eigen_francis_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = H->size1;
|
|
Packit |
67cb25 |
size_t i; /* looping */
|
|
Packit |
67cb25 |
gsl_matrix_view m;
|
|
Packit |
67cb25 |
double tau_i; /* householder coefficient */
|
|
Packit |
67cb25 |
double dat[3]; /* householder vector */
|
|
Packit |
67cb25 |
double scale; /* scale factor to avoid overflow */
|
|
Packit |
67cb25 |
gsl_vector_view v2, v3;
|
|
Packit |
67cb25 |
size_t q, r;
|
|
Packit |
67cb25 |
size_t top = 0; /* location of H in original matrix */
|
|
Packit |
67cb25 |
double s,
|
|
Packit |
67cb25 |
disc;
|
|
Packit |
67cb25 |
double h_nn, /* H(n,n) */
|
|
Packit |
67cb25 |
h_nm1nm1, /* H(n-1,n-1) */
|
|
Packit |
67cb25 |
h_cross, /* H(n,n-1) * H(n-1,n) */
|
|
Packit |
67cb25 |
h_tmp1,
|
|
Packit |
67cb25 |
h_tmp2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v2 = gsl_vector_view_array(dat, 2);
|
|
Packit |
67cb25 |
v3 = gsl_vector_view_array(dat, 3);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ((w->n_iter % 10) == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* exceptional shifts: we have gone 10 iterations
|
|
Packit |
67cb25 |
* without finding a new eigenvalue, try a new choice of shifts.
|
|
Packit |
67cb25 |
* See LAPACK routine DLAHQR
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
s = fabs(gsl_matrix_get(H, N - 1, N - 2)) +
|
|
Packit |
67cb25 |
fabs(gsl_matrix_get(H, N - 2, N - 3));
|
|
Packit |
67cb25 |
h_nn = gsl_matrix_get(H, N - 1, N - 1) + GSL_FRANCIS_COEFF1 * s;
|
|
Packit |
67cb25 |
h_nm1nm1 = h_nn;
|
|
Packit |
67cb25 |
h_cross = GSL_FRANCIS_COEFF2 * s * s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* normal shifts - compute Rayleigh quotient and use
|
|
Packit |
67cb25 |
* Wilkinson shift if possible
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
h_nn = gsl_matrix_get(H, N - 1, N - 1);
|
|
Packit |
67cb25 |
h_nm1nm1 = gsl_matrix_get(H, N - 2, N - 2);
|
|
Packit |
67cb25 |
h_cross = gsl_matrix_get(H, N - 1, N - 2) *
|
|
Packit |
67cb25 |
gsl_matrix_get(H, N - 2, N - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
disc = 0.5 * (h_nm1nm1 - h_nn);
|
|
Packit |
67cb25 |
disc = disc * disc + h_cross;
|
|
Packit |
67cb25 |
if (disc > 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ave;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* real roots - use Wilkinson's shift twice */
|
|
Packit |
67cb25 |
disc = sqrt(disc);
|
|
Packit |
67cb25 |
ave = 0.5 * (h_nm1nm1 + h_nn);
|
|
Packit |
67cb25 |
if (fabs(h_nm1nm1) - fabs(h_nn) > 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
h_nm1nm1 = h_nm1nm1 * h_nn - h_cross;
|
|
Packit |
67cb25 |
h_nn = h_nm1nm1 / (disc * GSL_SIGN(ave) + ave);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
h_nn = disc * GSL_SIGN(ave) + ave;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
h_nm1nm1 = h_nn;
|
|
Packit |
67cb25 |
h_cross = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
h_tmp1 = h_nm1nm1 - gsl_matrix_get(H, 0, 0);
|
|
Packit |
67cb25 |
h_tmp2 = h_nn - gsl_matrix_get(H, 0, 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* These formulas are equivalent to those in Golub & Van Loan
|
|
Packit |
67cb25 |
* for the normal shift case - the terms have been rearranged
|
|
Packit |
67cb25 |
* to reduce possible roundoff error when subdiagonal elements
|
|
Packit |
67cb25 |
* are small
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
dat[0] = (h_tmp1*h_tmp2 - h_cross) / gsl_matrix_get(H, 1, 0) +
|
|
Packit |
67cb25 |
gsl_matrix_get(H, 0, 1);
|
|
Packit |
67cb25 |
dat[1] = gsl_matrix_get(H, 1, 1) - gsl_matrix_get(H, 0, 0) - h_tmp1 -
|
|
Packit |
67cb25 |
h_tmp2;
|
|
Packit |
67cb25 |
dat[2] = gsl_matrix_get(H, 2, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
scale = fabs(dat[0]) + fabs(dat[1]) + fabs(dat[2]);
|
|
Packit |
67cb25 |
if (scale != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* scale to prevent overflow or underflow */
|
|
Packit |
67cb25 |
dat[0] /= scale;
|
|
Packit |
67cb25 |
dat[1] /= scale;
|
|
Packit |
67cb25 |
dat[2] /= scale;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->Z || w->compute_t)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* get absolute indices of this (sub)matrix relative to the
|
|
Packit |
67cb25 |
* original Hessenberg matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
top = francis_get_submatrix(w->H, H);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N - 2; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tau_i = gsl_linalg_householder_transform(&v3.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau_i != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* q = max(1, i - 1) */
|
|
Packit |
67cb25 |
q = (1 > ((int)i - 1)) ? 0 : (i - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* r = min(i + 3, N - 1) */
|
|
Packit |
67cb25 |
r = ((i + 3) < (N - 1)) ? (i + 3) : (N - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->compute_t)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* We are computing the Schur form T, so we
|
|
Packit |
67cb25 |
* need to transform the whole matrix H
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* H -> P_k^t H P_k
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where P_k is the current Householder matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply left householder matrix (I - tau_i v v') to H */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(w->H,
|
|
Packit |
67cb25 |
top + i,
|
|
Packit |
67cb25 |
top + q,
|
|
Packit |
67cb25 |
3,
|
|
Packit |
67cb25 |
w->size - top - q);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm(tau_i, &v3.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply right householder matrix (I - tau_i v v') to H */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(w->H,
|
|
Packit |
67cb25 |
0,
|
|
Packit |
67cb25 |
top + i,
|
|
Packit |
67cb25 |
top + r + 1,
|
|
Packit |
67cb25 |
3);
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &v3.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* We are not computing the Schur form T, so we
|
|
Packit |
67cb25 |
* only need to transform the active block
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply left householder matrix (I - tau_i v v') to H */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(H, i, q, 3, N - q);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm(tau_i, &v3.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply right householder matrix (I - tau_i v v') to H */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(H, 0, i, r + 1, 3);
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &v3.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->Z)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* accumulate the similarity transformation into Z */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(w->Z, 0, top + i, w->size, 3);
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &v3.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* if (tau_i != 0.0) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
dat[0] = gsl_matrix_get(H, i + 1, i);
|
|
Packit |
67cb25 |
dat[1] = gsl_matrix_get(H, i + 2, i);
|
|
Packit |
67cb25 |
if (i < (N - 3))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
dat[2] = gsl_matrix_get(H, i + 3, i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
scale = fabs(dat[0]) + fabs(dat[1]) + fabs(dat[2]);
|
|
Packit |
67cb25 |
if (scale != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* scale to prevent overflow or underflow */
|
|
Packit |
67cb25 |
dat[0] /= scale;
|
|
Packit |
67cb25 |
dat[1] /= scale;
|
|
Packit |
67cb25 |
dat[2] /= scale;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* for (i = 0; i < N - 2; ++i) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
scale = fabs(dat[0]) + fabs(dat[1]);
|
|
Packit |
67cb25 |
if (scale != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* scale to prevent overflow or underflow */
|
|
Packit |
67cb25 |
dat[0] /= scale;
|
|
Packit |
67cb25 |
dat[1] /= scale;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tau_i = gsl_linalg_householder_transform(&v2.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->compute_t)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(w->H,
|
|
Packit |
67cb25 |
top + N - 2,
|
|
Packit |
67cb25 |
top + N - 3,
|
|
Packit |
67cb25 |
2,
|
|
Packit |
67cb25 |
w->size - top - N + 3);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm(tau_i, &v2.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(w->H,
|
|
Packit |
67cb25 |
0,
|
|
Packit |
67cb25 |
top + N - 2,
|
|
Packit |
67cb25 |
top + N,
|
|
Packit |
67cb25 |
2);
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &v2.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm(tau_i, &v2.vector, &m.matrix);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &v2.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->Z)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* accumulate transformation into Z */
|
|
Packit |
67cb25 |
m = gsl_matrix_submatrix(w->Z, 0, top + N - 2, w->size, 2);
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh(tau_i, &v2.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
} /* francis_qrstep() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
francis_search_subdiag_small_elements()
|
|
Packit |
67cb25 |
Search for a small subdiagonal element starting from the bottom
|
|
Packit |
67cb25 |
of a matrix A. A small element is one that satisfies:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|A_{i,i-1}| <= eps * (|A_{i,i}| + |A_{i-1,i-1}|)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - matrix (must be at least 3-by-3)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: row index of small subdiagonal element or 0 if not found
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: the first small element that is found (starting from bottom)
|
|
Packit |
67cb25 |
is set to zero
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline size_t
|
|
Packit |
67cb25 |
francis_search_subdiag_small_elements(gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
double dpel = gsl_matrix_get(A, N - 2, N - 2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = N - 1; i > 0; --i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sel = gsl_matrix_get(A, i, i - 1);
|
|
Packit |
67cb25 |
double del = gsl_matrix_get(A, i, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ((sel == 0.0) ||
|
|
Packit |
67cb25 |
(fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel))))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(A, i, i - 1, 0.0);
|
|
Packit |
67cb25 |
return (i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
dpel = del;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (0);
|
|
Packit |
67cb25 |
} /* francis_search_subdiag_small_elements() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
francis_schur_standardize()
|
|
Packit |
67cb25 |
Convert a 2-by-2 diagonal block in the Schur form to standard form
|
|
Packit |
67cb25 |
and update the rest of T and Z matrices if required.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - 2-by-2 matrix
|
|
Packit |
67cb25 |
eval1 - where to store eigenvalue 1
|
|
Packit |
67cb25 |
eval2 - where to store eigenvalue 2
|
|
Packit |
67cb25 |
w - francis workspace
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline void
|
|
Packit |
67cb25 |
francis_schur_standardize(gsl_matrix *A, gsl_complex *eval1,
|
|
Packit |
67cb25 |
gsl_complex *eval2,
|
|
Packit |
67cb25 |
gsl_eigen_francis_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = w->size;
|
|
Packit |
67cb25 |
double cs, sn;
|
|
Packit |
67cb25 |
size_t top;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* figure out where the submatrix A resides in the
|
|
Packit |
67cb25 |
* original matrix H
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
top = francis_get_submatrix(w->H, A);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* convert 2-by-2 block to standard form */
|
|
Packit |
67cb25 |
francis_standard_form(A, &cs, &sn;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* set eigenvalues */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_SET_REAL(eval1, gsl_matrix_get(A, 0, 0));
|
|
Packit |
67cb25 |
GSL_SET_REAL(eval2, gsl_matrix_get(A, 1, 1));
|
|
Packit |
67cb25 |
if (gsl_matrix_get(A, 1, 0) == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_SET_IMAG(eval1, 0.0);
|
|
Packit |
67cb25 |
GSL_SET_IMAG(eval2, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double tmp = sqrt(fabs(gsl_matrix_get(A, 0, 1)) *
|
|
Packit |
67cb25 |
fabs(gsl_matrix_get(A, 1, 0)));
|
|
Packit |
67cb25 |
GSL_SET_IMAG(eval1, tmp);
|
|
Packit |
67cb25 |
GSL_SET_IMAG(eval2, -tmp);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->compute_t)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view xv, yv;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* The above call to francis_standard_form transformed a 2-by-2 block
|
|
Packit |
67cb25 |
* of T into upper triangular form via the transformation
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* U = [ CS -SN ]
|
|
Packit |
67cb25 |
* [ SN CS ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The original matrix T was
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* T = [ T_{11} | T_{12} | T_{13} ]
|
|
Packit |
67cb25 |
* [ 0* | A | T_{23} ]
|
|
Packit |
67cb25 |
* [ 0 | 0* | T_{33} ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where 0* indicates all zeros except for possibly
|
|
Packit |
67cb25 |
* one subdiagonal element next to A.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* After francis_standard_form, T looks like this:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* T = [ T_{11} | T_{12} | T_{13} ]
|
|
Packit |
67cb25 |
* [ 0* | U^t A U | T_{23} ]
|
|
Packit |
67cb25 |
* [ 0 | 0* | T_{33} ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* since only the 2-by-2 block of A was changed. However,
|
|
Packit |
67cb25 |
* in order to be able to back transform T at the end,
|
|
Packit |
67cb25 |
* we need to apply the U transformation to the rest
|
|
Packit |
67cb25 |
* of the matrix T since there is no way to apply a
|
|
Packit |
67cb25 |
* similarity transformation to T and change only the
|
|
Packit |
67cb25 |
* middle 2-by-2 block. In other words, let
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* M = [ I 0 0 ]
|
|
Packit |
67cb25 |
* [ 0 U 0 ]
|
|
Packit |
67cb25 |
* [ 0 0 I ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* and compute
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* M^t T M = [ T_{11} | T_{12} U | T_{13} ]
|
|
Packit |
67cb25 |
* [ U^t 0* | U^t A U | U^t T_{23} ]
|
|
Packit |
67cb25 |
* [ 0 | 0* U | T_{33} ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* So basically we need to apply the transformation U
|
|
Packit |
67cb25 |
* to the i x 2 matrix T_{12} and the 2 x (n - i + 2)
|
|
Packit |
67cb25 |
* matrix T_{23}, where i is the index of the top of A
|
|
Packit |
67cb25 |
* in T.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The BLAS routine drot() is suited for this.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (top < (N - 2))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* transform the 2 rows of T_{23} */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xv = gsl_matrix_subrow(w->H, top, top + 2, N - top - 2);
|
|
Packit |
67cb25 |
yv = gsl_matrix_subrow(w->H, top + 1, top + 2, N - top - 2);
|
|
Packit |
67cb25 |
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (top > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* transform the 2 columns of T_{12} */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xv = gsl_matrix_subcolumn(w->H, top, 0, top);
|
|
Packit |
67cb25 |
yv = gsl_matrix_subcolumn(w->H, top + 1, 0, top);
|
|
Packit |
67cb25 |
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* if (w->compute_t) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (w->Z)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view xv, yv;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Accumulate the transformation in Z. Here, Z -> Z * M
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* So:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Z -> [ Z_{11} | Z_{12} U | Z_{13} ]
|
|
Packit |
67cb25 |
* [ Z_{21} | Z_{22} U | Z_{23} ]
|
|
Packit |
67cb25 |
* [ Z_{31} | Z_{32} U | Z_{33} ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* So we just need to apply drot() to the 2 columns
|
|
Packit |
67cb25 |
* starting at index 'top'
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xv = gsl_matrix_column(w->Z, top);
|
|
Packit |
67cb25 |
yv = gsl_matrix_column(w->Z, top + 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
|
|
Packit |
67cb25 |
} /* if (w->Z) */
|
|
Packit |
67cb25 |
} /* francis_schur_standardize() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
francis_get_submatrix()
|
|
Packit |
67cb25 |
B is a submatrix of A. The goal of this function is to
|
|
Packit |
67cb25 |
compute the indices in A of where the matrix B resides
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline size_t
|
|
Packit |
67cb25 |
francis_get_submatrix(gsl_matrix *A, gsl_matrix *B)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t diff;
|
|
Packit |
67cb25 |
double ratio;
|
|
Packit |
67cb25 |
size_t top;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
diff = (size_t) (B->data - A->data);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ratio = (double)diff / ((double) (A->tda + 1));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
top = (size_t) floor(ratio);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return top;
|
|
Packit |
67cb25 |
} /* francis_get_submatrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
francis_standard_form()
|
|
Packit |
67cb25 |
Compute the Schur factorization of a real 2-by-2 matrix in
|
|
Packit |
67cb25 |
standard form:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ A B ] = [ CS -SN ] [ T11 T12 ] [ CS SN ]
|
|
Packit |
67cb25 |
[ C D ] [ SN CS ] [ T21 T22 ] [-SN CS ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where either:
|
|
Packit |
67cb25 |
1) T21 = 0 so that T11 and T22 are real eigenvalues of the matrix, or
|
|
Packit |
67cb25 |
2) T11 = T22 and T21*T12 < 0, so that T11 +/- sqrt(|T21*T12|) are
|
|
Packit |
67cb25 |
complex conjugate eigenvalues
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - 2-by-2 matrix
|
|
Packit |
67cb25 |
cs - where to store cosine parameter of rotation matrix
|
|
Packit |
67cb25 |
sn - where to store sine parameter of rotation matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: 1) based on LAPACK routine DLANV2
|
|
Packit |
67cb25 |
2) On output, A is modified to contain the matrix in standard form
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
francis_standard_form(gsl_matrix *A, double *cs, double *sn)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double a, b, c, d; /* input matrix values */
|
|
Packit |
67cb25 |
double tmp;
|
|
Packit |
67cb25 |
double p, z;
|
|
Packit |
67cb25 |
double bcmax, bcmis, scale;
|
|
Packit |
67cb25 |
double tau, sigma;
|
|
Packit |
67cb25 |
double cs1, sn1;
|
|
Packit |
67cb25 |
double aa, bb, cc, dd;
|
|
Packit |
67cb25 |
double sab, sac;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
a = gsl_matrix_get(A, 0, 0);
|
|
Packit |
67cb25 |
b = gsl_matrix_get(A, 0, 1);
|
|
Packit |
67cb25 |
c = gsl_matrix_get(A, 1, 0);
|
|
Packit |
67cb25 |
d = gsl_matrix_get(A, 1, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (c == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* matrix is already upper triangular - set rotation matrix
|
|
Packit |
67cb25 |
* to the identity
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
*cs = 1.0;
|
|
Packit |
67cb25 |
*sn = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (b == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* swap rows and columns to make it upper triangular */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*cs = 0.0;
|
|
Packit |
67cb25 |
*sn = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tmp = d;
|
|
Packit |
67cb25 |
d = a;
|
|
Packit |
67cb25 |
a = tmp;
|
|
Packit |
67cb25 |
b = -c;
|
|
Packit |
67cb25 |
c = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (((a - d) == 0.0) && (GSL_SIGN(b) != GSL_SIGN(c)))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* the matrix has complex eigenvalues with a == d */
|
|
Packit |
67cb25 |
*cs = 1.0;
|
|
Packit |
67cb25 |
*sn = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tmp = a - d;
|
|
Packit |
67cb25 |
p = 0.5 * tmp;
|
|
Packit |
67cb25 |
bcmax = GSL_MAX(fabs(b), fabs(c));
|
|
Packit |
67cb25 |
bcmis = GSL_MIN(fabs(b), fabs(c)) * GSL_SIGN(b) * GSL_SIGN(c);
|
|
Packit |
67cb25 |
scale = GSL_MAX(fabs(p), bcmax);
|
|
Packit |
67cb25 |
z = (p / scale) * p + (bcmax / scale) * bcmis;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (z >= 4.0 * GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* real eigenvalues, compute a and d */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
z = p + GSL_SIGN(p) * fabs(sqrt(scale) * sqrt(z));
|
|
Packit |
67cb25 |
a = d + z;
|
|
Packit |
67cb25 |
d -= (bcmax / z) * bcmis;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute b and the rotation matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tau = gsl_hypot(c, z);
|
|
Packit |
67cb25 |
*cs = z / tau;
|
|
Packit |
67cb25 |
*sn = c / tau;
|
|
Packit |
67cb25 |
b -= c;
|
|
Packit |
67cb25 |
c = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* complex eigenvalues, or real (almost) equal eigenvalues -
|
|
Packit |
67cb25 |
* make diagonal elements equal
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sigma = b + c;
|
|
Packit |
67cb25 |
tau = gsl_hypot(sigma, tmp);
|
|
Packit |
67cb25 |
*cs = sqrt(0.5 * (1.0 + fabs(sigma) / tau));
|
|
Packit |
67cb25 |
*sn = -(p / (tau * (*cs))) * GSL_SIGN(sigma);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Compute [ AA BB ] = [ A B ] [ CS -SN ]
|
|
Packit |
67cb25 |
* [ CC DD ] [ C D ] [ SN CS ]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
aa = a * (*cs) + b * (*sn);
|
|
Packit |
67cb25 |
bb = -a * (*sn) + b * (*cs);
|
|
Packit |
67cb25 |
cc = c * (*cs) + d * (*sn);
|
|
Packit |
67cb25 |
dd = -c * (*sn) + d * (*cs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Compute [ A B ] = [ CS SN ] [ AA BB ]
|
|
Packit |
67cb25 |
* [ C D ] [-SN CS ] [ CC DD ]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
a = aa * (*cs) + cc * (*sn);
|
|
Packit |
67cb25 |
b = bb * (*cs) + dd * (*sn);
|
|
Packit |
67cb25 |
c = -aa * (*sn) + cc * (*cs);
|
|
Packit |
67cb25 |
d = -bb * (*sn) + dd * (*cs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tmp = 0.5 * (a + d);
|
|
Packit |
67cb25 |
a = d = tmp;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (c != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (b != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (GSL_SIGN(b) == GSL_SIGN(c))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* real eigenvalues: reduce to upper triangular
|
|
Packit |
67cb25 |
* form
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
sab = sqrt(fabs(b));
|
|
Packit |
67cb25 |
sac = sqrt(fabs(c));
|
|
Packit |
67cb25 |
p = GSL_SIGN(c) * fabs(sab * sac);
|
|
Packit |
67cb25 |
tau = 1.0 / sqrt(fabs(b + c));
|
|
Packit |
67cb25 |
a = tmp + p;
|
|
Packit |
67cb25 |
d = tmp - p;
|
|
Packit |
67cb25 |
b -= c;
|
|
Packit |
67cb25 |
c = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
cs1 = sab * tau;
|
|
Packit |
67cb25 |
sn1 = sac * tau;
|
|
Packit |
67cb25 |
tmp = (*cs) * cs1 - (*sn) * sn1;
|
|
Packit |
67cb25 |
*sn = (*cs) * sn1 + (*sn) * cs1;
|
|
Packit |
67cb25 |
*cs = tmp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
b = -c;
|
|
Packit |
67cb25 |
c = 0.0;
|
|
Packit |
67cb25 |
tmp = *cs;
|
|
Packit |
67cb25 |
*cs = -(*sn);
|
|
Packit |
67cb25 |
*sn = tmp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* set new matrix elements */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set(A, 0, 0, a);
|
|
Packit |
67cb25 |
gsl_matrix_set(A, 0, 1, b);
|
|
Packit |
67cb25 |
gsl_matrix_set(A, 1, 0, c);
|
|
Packit |
67cb25 |
gsl_matrix_set(A, 1, 1, d);
|
|
Packit |
67cb25 |
} /* francis_standard_form() */
|