/* 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);

/*

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)

{

  RETURN_IF_NULL (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

            {

              /* 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, 2, 2, N, N);

          r = gsl_matrix_submatrix(&r.matrix, 2, 2, N, N);

        }

      else

        {

          /*

           * There is a zero element on the subdiagonal somewhere

           * in the middle of the matrix - we can now operate

           * separately on the two submatrices split by this

           * element. q is the row index of the zero element.

           */

          /* operate on lower right (N - q)-by-(N - q) block first */

          vh = gsl_matrix_submatrix(&h.matrix, q, q, N - q, N - q);

          vr = gsl_matrix_submatrix(&r.matrix, q, q, N - q, N - q);

          gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);

          /* operate on upper left q-by-q block */

          vh = gsl_matrix_submatrix(&h.matrix, 0, 0, q, q);

          vr = gsl_matrix_submatrix(&r.matrix, 0, 0, q, q);

          gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);

          N = 0;

        }

    } /* while ((N > 1) && (w->n_iter)++ < w->max_iterations) */

  /* handle special case of N = 1 */

  if (N == 1)

    {

      gen_schur_standardize1(&h.matrix, &r.matrix, &a, &b, w);

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

    }

} /* gen_schur_decomp() */

/*

gen_qzstep()

  This routine determines what type of QZ step to perform on

the generalized matrix pair (H, R). If the pair is 3-by-3 or bigger,

we look at the bottom right 2-by-2 block. If this block has complex

eigenvalues, we perform a Francis double shift QZ sweep. If it

has real eigenvalues, we perform an implicit single shift QZ sweep.

If the pair is 2-by-2 with real eigenvalues, we perform a single

shift sweep. If it has complex eigenvalues, we return GSL_CONTINUE

to notify the calling function that a 2-by-2 block with complex

eigenvalues has converged, so that it may then call

gen_schur_standardize2(). In the real eigenvalue case, we want to

continue doing QZ sweeps to break it up into two 1-by-1 blocks.

See LAPACK routine DHGEQZ and [1] for more information.

Inputs: H - upper Hessenberg matrix (at least 2-by-2)

        R - upper triangular matrix (at least 2-by-2)

        w - workspace

Return: GSL_SUCCESS on normal completion

        GSL_CONTINUE if we detect a 2-by-2 block with complex eigenvalues

*/

static inline int

gen_qzstep(gsl_matrix *H, gsl_matrix *R, gsl_eigen_gen_workspace *w)

{

  const size_t N = H->size1;

  gsl_matrix_view vh, vr; /* views of bottom right 2-by-2 block */

  double wr1, wr2, wi;

  double scale1, scale2, scale;

  double cs, sn;          /* givens rotation */

  double temp,            /* temporary variables */

         temp2;

  size_t j;               /* looping */

  gsl_vector_view xv, yv; /* temporary views */

  size_t top = 0;

  size_t rows;

  if (w->n_iter % 10 == 0)

    {

      /*

       * Exceptional shift - we have gone 10 iterations without finding

       * a new eigenvalue, do a single shift sweep with an

       * exceptional shift

       */

      if ((GSL_DBL_MIN * w->max_iterations) *

          fabs(gsl_matrix_get(H, N - 2, N - 1)) <

          fabs(gsl_matrix_get(R, N - 2, N - 2)))

        {

          w->eshift += gsl_matrix_get(H, N - 2, N - 1) /

                       gsl_matrix_get(R, N - 2, N - 2);

        }

      else

        w->eshift += 1.0 / (GSL_DBL_MIN * w->max_iterations);

      if ((w->eshift < GSL_DBL_EPSILON) &&

          (GSL_DBL_MIN * w->max_iterations) *

          fabs(gsl_matrix_get(H, N - 1, N - 2)) <

          fabs(gsl_matrix_get(R, N - 2, N - 2)))

        {

          w->eshift = GEN_ESHIFT_COEFF *

                      (w->ascale * gsl_matrix_get(H, N - 1, N - 2)) /

                      (w->bscale * gsl_matrix_get(R, N - 2, N - 2));

        }

      scale1 = 1.0;

      wr1 = w->eshift;

    }

  else

    {

      /*

       * Compute generalized eigenvalues of bottom right 2-by-2 block

       * to be used as shifts - wr1 is the Wilkinson shift

       */

      vh = gsl_matrix_submatrix(H, N - 2, N - 2, 2, 2);

      vr = gsl_matrix_submatrix(R, N - 2, N - 2, 2, 2);

      gsl_schur_gen_eigvals(&vh.matrix,

                            &vr.matrix,

                            &wr1,

                            &wr2,

                            &wi,

                            &scale1,

                            &scale2);

      if (wi != 0.0)

        {

          /* complex eigenvalues */

          if (N == 2)

            {

              /*

               * its a 2-by-2 block with complex eigenvalues - notify

               * the calling function to deflate

               */

              return (GSL_CONTINUE);

            }

          else

            {

              /* do a francis double shift sweep */

              gen_qzstep_d(H, R, w);

            }

          return GSL_SUCCESS;

        }

    }

  /* real eigenvalues - perform single shift QZ step */

  temp = GSL_MIN(w->ascale, 1.0) * (0.5 / GSL_DBL_MIN);

  if (scale1 > temp)

    scale = temp / scale1;

  else

    scale = 1.0;

  temp = GSL_MIN(w->bscale, 1.0) * (0.5 / GSL_DBL_MIN);

  if (fabs(wr1) > temp)

    scale = GSL_MIN(scale, temp / fabs(wr1));

  scale1 *= scale;

  wr1 *= scale;

  if (w->needtop)

    {

      /* get absolute index of this matrix relative to original matrix */

      top = gen_get_submatrix(w->H, H);

    }

  temp = scale1*gsl_matrix_get(H, 0, 0) - wr1*gsl_matrix_get(R, 0, 0);

  temp2 = scale1*gsl_matrix_get(H, 1, 0);

  gsl_linalg_givens(temp, temp2, &cs, &sn;;

  sn = -sn;

  for (j = 0; j < N - 1; ++j)

    {

      if (j > 0)

        {

          temp = gsl_matrix_get(H, j, j - 1);

          temp2 = gsl_matrix_get(H, j + 1, j - 1);

          gsl_linalg_givens(temp, temp2, &cs, &sn;;

          sn = -sn;

          /* apply to column (j - 1) */

          temp = cs * gsl_matrix_get(H, j, j - 1) +

                 sn * gsl_matrix_get(H, j + 1, j - 1);

          gsl_matrix_set(H, j, j - 1, temp);

          gsl_matrix_set(H, j + 1, j - 1, 0.0);

        }

      /* apply G to H(j:j+1,:) and T(j:j+1,:) */

      if (w->compute_s)

        {

          xv = gsl_matrix_subrow(w->H, top + j, top + j, w->size - top - j);

          yv = gsl_matrix_subrow(w->H, top + j + 1, top + j, w->size - top - j);

        }

      else

        {

          xv = gsl_matrix_subrow(H, j, j, N - j);

          yv = gsl_matrix_subrow(H, j + 1, j, N - j);

        }

      gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

      if (w->compute_t)

        {

          xv = gsl_matrix_subrow(w->R, top + j, top + j, w->size - top - j);

          yv = gsl_matrix_subrow(w->R, top + j + 1, top + j, w->size - top - j);

        }

      else

        {

          xv = gsl_matrix_subrow(R, j, j, N - j);

          yv = gsl_matrix_subrow(R, j + 1, j, N - j);

        }

      gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

      if (w->Q)

        {

          /* accumulate Q: Q -> QG */

          xv = gsl_matrix_column(w->Q, top + j);

          yv = gsl_matrix_column(w->Q, top + j + 1);

          gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

        }

      temp = gsl_matrix_get(R, j + 1, j + 1);

      temp2 = gsl_matrix_get(R, j + 1, j);

      gsl_linalg_givens(temp, temp2, &cs, &sn;;

      rows = GSL_MIN(j + 3, N);

      if (w->compute_s)

        {

          xv = gsl_matrix_subcolumn(w->H, top + j, 0, top + rows);

          yv = gsl_matrix_subcolumn(w->H, top + j + 1, 0, top + rows);

        }

      else

        {

          xv = gsl_matrix_subcolumn(H, j, 0, rows);

          yv = gsl_matrix_subcolumn(H, j + 1, 0, rows);

        }

      gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

      rows = GSL_MIN(j + 2, N);

      if (w->compute_t)

        {

          xv = gsl_matrix_subcolumn(w->R, top + j, 0, top + rows);

          yv = gsl_matrix_subcolumn(w->R, top + j + 1, 0, top + rows);

        }

      else

        {

          xv = gsl_matrix_subcolumn(R, j, 0, rows);

          yv = gsl_matrix_subcolumn(R, j + 1, 0, rows);

        }

      gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

      if (w->Z)

        {

          /* accumulate Z: Z -> ZG */

          xv = gsl_matrix_column(w->Z, top + j);

          yv = gsl_matrix_column(w->Z, top + j + 1);

          gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

        }

    } /* for (j = 0; j < N - 1; ++j) */

  return GSL_SUCCESS;

} /* gen_qzstep() */

/*

gen_qzstep_d()

  Perform an implicit double shift QZ step.

See Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm 7.7.2

Inputs: H - upper Hessenberg matrix (at least 3-by-3)

        R - upper triangular matrix (at least 3-by-3)

        w - workspace

*/

static inline void

gen_qzstep_d(gsl_matrix *H, gsl_matrix *R, gsl_eigen_gen_workspace *w)

{

  const size_t N = H->size1;

  size_t j;               /* looping */

  double dat[3];          /* householder vector */

  double tau;             /* householder coefficient */

  gsl_vector_view v2, v3; /* views into 'dat' */

  gsl_matrix_view m;      /* temporary view */

  double tmp;

  size_t q, r;

  size_t top = 0;         /* location of H in original matrix */

  double scale;

  double AB11,            /* various matrix element ratios */

         AB22,

         ABNN,

         ABMM,

         AMNBNN,

         ANMBMM,

         A21B11,

         A12B22,

         A32B22,

         B12B22,

         BMNBNN;

  v2 = gsl_vector_view_array(dat, 2);

  v3 = gsl_vector_view_array(dat, 3);

  if (w->needtop)

    {

      /* get absolute index of this matrix relative to original matrix */

      top = gen_get_submatrix(w->H, H);

    }

  /*

   * Similar to the QR method, we take the shifts to be the two

   * zeros of the problem

   *

   * det[H(n-1:n,n-1:n) - s*R(n-1:n,n-1:n)] = 0

   *

   * The initial householder vector elements are then given by

   * Eq. 4.1 of [1], which are designed to reduce errors when

   * off diagonal elements are small.

   */

  ABMM = (w->ascale * gsl_matrix_get(H, N - 2, N - 2)) /

         (w->bscale * gsl_matrix_get(R, N - 2, N - 2));

  ABNN = (w->ascale * gsl_matrix_get(H, N - 1, N - 1)) /

         (w->bscale * gsl_matrix_get(R, N - 1, N - 1));

  AB11 = (w->ascale * gsl_matrix_get(H, 0, 0)) /

         (w->bscale * gsl_matrix_get(R, 0, 0));

  AB22 = (w->ascale * gsl_matrix_get(H, 1, 1)) /

         (w->bscale * gsl_matrix_get(R, 1, 1));

  AMNBNN = (w->ascale * gsl_matrix_get(H, N - 2, N - 1)) /

           (w->bscale * gsl_matrix_get(R, N - 1, N - 1));

  ANMBMM = (w->ascale * gsl_matrix_get(H, N - 1, N - 2)) /

           (w->bscale * gsl_matrix_get(R, N - 2, N - 2));

  BMNBNN = gsl_matrix_get(R, N - 2, N - 1) /

           gsl_matrix_get(R, N - 1, N - 1);

  A21B11 = (w->ascale * gsl_matrix_get(H, 1, 0)) /

           (w->bscale * gsl_matrix_get(R, 0, 0));

  A12B22 = (w->ascale * gsl_matrix_get(H, 0, 1)) /

           (w->bscale * gsl_matrix_get(R, 1, 1));

  A32B22 = (w->ascale * gsl_matrix_get(H, 2, 1)) /