htswanalysis/MACS/lib/gsl/gsl-1.11/eigen/gen.c

   1 /* eigen/gen.c
   2  *
   3  * Copyright (C) 2006, 2007 Patrick Alken
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <stdlib.h>
  21 #include <math.h>
  22
  23 #include <config.h>
  24 #include <gsl/gsl_eigen.h>
  25 #include <gsl/gsl_linalg.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_blas.h>
  28 #include <gsl/gsl_vector.h>
  29 #include <gsl/gsl_vector_complex.h>
  30 #include <gsl/gsl_matrix.h>
  31
  32 /*
  33  * This module computes the eigenvalues of a real generalized
  34  * eigensystem A x = \lambda B x. Left and right Schur vectors
  35  * are optionally computed as well.
  36  *
  37  * Based on the algorithm from Moler and Stewart
  38  * [1] C. Moler, G. Stewart, "An Algorithm for Generalized Matrix
  39  *     Eigenvalue Problems", SIAM J. Numer. Anal., Vol 10, No 2, 1973.
  40  *
  41  * This algorithm is also described in the book
  42  * [2] Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm 7.7.3
  43  *
  44  * This file contains routines based on original code from LAPACK
  45  * which is distributed under the modified BSD license.
  46  */
  47
  48 #define GEN_ESHIFT_COEFF     (1.736)
  49
  50 static void gen_schur_decomp(gsl_matrix *H, gsl_matrix *R,
  51                              gsl_vector_complex *alpha, gsl_vector *beta,
  52                              gsl_eigen_gen_workspace *w);
  53 static inline int gen_qzstep(gsl_matrix *H, gsl_matrix *R,
  54                              gsl_eigen_gen_workspace *w);
  55 static inline void gen_qzstep_d(gsl_matrix *H, gsl_matrix *R,
  56                                 gsl_eigen_gen_workspace *w);
  57 static void gen_tri_split_top(gsl_matrix *H, gsl_matrix *R,
  58                               gsl_eigen_gen_workspace *w);
  59 static inline void gen_tri_chase_zero(gsl_matrix *H, gsl_matrix *R,
  60                                       size_t q,
  61                                       gsl_eigen_gen_workspace *w);
  62 static inline void gen_tri_zero_H(gsl_matrix *H, gsl_matrix *R,
  63                                   gsl_eigen_gen_workspace *w);
  64 static inline size_t gen_search_small_elements(gsl_matrix *H,
  65                                                gsl_matrix *R,
  66                                                int *flag,
  67                                                gsl_eigen_gen_workspace *w);
  68 static int gen_schur_standardize1(gsl_matrix *A, gsl_matrix *B,
  69                                   double *alphar, double *beta,
  70                                   gsl_eigen_gen_workspace *w);
  71 static int gen_schur_standardize2(gsl_matrix *A, gsl_matrix *B,
  72                                   gsl_complex *alpha1,
  73                                   gsl_complex *alpha2,
  74                                   double *beta1, double *beta2,
  75                                   gsl_eigen_gen_workspace *w);
  76 static int gen_compute_eigenvals(gsl_matrix *A, gsl_matrix *B,
  77                                  gsl_complex *alpha1,
  78                                  gsl_complex *alpha2, double *beta1,
  79                                  double *beta2);
  80 static void gen_store_eigval1(const gsl_matrix *H, const double a,
  81                               const double b, gsl_vector_complex *alpha,
  82                               gsl_vector *beta,
  83                               gsl_eigen_gen_workspace *w);
  84 static void gen_store_eigval2(const gsl_matrix *H,
  85                               const gsl_complex *alpha1,
  86                               const double beta1,
  87                               const gsl_complex *alpha2,
  88                               const double beta2,
  89                               gsl_vector_complex *alpha,
  90                               gsl_vector *beta,
  91                               gsl_eigen_gen_workspace *w);
  92 static inline size_t gen_get_submatrix(const gsl_matrix *A,
  93                                        const gsl_matrix *B);
  94
  95 /*FIX**/
  96 inline static double normF (gsl_matrix * A);
  97 inline static void create_givens (const double a, const double b, double *c, double *s);
  98
  99 /*
 100 gsl_eigen_gen_alloc()
 101
 102 Allocate a workspace for solving the generalized eigenvalue problem.
 103 The size of this workspace is O(n)
 104
 105 Inputs: n - size of matrices
 106
 107 Return: pointer to workspace
 108 */
 109
 110 gsl_eigen_gen_workspace *
 111 gsl_eigen_gen_alloc(const size_t n)
 112 {
 113   gsl_eigen_gen_workspace *w;
 114
 115   if (n == 0)
 116     {
 117       GSL_ERROR_NULL ("matrix dimension must be positive integer",
 118                       GSL_EINVAL);
 119     }
 120
 121   w = calloc (1, sizeof (gsl_eigen_gen_workspace));
 122
 123   if (w == 0)
 124     {
 125       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
 126     }
 127
 128   w->size = n;
 129   w->max_iterations = 30 * n;
 130   w->n_evals = 0;
 131   w->n_iter = 0;
 132   w->needtop = 0;
 133   w->atol = 0.0;
 134   w->btol = 0.0;
 135   w->ascale = 0.0;
 136   w->bscale = 0.0;
 137   w->eshift = 0.0;
 138   w->H = NULL;
 139   w->R = NULL;
 140   w->compute_s = 0;
 141   w->compute_t = 0;
 142   w->Q = NULL;
 143   w->Z = NULL;
 144
 145   w->work = gsl_vector_alloc(n);
 146
 147   if (w->work == 0)
 148     {
 149       gsl_eigen_gen_free(w);
 150       GSL_ERROR_NULL ("failed to allocate space for additional workspace", GSL_ENOMEM);
 151     }
 152
 153   return (w);
 154 } /* gsl_eigen_gen_alloc() */
 155
 156 /*
 157 gsl_eigen_gen_free()
 158   Free workspace w
 159 */
 160
 161 void
 162 gsl_eigen_gen_free (gsl_eigen_gen_workspace * w)
 163 {
 164   if (w->work)
 165     gsl_vector_free(w->work);
 166
 167   free(w);
 168 } /* gsl_eigen_gen_free() */
 169
 170 /*
 171 gsl_eigen_gen_params()
 172   Set parameters which define how we solve the eigenvalue problem
 173
 174 Inputs: compute_s - 1 if we want to compute S, 0 if not
 175         compute_t - 1 if we want to compute T, 0 if not
 176         balance   - 1 if we want to balance matrices, 0 if not
 177         w         - gen workspace
 178
 179 Return: none
 180 */
 181
 182 void
 183 gsl_eigen_gen_params (const int compute_s, const int compute_t,
 184                       const int balance, gsl_eigen_gen_workspace *w)
 185 {
 186   w->compute_s = compute_s;
 187   w->compute_t = compute_t;
 188 } /* gsl_eigen_gen_params() */
 189
 190 /*
 191 gsl_eigen_gen()
 192
 193 Solve the generalized eigenvalue problem
 194
 195 A x = \lambda B x
 196
 197 for the eigenvalues \lambda.
 198
 199 Inputs: A     - general real matrix
 200         B     - general real matrix
 201         alpha - where to store eigenvalue numerators
 202         beta  - where to store eigenvalue denominators
 203         w     - workspace
 204
 205 Return: success or error
 206 */
 207
 208 int
 209 gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B, gsl_vector_complex * alpha,
 210                gsl_vector * beta, gsl_eigen_gen_workspace * w)
 211 {
 212   const size_t N = A->size1;
 213
 214   /* check matrix and vector sizes */
 215
 216   if (N != A->size2)
 217     {
 218       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 219     }
 220   else if ((N != B->size1) || (N != B->size2))
 221     {
 222       GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
 223     }
 224   else if (alpha->size != N || beta->size != N)
 225     {
 226       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 227     }
 228   else if (w->size != N)
 229     {
 230       GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
 231     }
 232   else
 233     {
 234       double anorm, bnorm;
 235
 236       /* compute the Hessenberg-Triangular reduction of (A, B) */
 237       gsl_linalg_hesstri_decomp(A, B, w->Q, w->Z, w->work);
 238
 239       /* save pointers to original matrices */
 240       w->H = A;
 241       w->R = B;
 242
 243       w->n_evals = 0;
 244       w->n_iter = 0;
 245       w->eshift = 0.0;
 246
 247       /* determine if we need to compute top indices in QZ step */
 248       w->needtop = w->Q != 0 || w->Z != 0 || w->compute_t || w->compute_s;
 249
 250       /* compute matrix norms */
 251       anorm = normF(A);
 252       bnorm = normF(B);
 253
 254       /* compute tolerances and scaling factors */
 255       w->atol = GSL_MAX(GSL_DBL_MIN, GSL_DBL_EPSILON * anorm);
 256       w->btol = GSL_MAX(GSL_DBL_MIN, GSL_DBL_EPSILON * bnorm);
 257       w->ascale = 1.0 / GSL_MAX(GSL_DBL_MIN, anorm);
 258       w->bscale = 1.0 / GSL_MAX(GSL_DBL_MIN, bnorm);
 259
 260       /* compute the generalized Schur decomposition and eigenvalues */
 261       gen_schur_decomp(A, B, alpha, beta, w);
 262
 263       if (w->n_evals != N)
 264         return GSL_EMAXITER;
 265
 266       return GSL_SUCCESS;
 267     }
 268 } /* gsl_eigen_gen() */
 269
 270 /*
 271 gsl_eigen_gen_QZ()
 272
 273 Solve the generalized eigenvalue problem
 274
 275 A x = \lambda B x
 276
 277 for the eigenvalues \lambda. Optionally compute left and/or right
 278 Schur vectors Q and Z which satisfy:
 279
 280 A = Q S Z^t
 281 B = Q T Z^t
 282
 283 where (S, T) is the generalized Schur form of (A, B)
 284
 285 Inputs: A     - general real matrix
 286         B     - general real matrix
 287         alpha - where to store eigenvalue numerators
 288         beta  - where to store eigenvalue denominators
 289         Q     - if non-null, where to store left Schur vectors
 290         Z     - if non-null, where to store right Schur vectors
 291         w     - workspace
 292
 293 Return: success or error
 294 */
 295
 296 int
 297 gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
 298                   gsl_vector_complex * alpha, gsl_vector * beta,
 299                   gsl_matrix * Q, gsl_matrix * Z,
 300                   gsl_eigen_gen_workspace * w)
 301 {
 302   if (Q && (A->size1 != Q->size1 || A->size1 != Q->size2))
 303     {
 304       GSL_ERROR("Q matrix has wrong dimensions", GSL_EBADLEN);
 305     }
 306   else if (Z && (A->size1 != Z->size1 || A->size1 != Z->size2))
 307     {
 308       GSL_ERROR("Z matrix has wrong dimensions", GSL_EBADLEN);
 309     }
 310   else
 311     {
 312       int s;
 313
 314       w->Q = Q;
 315       w->Z = Z;
 316
 317       s = gsl_eigen_gen(A, B, alpha, beta, w);
 318
 319       w->Q = NULL;
 320       w->Z = NULL;
 321
 322       return s;
 323     }
 324 } /* gsl_eigen_gen_QZ() */
 325
 326 /********************************************
 327  *           INTERNAL ROUTINES              *
 328  ********************************************/
 329
 330 /*
 331 gen_schur_decomp()
 332   Compute the generalized Schur decomposition of the matrix pencil
 333 (H, R) which is in Hessenberg-Triangular form
 334
 335 Inputs: H     - upper hessenberg matrix
 336         R     - upper triangular matrix
 337         alpha - (output) where to store eigenvalue numerators
 338         beta  - (output) where to store eigenvalue denominators
 339         w     - workspace
 340
 341 Return: none
 342
 343 Notes: 1) w->n_evals is updated to keep track of how many eigenvalues
 344           are found
 345 */
 346
 347 static void
 348 gen_schur_decomp(gsl_matrix *H, gsl_matrix *R, gsl_vector_complex *alpha,
 349                  gsl_vector *beta, gsl_eigen_gen_workspace *w)
 350 {
 351   size_t N;
 352   gsl_matrix_view h, r;
 353   gsl_matrix_view vh, vr;
 354   size_t q;             /* index of small subdiagonal element */
 355   gsl_complex z1, z2;   /* complex values */
 356   double a, b;
 357   int s;
 358   int flag;
 359
 360   N = H->size1;
 361
 362   h = gsl_matrix_submatrix(H, 0, 0, N, N);
 363   r = gsl_matrix_submatrix(R, 0, 0, N, N);
 364
 365   while ((N > 1) && (w->n_iter)++ < w->max_iterations)
 366     {
 367       q = gen_search_small_elements(&h.matrix, &r.matrix, &flag, w);
 368
 369       if (flag == 0)
 370         {
 371           /* no small elements found - do a QZ sweep */
 372           s = gen_qzstep(&h.matrix, &r.matrix, w);
 373
 374           if (s == GSL_CONTINUE)
 375             {
 376               /*
 377                * (h, r) is a 2-by-2 block with complex eigenvalues -
 378                * standardize and read off eigenvalues
 379                */
 380               s = gen_schur_standardize2(&h.matrix,
 381                                          &r.matrix,
 382                                          &z1,
 383                                          &z2,
 384                                          &a,
 385                                          &b,
 386                                          w);
 387               if (s != GSL_SUCCESS)
 388                 {
 389                   /*
 390                    * if we get here, then the standardization process
 391                    * perturbed the eigenvalues onto the real line -
 392                    * continue QZ iteration to break them into 1-by-1
 393                    * blocks
 394                    */
 395                   continue;
 396                 }
 397
 398               gen_store_eigval2(&h.matrix, &z1, a, &z2, b, alpha, beta, w);
 399               N = 0;
 400             } /* if (s) */
 401
 402           continue;
 403         } /* if (flag == 0) */
 404       else if (flag == 2)
 405         {
 406           if (q == 0)
 407             {
 408               /*
 409                * the leading element of R is zero, split off a block
 410                * at the top
 411                */
 412               gen_tri_split_top(&h.matrix, &r.matrix, w);
 413             }
 414           else
 415             {
 416               /*
 417                * we found a small element on the diagonal of R - chase the
 418                * zero to the bottom of the active block and then zero
 419                * H(n, n - 1) to split off a 1-by-1 block
 420                */
 421
 422               if (q != N - 1)
 423                 gen_tri_chase_zero(&h.matrix, &r.matrix, q, w);
 424
 425               /* now zero H(n, n - 1) */
 426               gen_tri_zero_H(&h.matrix, &r.matrix, w);
 427             }
 428
 429           /* continue so the next iteration detects the zero in H */
 430           continue;
 431         }
 432
 433       /*
 434        * a small subdiagonal element of H was found - one or two
 435        * eigenvalues have converged or the matrix has split into
 436        * two smaller matrices
 437        */
 438
 439       if (q == (N - 1))
 440         {
 441           /*
 442            * the last subdiagonal element of the hessenberg matrix is 0 -
 443            * H_{NN} / R_{NN} is a real eigenvalue - standardize so
 444            * R_{NN} > 0
 445            */
 446
 447           vh = gsl_matrix_submatrix(&h.matrix, q, q, 1, 1);
 448           vr = gsl_matrix_submatrix(&r.matrix, q, q, 1, 1);
 449           gen_schur_standardize1(&vh.matrix, &vr.matrix, &a, &b, w);
 450
 451           gen_store_eigval1(&vh.matrix, a, b, alpha, beta, w);
 452
 453           --N;
 454           h = gsl_matrix_submatrix(&h.matrix, 0, 0, N, N);
 455           r = gsl_matrix_submatrix(&r.matrix, 0, 0, N, N);
 456         }
 457       else if (q == (N - 2))
 458         {
 459           /* bottom right 2-by-2 block may have converged */
 460
 461           vh = gsl_matrix_submatrix(&h.matrix, q, q, 2, 2);
 462           vr = gsl_matrix_submatrix(&r.matrix, q, q, 2, 2);
 463           s = gen_schur_standardize2(&vh.matrix,
 464                                      &vr.matrix,
 465                                      &z1,
 466                                      &z2,
 467                                      &a,
 468                                      &b,
 469                                      w);
 470           if (s != GSL_SUCCESS)
 471             {
 472               /*
 473                * this 2-by-2 block contains real eigenvalues that
 474                * have not yet separated into 1-by-1 blocks -
 475                * recursively call gen_schur_decomp() to finish off
 476                * this block
 477                */
 478               gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);
 479             }
 480           else
 481             {
 482               /* we got 2 complex eigenvalues */
 483
 484               gen_store_eigval2(&vh.matrix, &z1, a, &z2, b, alpha, beta, w);
 485             }
 486
 487           N -= 2;
 488           h = gsl_matrix_submatrix(&h.matrix, 0, 0, N, N);
 489           r = gsl_matrix_submatrix(&r.matrix, 0, 0, N, N);
 490         }
 491       else if (q == 1)
 492         {
 493           /* H_{11} / R_{11} is an eigenvalue */
 494
 495           vh = gsl_matrix_submatrix(&h.matrix, 0, 0, 1, 1);
 496           vr = gsl_matrix_submatrix(&r.matrix, 0, 0, 1, 1);
 497           gen_schur_standardize1(&vh.matrix, &vr.matrix, &a, &b, w);
 498
 499           gen_store_eigval1(&vh.matrix, a, b, alpha, beta, w);
 500
 501           --N;
 502           h = gsl_matrix_submatrix(&h.matrix, 1, 1, N, N);
 503           r = gsl_matrix_submatrix(&r.matrix, 1, 1, N, N);
 504         }
 505       else if (q == 2)
 506         {
 507           /* upper left 2-by-2 block may have converged */
 508
 509           vh = gsl_matrix_submatrix(&h.matrix, 0, 0, 2, 2);
 510           vr = gsl_matrix_submatrix(&r.matrix, 0, 0, 2, 2);
 511           s = gen_schur_standardize2(&vh.matrix,
 512                                      &vr.matrix,
 513                                      &z1,
 514                                      &z2,
 515                                      &a,
 516                                      &b,
 517                                      w);
 518           if (s != GSL_SUCCESS)
 519             {
 520               /*
 521                * this 2-by-2 block contains real eigenvalues that
 522                * have not yet separated into 1-by-1 blocks -
 523                * recursively call gen_schur_decomp() to finish off
 524                * this block
 525                */
 526               gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);
 527             }
 528           else
 529             {
 530               /* we got 2 complex eigenvalues */
 531               gen_store_eigval2(&vh.matrix, &z1, a, &z2, b, alpha, beta, w);
 532             }
 533
 534           N -= 2;
 535           h = gsl_matrix_submatrix(&h.matrix, 2, 2, N, N);
 536           r = gsl_matrix_submatrix(&r.matrix, 2, 2, N, N);
 537         }
 538       else
 539         {
 540           /*
 541            * There is a zero element on the subdiagonal somewhere
 542            * in the middle of the matrix - we can now operate
 543            * separately on the two submatrices split by this
 544            * element. q is the row index of the zero element.
 545            */
 546
 547           /* operate on lower right (N - q)-by-(N - q) block first */
 548           vh = gsl_matrix_submatrix(&h.matrix, q, q, N - q, N - q);
 549           vr = gsl_matrix_submatrix(&r.matrix, q, q, N - q, N - q);
 550           gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);
 551
 552           /* operate on upper left q-by-q block */
 553           vh = gsl_matrix_submatrix(&h.matrix, 0, 0, q, q);
 554           vr = gsl_matrix_submatrix(&r.matrix, 0, 0, q, q);
 555           gen_schur_decomp(&vh.matrix, &vr.matrix, alpha, beta, w);
 556
 557           N = 0;
 558         }
 559     } /* while ((N > 1) && (w->n_iter)++ < w->max_iterations) */
 560
 561   /* handle special case of N = 1 */
 562
 563   if (N == 1)
 564     {
 565       gen_schur_standardize1(&h.matrix, &r.matrix, &a, &b, w);
 566       gen_store_eigval1(&h.matrix, a, b, alpha, beta, w);
 567     }
 568 } /* gen_schur_decomp() */
 569
 570 /*
 571 gen_qzstep()
 572   This routine determines what type of QZ step to perform on
 573 the generalized matrix pair (H, R). If the pair is 3-by-3 or bigger,
 574 we look at the bottom right 2-by-2 block. If this block has complex
 575 eigenvalues, we perform a Francis double shift QZ sweep. If it
 576 has real eigenvalues, we perform an implicit single shift QZ sweep.
 577
 578 If the pair is 2-by-2 with real eigenvalues, we perform a single
 579 shift sweep. If it has complex eigenvalues, we return GSL_CONTINUE
 580 to notify the calling function that a 2-by-2 block with complex
 581 eigenvalues has converged, so that it may then call
 582 gen_schur_standardize2(). In the real eigenvalue case, we want to
 583 continue doing QZ sweeps to break it up into two 1-by-1 blocks.
 584
 585 See LAPACK routine DHGEQZ and [1] for more information.
 586
 587 Inputs: H - upper Hessenberg matrix (at least 2-by-2)
 588         R - upper triangular matrix (at least 2-by-2)
 589         w - workspace
 590
 591 Return: GSL_SUCCESS on normal completion
 592         GSL_CONTINUE if we detect a 2-by-2 block with complex eigenvalues
 593 */
 594
 595 static inline int
 596 gen_qzstep(gsl_matrix *H, gsl_matrix *R, gsl_eigen_gen_workspace *w)
 597 {
 598   const size_t N = H->size1;
 599   gsl_matrix_view vh, vr; /* views of bottom right 2-by-2 block */
 600   double wr1, wr2, wi;
 601   double scale1, scale2, scale;
 602   double cs, sn;          /* givens rotation */
 603   double temp,            /* temporary variables */
 604          temp2;
 605   size_t j;               /* looping */
 606   gsl_vector_view xv, yv; /* temporary views */
 607   size_t top;
 608   size_t rows;
 609
 610   if (w->n_iter % 10 == 0)
 611     {
 612       /*
 613        * Exceptional shift - we have gone 10 iterations without finding
 614        * a new eigenvalue, do a single shift sweep with an
 615        * exceptional shift
 616        */
 617
 618       if ((GSL_DBL_MIN * w->max_iterations) *
 619           fabs(gsl_matrix_get(H, N - 2, N - 1)) <
 620           fabs(gsl_matrix_get(R, N - 2, N - 2)))
 621         {
 622           w->eshift += gsl_matrix_get(H, N - 2, N - 1) /
 623                        gsl_matrix_get(R, N - 2, N - 2);
 624         }
 625       else
 626         w->eshift += 1.0 / (GSL_DBL_MIN * w->max_iterations);
 627
 628       if ((w->eshift < GSL_DBL_EPSILON) &&
 629           (GSL_DBL_MIN * w->max_iterations) *
 630           fabs(gsl_matrix_get(H, N - 1, N - 2)) <
 631           fabs(gsl_matrix_get(R, N - 2, N - 2)))
 632         {
 633           w->eshift = GEN_ESHIFT_COEFF *
 634                       (w->ascale * gsl_matrix_get(H, N - 1, N - 2)) /
 635                       (w->bscale * gsl_matrix_get(R, N - 2, N - 2));
 636         }
 637
 638       scale1 = 1.0;
 639       wr1 = w->eshift;
 640     }
 641   else
 642     {
 643       /*
 644        * Compute generalized eigenvalues of bottom right 2-by-2 block
 645        * to be used as shifts - wr1 is the Wilkinson shift
 646        */
 647
 648       vh = gsl_matrix_submatrix(H, N - 2, N - 2, 2, 2);
 649       vr = gsl_matrix_submatrix(R, N - 2, N - 2, 2, 2);
 650       gsl_schur_gen_eigvals(&vh.matrix,
 651                             &vr.matrix,
 652                             &wr1,
 653                             &wr2,
 654                             &wi,
 655                             &scale1,
 656                             &scale2);
 657
 658       if (wi != 0.0)
 659         {
 660           /* complex eigenvalues */
 661
 662           if (N == 2)
 663             {
 664               /*
 665                * its a 2-by-2 block with complex eigenvalues - notify
 666                * the calling function to deflate
 667                */
 668               return (GSL_CONTINUE);
 669             }
 670           else
 671             {
 672               /* do a francis double shift sweep */
 673               gen_qzstep_d(H, R, w);
 674             }
 675
 676           return GSL_SUCCESS;
 677         }
 678     }
 679
 680   /* real eigenvalues - perform single shift QZ step */
 681
 682   temp = GSL_MIN(w->ascale, 1.0) * (0.5 / GSL_DBL_MIN);
 683   if (scale1 > temp)
 684     scale = temp / scale1;
 685   else
 686     scale = 1.0;
 687
 688   temp = GSL_MIN(w->bscale, 1.0) * (0.5 / GSL_DBL_MIN);
 689   if (fabs(wr1) > temp)
 690     scale = GSL_MIN(scale, temp / fabs(wr1));
 691
 692   scale1 *= scale;
 693   wr1 *= scale;
 694
 695   if (w->needtop)
 696     {
 697       /* get absolute index of this matrix relative to original matrix */
 698       top = gen_get_submatrix(w->H, H);
 699     }
 700
 701   temp = scale1*gsl_matrix_get(H, 0, 0) - wr1*gsl_matrix_get(R, 0, 0);
 702   temp2 = scale1*gsl_matrix_get(H, 1, 0);
 703
 704   create_givens(temp, temp2, &cs, &sn);
 705   sn = -sn;
 706
 707   for (j = 0; j < N - 1; ++j)
 708     {
 709       if (j > 0)
 710         {
 711           temp = gsl_matrix_get(H, j, j - 1);
 712           temp2 = gsl_matrix_get(H, j + 1, j - 1);
 713           create_givens(temp, temp2, &cs, &sn);
 714           sn = -sn;
 715
 716           /* apply to column (j - 1) */
 717           temp = cs * gsl_matrix_get(H, j, j - 1) +
 718                  sn * gsl_matrix_get(H, j + 1, j - 1);
 719           gsl_matrix_set(H, j, j - 1, temp);
 720           gsl_matrix_set(H, j + 1, j - 1, 0.0);
 721         }
 722
 723       /* apply G to H(j:j+1,:) and T(j:j+1,:) */
 724
 725       if (w->compute_s)
 726         {
 727           xv = gsl_matrix_subrow(w->H, top + j, top + j, w->size - top - j);
 728           yv = gsl_matrix_subrow(w->H, top + j + 1, top + j, w->size - top - j);
 729         }
 730       else
 731         {
 732           xv = gsl_matrix_subrow(H, j, j, N - j);
 733           yv = gsl_matrix_subrow(H, j + 1, j, N - j);
 734         }
 735
 736       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 737
 738       if (w->compute_t)
 739         {
 740           xv = gsl_matrix_subrow(w->R, top + j, top + j, w->size - top - j);
 741           yv = gsl_matrix_subrow(w->R, top + j + 1, top + j, w->size - top - j);
 742         }
 743       else
 744         {
 745           xv = gsl_matrix_subrow(R, j, j, N - j);
 746           yv = gsl_matrix_subrow(R, j + 1, j, N - j);
 747         }
 748
 749       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 750
 751       if (w->Q)
 752         {
 753           /* accumulate Q: Q -> QG */
 754           xv = gsl_matrix_column(w->Q, top + j);
 755           yv = gsl_matrix_column(w->Q, top + j + 1);
 756           gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 757         }
 758
 759       temp = gsl_matrix_get(R, j + 1, j + 1);
 760       temp2 = gsl_matrix_get(R, j + 1, j);
 761       create_givens(temp, temp2, &cs, &sn);
 762
 763       rows = GSL_MIN(j + 3, N);
 764
 765       if (w->compute_s)
 766         {
 767           xv = gsl_matrix_subcolumn(w->H, top + j, 0, top + rows);
 768           yv = gsl_matrix_subcolumn(w->H, top + j + 1, 0, top + rows);
 769         }
 770       else
 771         {
 772           xv = gsl_matrix_subcolumn(H, j, 0, rows);
 773           yv = gsl_matrix_subcolumn(H, j + 1, 0, rows);
 774         }
 775
 776       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 777
 778       rows = GSL_MIN(j + 2, N);
 779
 780       if (w->compute_t)
 781         {
 782           xv = gsl_matrix_subcolumn(w->R, top + j, 0, top + rows);
 783           yv = gsl_matrix_subcolumn(w->R, top + j + 1, 0, top + rows);
 784         }
 785       else
 786         {
 787           xv = gsl_matrix_subcolumn(R, j, 0, rows);
 788           yv = gsl_matrix_subcolumn(R, j + 1, 0, rows);
 789         }
 790
 791       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 792
 793       if (w->Z)
 794         {
 795           /* accumulate Z: Z -> ZG */
 796           xv = gsl_matrix_column(w->Z, top + j);
 797           yv = gsl_matrix_column(w->Z, top + j + 1);
 798           gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 799         }
 800     } /* for (j = 0; j < N - 1; ++j) */
 801
 802   return GSL_SUCCESS;
 803 } /* gen_qzstep() */
 804
 805 /*
 806 gen_qzstep_d()
 807   Perform an implicit double shift QZ step.
 808
 809 See Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm 7.7.2
 810
 811 Inputs: H - upper Hessenberg matrix (at least 3-by-3)
 812         R - upper triangular matrix (at least 3-by-3)
 813         w - workspace
 814 */
 815
 816 static inline void
 817 gen_qzstep_d(gsl_matrix *H, gsl_matrix *R, gsl_eigen_gen_workspace *w)
 818 {
 819   const size_t N = H->size1;
 820   size_t j;               /* looping */
 821   double dat[3];          /* householder vector */
 822   double tau;             /* householder coefficient */
 823   gsl_vector_view v2, v3; /* views into 'dat' */
 824   gsl_matrix_view m;      /* temporary view */
 825   double tmp;
 826   size_t q, r;
 827   size_t top;             /* location of H in original matrix */
 828   double scale;
 829   double AB11,            /* various matrix element ratios */
 830          AB22,
 831          ABNN,
 832          ABMM,
 833          AMNBNN,
 834          ANMBMM,
 835          A21B11,
 836          A12B22,
 837          A32B22,
 838          B12B22,
 839          BMNBNN;
 840
 841   v2 = gsl_vector_view_array(dat, 2);
 842   v3 = gsl_vector_view_array(dat, 3);
 843
 844   if (w->needtop)
 845     {
 846       /* get absolute index of this matrix relative to original matrix */
 847       top = gen_get_submatrix(w->H, H);
 848     }
 849
 850   /*
 851    * Similar to the QR method, we take the shifts to be the two
 852    * zeros of the problem
 853    *
 854    * det[H(n-1:n,n-1:n) - s*R(n-1:n,n-1:n)] = 0
 855    *
 856    * The initial householder vector elements are then given by
 857    * Eq. 4.1 of [1], which are designed to reduce errors when
 858    * off diagonal elements are small.
 859    */
 860
 861   ABMM = (w->ascale * gsl_matrix_get(H, N - 2, N - 2)) /
 862          (w->bscale * gsl_matrix_get(R, N - 2, N - 2));
 863   ABNN = (w->ascale * gsl_matrix_get(H, N - 1, N - 1)) /
 864          (w->bscale * gsl_matrix_get(R, N - 1, N - 1));
 865   AB11 = (w->ascale * gsl_matrix_get(H, 0, 0)) /
 866          (w->bscale * gsl_matrix_get(R, 0, 0));
 867   AB22 = (w->ascale * gsl_matrix_get(H, 1, 1)) /
 868          (w->bscale * gsl_matrix_get(R, 1, 1));
 869   AMNBNN = (w->ascale * gsl_matrix_get(H, N - 2, N - 1)) /
 870            (w->bscale * gsl_matrix_get(R, N - 1, N - 1));
 871   ANMBMM = (w->ascale * gsl_matrix_get(H, N - 1, N - 2)) /
 872            (w->bscale * gsl_matrix_get(R, N - 2, N - 2));
 873   BMNBNN = gsl_matrix_get(R, N - 2, N - 1) /
 874            gsl_matrix_get(R, N - 1, N - 1);
 875   A21B11 = (w->ascale * gsl_matrix_get(H, 1, 0)) /
 876            (w->bscale * gsl_matrix_get(R, 0, 0));
 877   A12B22 = (w->ascale * gsl_matrix_get(H, 0, 1)) /
 878            (w->bscale * gsl_matrix_get(R, 1, 1));
 879   A32B22 = (w->ascale * gsl_matrix_get(H, 2, 1)) /
 880            (w->bscale * gsl_matrix_get(R, 1, 1));
 881   B12B22 = gsl_matrix_get(R, 0, 1) / gsl_matrix_get(R, 1, 1);
 882
 883   /*
 884    * These are the Eqs (4.1) of [1], just multiplied by the factor
 885    * (A_{21} / B_{11})
 886    */
 887   dat[0] = (ABMM - AB11) * (ABNN - AB11) - (AMNBNN * ANMBMM) +
 888            (ANMBMM * BMNBNN * AB11) + (A12B22 - (AB11 * B12B22)) * A21B11;
 889   dat[1] = ((AB22 - AB11) - (A21B11 * B12B22) - (ABMM - AB11) -
 890             (ABNN - AB11) + (ANMBMM * BMNBNN)) * A21B11;
 891   dat[2] = A32B22 * A21B11;
 892
 893   scale = fabs(dat[0]) + fabs(dat[1]) + fabs(dat[2]);
 894   if (scale != 0.0)
 895     {
 896       dat[0] /= scale;
 897       dat[1] /= scale;
 898       dat[2] /= scale;
 899     }
 900
 901   for (j = 0; j < N - 2; ++j)
 902     {
 903       r = GSL_MIN(j + 4, N);
 904
 905       /*
 906        * Find householder Q so that
 907        *
 908        * Q [x y z]^t = [ * 0 0 ]^t
 909        */
 910
 911       tau = gsl_linalg_householder_transform(&v3.vector);
 912
 913       if (tau != 0.0)
 914         {
 915           /*
 916            * q is the initial column to start applying the householder
 917            * transformation. The GSL_MAX() simply ensures we don't
 918            * try to apply it to column (-1), since we are zeroing out
 919            * column (j - 1) except for the first iteration which
 920            * introduces the bulge.
 921            */
 922           q = (size_t) GSL_MAX(0, (int)j - 1);
 923
 924           /* H -> QH, R -> QR */
 925
 926           if (w->compute_s)
 927             {
 928               /*
 929                * We are computing the Schur form S, so we need to
 930                * transform the whole matrix H
 931                */
 932               m = gsl_matrix_submatrix(w->H,
 933                                        top + j,
 934                                        top + q,
 935                                        3,
 936                                        w->size - top - q);
 937               gsl_linalg_householder_hm(tau, &v3.vector, &m.matrix);
 938             }
 939           else
 940             {
 941               /* just transform the active block */
 942               m = gsl_matrix_submatrix(H, j, q, 3, N - q);
 943               gsl_linalg_householder_hm(tau, &v3.vector, &m.matrix);
 944             }
 945
 946           if (w->compute_t)
 947             {
 948               /*
 949                * We are computing the Schur form T, so we need to
 950                * transform the whole matrix R
 951                */
 952               m = gsl_matrix_submatrix(w->R,
 953                                        top + j,
 954                                        top + j,
 955                                        3,
 956                                        w->size - top - j);
 957               gsl_linalg_householder_hm(tau, &v3.vector, &m.matrix);
 958             }
 959           else
 960             {
 961               /* just transform the active block */
 962               m = gsl_matrix_submatrix(R, j, j, 3, N - j);
 963               gsl_linalg_householder_hm(tau, &v3.vector, &m.matrix);
 964             }
 965
 966           if (w->Q)
 967             {
 968               /* accumulate the transformation into Q */
 969               m = gsl_matrix_submatrix(w->Q, 0, top + j, w->size, 3);
 970               gsl_linalg_householder_mh(tau, &v3.vector, &m.matrix);
 971             }
 972         } /* if (tau != 0.0) */
 973
 974       /*
 975        * Find householder Z so that
 976        *
 977        * [ r_{j+2,j} r_{j+2, j+1}, r_{j+2, j+2} ] Z = [ 0 0 * ]
 978        *
 979        * This isn't exactly what gsl_linalg_householder_transform
 980        * does, so we need to rotate the input vector so it preserves
 981        * the last element, and then rotate it back afterwards.
 982        *
 983        * So instead of transforming [x y z], we transform [z x y],
 984        * and the resulting HH vector [1 v2 v3] -> [v2 v3 1] but
 985        * then needs to be scaled to have the first element = 1, so
 986        * it becomes [1 v3/v2 1/v2] (tau must also be scaled accordingly).
 987        */
 988
 989       dat[0] = gsl_matrix_get(R, j + 2, j + 2);
 990       dat[1] = gsl_matrix_get(R, j + 2, j);
 991       dat[2] = gsl_matrix_get(R, j + 2, j + 1);
 992       scale = fabs(dat[0]) + fabs(dat[1]) + fabs(dat[2]);
 993       if (scale != 0.0)
 994         {
 995           dat[0] /= scale;
 996           dat[1] /= scale;
 997           dat[2] /= scale;
 998         }
 999
1000       tau = gsl_linalg_householder_transform(&v3.vector);
1001
1002       if (tau != 0.0)
1003         {
1004           /* rotate back */
1005           tmp = gsl_vector_get(&v3.vector, 1);
1006           gsl_vector_set(&v3.vector, 1, gsl_vector_get(&v3.vector, 2)/tmp);
1007           gsl_vector_set(&v3.vector, 2, 1.0 / tmp);
1008           tau *= tmp * tmp;
1009
1010           /* H -> HZ, R -> RZ */
1011
1012           if (w->compute_s)
1013             {
1014               m = gsl_matrix_submatrix(w->H, 0, top + j, top + r, 3);
1015               gsl_linalg_householder_mh(tau, &v3.vector, &m.matrix);
1016             }
1017           else
1018             {
1019               m = gsl_matrix_submatrix(H, 0, j, r, 3);
1020               gsl_linalg_householder_mh(tau, &v3.vector, &m.matrix);
1021             }
1022
1023           if (w->compute_t)
1024             {
1025               m = gsl_matrix_submatrix(w->R, 0, top + j, top + j + 3, 3);
1026               gsl_linalg_householder_mh(tau, &v3.vector, &m.matrix);
1027             }
1028           else
1029             {
1030               m = gsl_matrix_submatrix(R, 0, j, j + 3, 3);
1031               gsl_linalg_householder_mh(tau, &v3.vector, &m.matrix);
1032             }
1033
1034           if (w->Z)
1035             {
1036               /* accumulate transformation into Z */
1037               m = gsl_matrix_submatrix(w->Z, 0, top + j, w->size, 3);
1038               gsl_linalg_householder_mh(tau, &v3.vector, &m.matrix);
1039             }
1040         } /* if (tau != 0.0) */
1041
1042       /*
1043        * Find householder Z so that
1044        *
1045        * [ r_{j+1,j} r_{j+1, j+1} ] Z = [ 0 * ]
1046        */
1047
1048       dat[0] = gsl_matrix_get(R, j + 1, j + 1);
1049       dat[1] = gsl_matrix_get(R, j + 1, j);
1050       scale = fabs(dat[0]) + fabs(dat[1]);
1051       if (scale != 0.0)
1052         {
1053           dat[0] /= scale;
1054           dat[1] /= scale;
1055         }
1056
1057       tau = gsl_linalg_householder_transform(&v2.vector);
1058
1059       if (tau != 0.0)
1060         {
1061           /* rotate back */
1062           tmp = gsl_vector_get(&v2.vector, 1);
1063           gsl_vector_set(&v2.vector, 1, 1.0 / tmp);
1064           tau *= tmp * tmp;
1065
1066           /* H -> HZ, R -> RZ */
1067
1068           if (w->compute_s)
1069             {
1070               m = gsl_matrix_submatrix(w->H, 0, top + j, top + r, 2);
1071               gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1072             }
1073           else
1074             {
1075               m = gsl_matrix_submatrix(H, 0, j, r, 2);
1076               gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1077             }
1078
1079           if (w->compute_t)
1080             {
1081               m = gsl_matrix_submatrix(w->R, 0, top + j, top + j + 3, 2);
1082               gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1083             }
1084           else
1085             {
1086               m = gsl_matrix_submatrix(R, 0, j, j + 3, 2);
1087               gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1088             }
1089
1090           if (w->Z)
1091             {
1092               /* accumulate transformation into Z */
1093               m = gsl_matrix_submatrix(w->Z, 0, top + j, w->size, 2);
1094               gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1095             }
1096         } /* if (tau != 0.0) */
1097
1098       dat[0] = gsl_matrix_get(H, j + 1, j);
1099       dat[1] = gsl_matrix_get(H, j + 2, j);
1100       if (j < N - 3)
1101         dat[2] = gsl_matrix_get(H, j + 3, j);
1102
1103       scale = fabs(dat[0]) + fabs(dat[1]) + fabs(dat[2]);
1104       if (scale != 0.0)
1105         {
1106           dat[0] /= scale;
1107           dat[1] /= scale;
1108           dat[2] /= scale;
1109         }
1110     } /* for (j = 0; j < N - 2; ++j) */
1111
1112   /*
1113    * Find Householder Q so that
1114    *
1115    * Q [ x y ]^t = [ * 0 ]^t
1116    */
1117
1118   scale = fabs(dat[0]) + fabs(dat[1]);
1119   if (scale != 0.0)
1120     {
1121       dat[0] /= scale;
1122       dat[1] /= scale;
1123     }
1124
1125   tau = gsl_linalg_householder_transform(&v2.vector);
1126
1127   if (w->compute_s)
1128     {
1129       m = gsl_matrix_submatrix(w->H,
1130                                top + N - 2,
1131                                top + N - 3,
1132                                2,
1133                                w->size - top - N + 3);
1134       gsl_linalg_householder_hm(tau, &v2.vector, &m.matrix);
1135     }
1136   else
1137     {
1138       m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
1139       gsl_linalg_householder_hm(tau, &v2.vector, &m.matrix);
1140     }
1141
1142   if (w->compute_t)
1143     {
1144       m = gsl_matrix_submatrix(w->R,
1145                                top + N - 2,
1146                                top + N - 2,
1147                                2,
1148                                w->size - top - N + 2);
1149       gsl_linalg_householder_hm(tau, &v2.vector, &m.matrix);
1150     }
1151   else
1152     {
1153       m = gsl_matrix_submatrix(R, N - 2, N - 2, 2, 2);
1154       gsl_linalg_householder_hm(tau, &v2.vector, &m.matrix);
1155     }
1156
1157   if (w->Q)
1158     {
1159       /* accumulate the transformation into Q */
1160       m = gsl_matrix_submatrix(w->Q, 0, top + N - 2, w->size, 2);
1161       gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1162     }
1163
1164   /*
1165    * Find Householder Z so that
1166    *
1167    * [ b_{n,n-1} b_{nn} ] Z = [ 0 * ]
1168    */
1169
1170   dat[0] = gsl_matrix_get(R, N - 1, N - 1);
1171   dat[1] = gsl_matrix_get(R, N - 1, N - 2);
1172   scale = fabs(dat[0]) + fabs(dat[1]);
1173   if (scale != 0.0)
1174     {
1175       dat[0] /= scale;
1176       dat[1] /= scale;
1177     }
1178
1179   tau = gsl_linalg_householder_transform(&v2.vector);
1180
1181   /* rotate back */
1182   tmp = gsl_vector_get(&v2.vector, 1);
1183   gsl_vector_set(&v2.vector, 1, 1.0 / tmp);
1184   tau *= tmp * tmp;
1185
1186   if (w->compute_s)
1187     {
1188       m = gsl_matrix_submatrix(w->H, 0, top + N - 2, top + N, 2);
1189       gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1190     }
1191   else
1192     {
1193       m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
1194       gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1195     }
1196
1197   if (w->compute_t)
1198     {
1199       m = gsl_matrix_submatrix(w->R, 0, top + N - 2, top + N, 2);
1200       gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1201     }
1202   else
1203     {
1204       m = gsl_matrix_submatrix(R, 0, N - 2, N, 2);
1205       gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1206     }
1207
1208   if (w->Z)
1209     {
1210       /* accumulate the transformation into Z */
1211       m = gsl_matrix_submatrix(w->Z, 0, top + N - 2, w->size, 2);
1212       gsl_linalg_householder_mh(tau, &v2.vector, &m.matrix);
1213     }
1214 } /* gen_qzstep_d() */
1215
1216 /*
1217 gen_tri_split_top()
1218   This routine is called when the leading element on the diagonal of R
1219 has become negligible. Split off a 1-by-1 block at the top.
1220
1221 Inputs: H - upper hessenberg matrix
1222         R - upper triangular matrix
1223         w - workspace
1224 */
1225
1226 static void
1227 gen_tri_split_top(gsl_matrix *H, gsl_matrix *R, gsl_eigen_gen_workspace *w)
1228 {
1229   const size_t N = H->size1;
1230   size_t j, top;
1231   double cs, sn;
1232   gsl_vector_view xv, yv;
1233
1234   if (w->needtop)
1235     top = gen_get_submatrix(w->H, H);
1236
1237   j = 0;
1238
1239   create_givens(gsl_matrix_get(H, j, j),
1240                 gsl_matrix_get(H, j + 1, j),
1241                 &cs,
1242                 &sn);
1243   sn = -sn;
1244
1245   if (w->compute_s)
1246     {
1247       xv = gsl_matrix_subrow(w->H, top + j, top, w->size - top);
1248       yv = gsl_matrix_subrow(w->H, top + j + 1, top, w->size - top);
1249     }
1250   else
1251     {
1252       xv = gsl_matrix_row(H, j);
1253       yv = gsl_matrix_row(H, j + 1);
1254     }
1255
1256   gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1257   gsl_matrix_set(H, j + 1, j, 0.0);
1258
1259   if (w->compute_t)
1260     {
1261       xv = gsl_matrix_subrow(w->R, top + j, top + 1, w->size - top - 1);
1262       yv = gsl_matrix_subrow(w->R, top + j + 1, top + 1, w->size - top - 1);
1263     }
1264   else
1265     {
1266       xv = gsl_matrix_subrow(R, j, 1, N - 1);
1267       yv = gsl_matrix_subrow(R, j + 1, 1, N - 1);
1268     }
1269
1270   gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1271
1272   if (w->Q)
1273     {
1274       xv = gsl_matrix_column(w->Q, top + j);
1275       yv = gsl_matrix_column(w->Q, top + j + 1);
1276       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1277     }
1278 } /* gen_tri_split_top() */
1279
1280 /*
1281 gen_tri_chase_zero()
1282   This routine is called when an element on the diagonal of R
1283 has become negligible. Chase the zero to the bottom of the active
1284 block so we can split off a 1-by-1 block.
1285
1286 Inputs: H - upper hessenberg matrix
1287         R - upper triangular matrix
1288         q - index such that R(q,q) = 0 (q must be > 0)
1289         w - workspace
1290 */
1291
1292 static inline void
1293 gen_tri_chase_zero(gsl_matrix *H, gsl_matrix *R, size_t q,
1294                    gsl_eigen_gen_workspace *w)
1295 {
1296   const size_t N = H->size1;
1297   size_t j, top;
1298   double cs, sn;
1299   gsl_vector_view xv, yv;
1300
1301   if (w->needtop)
1302     top = gen_get_submatrix(w->H, H);
1303
1304   for (j = q; j < N - 1; ++j)
1305     {
1306       create_givens(gsl_matrix_get(R, j, j + 1),
1307                     gsl_matrix_get(R, j + 1, j + 1),
1308                     &cs,
1309                     &sn);
1310       sn = -sn;
1311
1312       if (w->compute_t)
1313         {
1314           xv = gsl_matrix_subrow(w->R, top + j, top + j + 1, w->size - top - j - 1);
1315           yv = gsl_matrix_subrow(w->R, top + j + 1, top + j + 1, w->size - top - j - 1);
1316         }
1317       else
1318         {
1319           xv = gsl_matrix_subrow(R, j, j + 1, N - j - 1);
1320           yv = gsl_matrix_subrow(R, j + 1, j + 1, N - j - 1);
1321         }
1322
1323       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1324       gsl_matrix_set(R, j + 1, j + 1, 0.0);
1325
1326       if (w->compute_s)
1327         {
1328           xv = gsl_matrix_subrow(w->H, top + j, top + j - 1, w->size - top - j + 1);
1329           yv = gsl_matrix_subrow(w->H, top + j + 1, top + j - 1, w->size - top - j + 1);
1330         }
1331       else
1332         {
1333           xv = gsl_matrix_subrow(H, j, j - 1, N - j + 1);
1334           yv = gsl_matrix_subrow(H, j + 1, j - 1, N - j + 1);
1335         }
1336
1337       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1338
1339       if (w->Q)
1340         {
1341           /* accumulate Q */
1342           xv = gsl_matrix_column(w->Q, top + j);
1343           yv = gsl_matrix_column(w->Q, top + j + 1);
1344           gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1345         }
1346
1347       create_givens(gsl_matrix_get(H, j + 1, j),
1348                     gsl_matrix_get(H, j + 1, j - 1),
1349                     &cs,
1350                     &sn);
1351       sn = -sn;
1352
1353       if (w->compute_s)
1354         {
1355           xv = gsl_matrix_subcolumn(w->H, top + j, 0, top + j + 2);
1356           yv = gsl_matrix_subcolumn(w->H, top + j - 1, 0, top + j + 2);
1357         }
1358       else
1359         {
1360           xv = gsl_matrix_subcolumn(H, j, 0, j + 2);
1361           yv = gsl_matrix_subcolumn(H, j - 1, 0, j + 2);
1362         }
1363
1364       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1365       gsl_matrix_set(H, j + 1, j - 1, 0.0);
1366
1367       if (w->compute_t)
1368         {
1369           xv = gsl_matrix_subcolumn(w->R, top + j, 0, top + j + 1);
1370           yv = gsl_matrix_subcolumn(w->R, top + j - 1, 0, top + j + 1);
1371         }
1372       else
1373         {
1374           xv = gsl_matrix_subcolumn(R, j, 0, j + 1);
1375           yv = gsl_matrix_subcolumn(R, j - 1, 0, j + 1);
1376         }
1377
1378       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1379
1380       if (w->Z)
1381         {
1382           /* accumulate Z */
1383           xv = gsl_matrix_column(w->Z, top + j);
1384           yv = gsl_matrix_column(w->Z, top + j - 1);
1385           gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1386         }
1387     }
1388 } /* gen_tri_chase_zero() */
1389
1390 /*
1391 gen_tri_zero_H()
1392   Companion function to get_tri_chase_zero(). After the zero on
1393 the diagonal of R has been chased to the bottom, we zero the element
1394 H(n, n - 1) in order to split off a 1-by-1 block.
1395 */
1396
1397 static inline void
1398 gen_tri_zero_H(gsl_matrix *H, gsl_matrix *R, gsl_eigen_gen_workspace *w)
1399 {
1400   const size_t N = H->size1;
1401   size_t top;
1402   double cs, sn;
1403   gsl_vector_view xv, yv;
1404
1405   if (w->needtop)
1406     top = gen_get_submatrix(w->H, H);
1407
1408   create_givens(gsl_matrix_get(H, N - 1, N - 1),
1409                 gsl_matrix_get(H, N - 1, N - 2),
1410                 &cs,
1411                 &sn);
1412   sn = -sn;
1413
1414   if (w->compute_s)
1415     {
1416       xv = gsl_matrix_subcolumn(w->H, top + N - 1, 0, top + N);
1417       yv = gsl_matrix_subcolumn(w->H, top + N - 2, 0, top + N);
1418     }
1419   else
1420     {
1421       xv = gsl_matrix_column(H, N - 1);
1422       yv = gsl_matrix_column(H, N - 2);
1423     }
1424
1425   gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1426
1427   gsl_matrix_set(H, N - 1, N - 2, 0.0);
1428
1429   if (w->compute_t)
1430     {
1431       xv = gsl_matrix_subcolumn(w->R, top + N - 1, 0, top + N - 1);
1432       yv = gsl_matrix_subcolumn(w->R, top + N - 2, 0, top + N - 1);
1433     }
1434   else
1435     {
1436       xv = gsl_matrix_subcolumn(R, N - 1, 0, N - 1);
1437       yv = gsl_matrix_subcolumn(R, N - 2, 0, N - 1);
1438     }
1439
1440   gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1441
1442   if (w->Z)
1443     {
1444       /* accumulate Z */
1445       xv = gsl_matrix_column(w->Z, top + N - 1);
1446       yv = gsl_matrix_column(w->Z, top + N - 2);
1447       gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
1448     }
1449 } /* gen_tri_zero_H() */
1450
1451 /*
1452 gen_search_small_elements()
1453   This routine searches for small elements in the matrix pencil
1454 (H, R) to determine if any eigenvalues have converged.
1455
1456 Tests:
1457
1458 1. Test if the Hessenberg matrix has a small subdiagonal element:
1459    H(i, i - 1) < tolerance
1460
1461 2. Test if the Triangular matrix has a small diagonal element:
1462    R(i, i) < tolerance
1463
1464 Possible outcomes:
1465
1466 (A) Neither test passed: in this case 'flag' is set to 0 and the
1467     function returns 0
1468
1469 (B) Test 1 passes and 2 does not: in this case 'flag' is set to 1
1470     and we return the row index i such that H(i, i - 1) < tol
1471
1472 (C) Test 2 passes and 1 does not: in this case 'flag' is set to 2
1473     and we return the index i such that R(i, i) < tol
1474
1475 (D) Tests 1 and 2 both pass: in this case 'flag' is set to 3 and
1476     we return the index i such that H(i, i - 1) < tol and R(i, i) < tol
1477
1478 Inputs: H    - upper Hessenberg matrix
1479         R    - upper Triangular matrix
1480         flag - (output) flag set on output (see above)
1481         w    - workspace
1482
1483 Return: see above
1484 */
1485
1486 static inline size_t
1487 gen_search_small_elements(gsl_matrix *H, gsl_matrix *R,
1488                           int *flag, gsl_eigen_gen_workspace *w)
1489 {
1490   const size_t N = H->size1;
1491   int k;
1492   size_t i;
1493   int pass1 = 0;
1494   int pass2 = 0;
1495
1496   for (k = (int) N - 1; k >= 0; --k)
1497     {
1498       i = (size_t) k;
1499
1500       if (i != 0 && fabs(gsl_matrix_get(H, i, i - 1)) <= w->atol)
1501         {
1502           gsl_matrix_set(H, i, i - 1, 0.0);
1503           pass1 = 1;
1504         }
1505
1506       if (fabs(gsl_matrix_get(R, i, i)) < w->btol)
1507         {
1508           gsl_matrix_set(R, i, i, 0.0);
1509           pass2 = 1;
1510         }
1511
1512       if (pass1 && !pass2)      /* case B */
1513         {
1514           *flag = 1;
1515           return (i);
1516         }
1517       else if (!pass1 && pass2) /* case C */
1518         {
1519           *flag = 2;
1520           return (i);
1521         }
1522       else if (pass1 && pass2)  /* case D */
1523         {
1524           *flag = 3;
1525           return (i);
1526         }
1527     }
1528
1529   /* neither test passed: case A */
1530
1531   *flag = 0;
1532   return (0);
1533 } /* gen_search_subdiag_small_elements() */
1534
1535 /*
1536 gen_schur_standardize1()
1537   This function is called when a 1-by-1 block has converged -
1538 convert the block to standard form and update the Schur forms and
1539 vectors if required. Standard form here means that the diagonal
1540 element of B is positive.
1541
1542 Inputs: A      - 1-by-1 matrix in Schur form S
1543         B      - 1-by-1 matrix in Schur form T
1544         alphar - where to store real part of eigenvalue numerator
1545         beta   - where to store eigenvalue denominator
1546         w      - workspace
1547
1548 Return: success
1549 */
1550
1551 static int
1552 gen_schur_standardize1(gsl_matrix *A, gsl_matrix *B, double *alphar,
1553                        double *beta, gsl_eigen_gen_workspace *w)
1554 {
1555   size_t i;
1556   size_t top;
1557
1558   /*
1559    * it is a 1-by-1 block - the only requirement is that
1560    * B_{00} is > 0, so if it isn't apply a -I transformation
1561    */
1562   if (gsl_matrix_get(B, 0, 0) < 0.0)
1563     {
1564       if (w->needtop)
1565         top = gen_get_submatrix(w->H, A);
1566
1567       if (w->compute_t)
1568         {
1569           for (i = 0; i <= top; ++i)
1570             gsl_matrix_set(w->R, i, top, -gsl_matrix_get(w->R, i, top));
1571         }
1572       else
1573         gsl_matrix_set(B, 0, 0, -gsl_matrix_get(B, 0, 0));
1574
1575       if (w->compute_s)
1576         {
1577           for (i = 0; i <= top; ++i)
1578             gsl_matrix_set(w->H, i, top, -gsl_matrix_get(w->H, i, top));
1579         }
1580       else
1581         gsl_matrix_set(A, 0, 0, -gsl_matrix_get(A, 0, 0));
1582
1583       if (w->Z)
1584         {
1585           for (i = 0; i < w->size; ++i)
1586             gsl_matrix_set(w->Z, i, top, -gsl_matrix_get(w->Z, i, top));
1587         }
1588     }
1589
1590   *alphar = gsl_matrix_get(A, 0, 0);
1591   *beta = gsl_matrix_get(B, 0, 0);
1592
1593   return GSL_SUCCESS;
1594 } /* gen_schur_standardize1() */
1595
1596 /*
1597 gen_schur_standardize2()
1598   This function is called when a 2-by-2 generalized block has
1599 converged. Convert the block to standard form, which means B
1600 is rotated so that
1601
1602 B = [ B11  0  ] with B11, B22 non-negative
1603     [  0  B22 ]
1604
1605 If the resulting block (A, B) has complex eigenvalues, they are
1606 computed. Otherwise, the function will return GSL_CONTINUE to
1607 notify caller that we need to do more single shift sweeps to
1608 convert the 2-by-2 block into two 1-by-1 blocks.
1609
1610 Inputs: A      - 2-by-2 submatrix of schur form S
1611         B      - 2-by-2 submatrix of schur form T
1612         alpha1 - (output) where to store eigenvalue 1 numerator
1613         alpha2 - (output) where to store eigenvalue 2 numerator
1614         beta1  - (output) where to store eigenvalue 1 denominator
1615         beta2  - (output) where to store eigenvalue 2 denominator
1616         w      - workspace
1617
1618 Return: GSL_SUCCESS if block has complex eigenvalues (they are computed)
1619         GSL_CONTINUE if block has real eigenvalues (they are not computed)
1620 */
1621
1622 static int
1623 gen_schur_standardize2(gsl_matrix *A, gsl_matrix *B, gsl_complex *alpha1,
1624                        gsl_complex *alpha2, double *beta1, double *beta2,
1625                        gsl_eigen_gen_workspace *w)
1626 {
1627   double datB[4],
1628          datV[4],
1629          datS[2],
1630          work[2];
1631   gsl_matrix_view uv = gsl_matrix_view_array(datB, 2, 2);
1632   gsl_matrix_view vv = gsl_matrix_view_array(datV, 2, 2);
1633   gsl_vector_view sv = gsl_vector_view_array(datS, 2);
1634   gsl_vector_view wv = gsl_vector_view_array(work, 2);
1635   double B11, B22;
1636   size_t top;
1637   double det;
1638   double cr, sr, cl, sl;
1639   gsl_vector_view xv, yv;
1640   int s;
1641
1642   if (w->needtop)
1643     top = gen_get_submatrix(w->H, A);
1644
1645   /*
1646    * Rotate B so that
1647    *
1648    * B = [ B11  0  ]
1649    *     [  0  B22 ]
1650    *
1651    * with B11 non-negative
1652    */
1653
1654   gsl_matrix_memcpy(&uv.matrix, B);
1655   gsl_linalg_SV_decomp(&uv.matrix, &vv.matrix, &sv.vector, &wv.vector);
1656
1657   /*
1658    * Right now, B = U S V^t, where S = diag(s)
1659    *
1660    * The SVD routine may have computed reflection matrices U and V,
1661    * but it would be much nicer to have rotations since we won't have
1662    * to use BLAS mat-mat multiplications to update our matrices,
1663    * and can instead use drot. So convert them to rotations if
1664    * necessary
1665    */
1666
1667   det = gsl_matrix_get(&vv.matrix, 0, 0) * gsl_matrix_get(&vv.matrix, 1, 1) -
1668         gsl_matrix_get(&vv.matrix, 0, 1) * gsl_matrix_get(&vv.matrix, 1, 0);
1669   if (det < 0.0)
1670     {
1671       /* V is a reflection, convert it to a rotation by inserting
1672        * F = [1 0; 0 -1] so that:
1673        *
1674        * B = U S [1  0] [1  0] V^t
1675        *         [0 -1] [0 -1]
1676        *
1677        * so S -> S F and V -> V F where F is the reflection matrix
1678        * We just need to invert S22 since the first column of V
1679        * will remain unchanged and we can just read off the CS and SN
1680        * parameters.
1681        */
1682       datS[1] = -datS[1];
1683     }
1684
1685   cr = gsl_matrix_get(&vv.matrix, 0, 0);
1686   sr = gsl_matrix_get(&vv.matrix, 1, 0);
1687
1688   /* same for U */
1689   det = gsl_matrix_get(&uv.matrix, 0, 0) * gsl_matrix_get(&uv.matrix, 1, 1) -
1690         gsl_matrix_get(&uv.matrix, 0, 1) * gsl_matrix_get(&uv.matrix, 1, 0);
1691   if (det < 0.0)
1692     datS[1] = -datS[1];
1693
1694   cl = gsl_matrix_get(&uv.matrix, 0, 0);
1695   sl = gsl_matrix_get(&uv.matrix, 1, 0);
1696
1697   B11 = gsl_vector_get(&sv.vector, 0);
1698   B22 = gsl_vector_get(&sv.vector, 1);
1699
1700   /* make sure B11 is positive */
1701   if (B11 < 0.0)
1702     {
1703       B11 = -B11;
1704       B22 = -B22;
1705       cr = -cr;
1706       sr = -sr;
1707     }
1708
1709   /*
1710    * At this point,
1711    *
1712    * [ S11  0  ] = [ CSL  SNL ] B [ CSR -SNR ]
1713    * [  0  S22 ]   [-SNL  CSL ]   [ SNR  CSR ]
1714    *
1715    * apply rotations to H and rest of R
1716    */
1717
1718   if (w->compute_s)
1719     {
1720       xv = gsl_matrix_subrow(w->H, top, top, w->size - top);
1721       yv = gsl_matrix_subrow(w->H, top + 1, top, w->size - top);
1722       gsl_blas_drot(&xv.vector, &yv.vector, cl, sl);
1723
1724       xv = gsl_matrix_subcolumn(w->H, top, 0, top + 2);
1725       yv = gsl_matrix_subcolumn(w->H, top + 1, 0, top + 2);
1726       gsl_blas_drot(&xv.vector, &yv.vector, cr, sr);
1727     }
1728   else
1729     {
1730       xv = gsl_matrix_row(A, 0);
1731       yv = gsl_matrix_row(A, 1);
1732       gsl_blas_drot(&xv.vector, &yv.vector, cl, sl);
1733
1734       xv = gsl_matrix_column(A, 0);
1735       yv = gsl_matrix_column(A, 1);
1736       gsl_blas_drot(&xv.vector, &yv.vector, cr, sr);
1737     }
1738
1739   if (w->compute_t)
1740     {
1741       if (top != (w->size - 2))
1742         {
1743           xv = gsl_matrix_subrow(w->R, top, top + 2, w->size - top - 2);
1744           yv = gsl_matrix_subrow(w->R, top + 1, top + 2, w->size - top - 2);
1745           gsl_blas_drot(&xv.vector, &yv.vector, cl, sl);
1746         }
1747
1748       if (top != 0)
1749         {
1750           xv = gsl_matrix_subcolumn(w->R, top, 0, top);
1751           yv = gsl_matrix_subcolumn(w->R, top + 1, 0, top);
1752           gsl_blas_drot(&xv.vector, &yv.vector, cr, sr);
1753         }
1754     }
1755
1756   if (w->Q)
1757     {
1758       xv = gsl_matrix_column(w->Q, top);
1759       yv = gsl_matrix_column(w->Q, top + 1);
1760       gsl_blas_drot(&xv.vector, &yv.vector, cl, sl);
1761     }
1762
1763   if (w->Z)
1764     {
1765       xv = gsl_matrix_column(w->Z, top);
1766       yv = gsl_matrix_column(w->Z, top + 1);
1767       gsl_blas_drot(&xv.vector, &yv.vector, cr, sr);
1768     }
1769
1770   gsl_matrix_set(B, 0, 0, B11);
1771   gsl_matrix_set(B, 0, 1, 0.0);
1772   gsl_matrix_set(B, 1, 0, 0.0);
1773   gsl_matrix_set(B, 1, 1, B22);
1774
1775   /* if B22 is < 0, make it positive by negating its column */
1776   if (B22 < 0.0)
1777     {
1778       size_t i;
1779
1780       if (w->compute_s)
1781         {
1782           for (i = 0; i < top + 2; ++i)
1783             gsl_matrix_set(w->H, i, top + 1, -gsl_matrix_get(w->H, i, top + 1));
1784         }
1785       else
1786         {
1787           gsl_matrix_set(A, 0, 1, -gsl_matrix_get(A, 0, 1));
1788           gsl_matrix_set(A, 1, 1, -gsl_matrix_get(A, 1, 1));
1789         }
1790
1791       if (w->compute_t)
1792         {
1793           for (i = 0; i < top + 2; ++i)
1794             gsl_matrix_set(w->R, i, top + 1, -gsl_matrix_get(w->R, i, top + 1));
1795         }
1796       else
1797         {
1798           gsl_matrix_set(B, 0, 1, -gsl_matrix_get(B, 0, 1));
1799           gsl_matrix_set(B, 1, 1, -gsl_matrix_get(B, 1, 1));
1800         }
1801
1802       if (w->Z)
1803         {
1804           xv = gsl_matrix_column(w->Z, top + 1);
1805           gsl_vector_scale(&xv.vector, -1.0);
1806         }
1807     }
1808
1809   /* our block is now in standard form - compute eigenvalues */
1810
1811   s = gen_compute_eigenvals(A, B, alpha1, alpha2, beta1, beta2);
1812
1813   return s;
1814 } /* gen_schur_standardize2() */
1815
1816 /*
1817 gen_compute_eigenvals()
1818   Compute the complex eigenvalues of a 2-by-2 block
1819
1820 Return: GSL_CONTINUE if block contains real eigenvalues (they are not
1821         computed)
1822         GSL_SUCCESS on normal completion
1823 */
1824
1825 static int
1826 gen_compute_eigenvals(gsl_matrix *A, gsl_matrix *B, gsl_complex *alpha1,
1827                       gsl_complex *alpha2, double *beta1, double *beta2)
1828 {
1829   double wr1, wr2, wi, scale1, scale2;
1830   double s1inv;
1831   double A11, A12, A21, A22;
1832   double B11, B22;
1833   double c11r, c11i, c12, c21, c22r, c22i;
1834   double cz, cq;
1835   double szr, szi, sqr, sqi;
1836   double a1r, a1i, a2r, a2i, b1r, b1i, b1a, b2r, b2i, b2a;
1837   double alphar, alphai;
1838   double t1, an, bn, tempr, tempi, wabs;
1839
1840   /*
1841    * This function is called from gen_schur_standardize2() and
1842    * its possible the standardization has perturbed the eigenvalues
1843    * onto the real line - so check for this before computing them
1844    */
1845
1846   gsl_schur_gen_eigvals(A, B, &wr1, &wr2, &wi, &scale1, &scale2);
1847   if (wi == 0.0)
1848     return GSL_CONTINUE; /* real eigenvalues - continue QZ iteration */
1849
1850   /* complex eigenvalues - compute alpha and beta */
1851
1852   s1inv = 1.0 / scale1;
1853
1854   A11 = gsl_matrix_get(A, 0, 0);
1855   A12 = gsl_matrix_get(A, 0, 1);
1856   A21 = gsl_matrix_get(A, 1, 0);
1857   A22 = gsl_matrix_get(A, 1, 1);
1858
1859   B11 = gsl_matrix_get(B, 0, 0);
1860   B22 = gsl_matrix_get(B, 1, 1);
1861
1862   c11r = scale1 * A11 - wr1 * B11;
1863   c11i = -wi * B11;
1864   c12 = scale1 * A12;
1865   c21 = scale1 * A21;
1866   c22r = scale1 * A22 - wr1 * B22;
1867   c22i = -wi * B22;
1868
1869   if (fabs(c11r) + fabs(c11i) + fabs(c12) >
1870       fabs(c21) + fabs(c22r) + fabs(c22i))
1871     {
1872       t1 = gsl_hypot3(c12, c11r, c11i);
1873       if (t1 != 0.0)
1874         {
1875           cz = c12 / t1;
1876           szr = -c11r / t1;
1877           szi = -c11i / t1;
1878         }
1879       else
1880         {
1881           cz = 0.0;
1882           szr = 1.0;
1883           szi = 0.0;
1884         }
1885     }
1886   else
1887     {
1888       cz = hypot(c22r, c22i);
1889       if (cz <= GSL_DBL_MIN)
1890         {
1891           cz = 0.0;
1892           szr = 1.0;
1893           szi = 0.0;
1894         }
1895       else
1896         {
1897           tempr = c22r / cz;
1898           tempi = c22i / cz;
1899           t1 = hypot(cz, c21);
1900           cz /= t1;
1901           szr = -c21*tempr / t1;
1902           szi = c21*tempi / t1;
1903         }
1904     }
1905
1906   an = fabs(A11) + fabs(A12) + fabs(A21) + fabs(A22);
1907   bn = fabs(B11) + fabs(B22);
1908   wabs = fabs(wr1) + fabs(wi);
1909   if (scale1*an > wabs*bn)
1910     {
1911       cq = cz * B11;
1912       if (cq <= GSL_DBL_MIN)
1913         {
1914           cq = 0.0;
1915           sqr = 1.0;
1916           sqi = 0.0;
1917         }
1918       else
1919         {
1920           sqr = szr * B22;
1921           sqi = -szi * B22;
1922         }
1923     }
1924   else
1925     {
1926       a1r = cz * A11 + szr * A12;
1927       a1i = szi * A12;
1928       a2r = cz * A21 + szr * A22;
1929       a2i = szi * A22;
1930       cq = hypot(a1r, a1i);
1931       if (cq <= GSL_DBL_MIN)
1932         {
1933           cq = 0.0;
1934           sqr = 1.0;
1935           sqi = 0.0;
1936         }
1937       else
1938         {
1939           tempr = a1r / cq;
1940           tempi = a1i / cq;
1941           sqr = tempr * a2r + tempi * a2i;
1942           sqi = tempi * a2r - tempr * a2i;
1943         }
1944     }
1945
1946   t1 = gsl_hypot3(cq, sqr, sqi);
1947   cq /= t1;
1948   sqr /= t1;
1949   sqi /= t1;
1950
1951   tempr = sqr*szr - sqi*szi;
1952   tempi = sqr*szi + sqi*szr;
1953   b1r = cq*cz*B11 + tempr*B22;
1954   b1i = tempi*B22;
1955   b1a = hypot(b1r, b1i);
1956   b2r = cq*cz*B22 + tempr*B11;
1957   b2i = -tempi*B11;
1958   b2a = hypot(b2r, b2i);
1959
1960   *beta1 = b1a;
1961   *beta2 = b2a;
1962
1963   alphar = (wr1 * b1a) * s1inv;
1964   alphai = (wi * b1a) * s1inv;
1965   GSL_SET_COMPLEX(alpha1, alphar, alphai);
1966
1967   alphar = (wr1 * b2a) * s1inv;
1968   alphai = -(wi * b2a) * s1inv;
1969   GSL_SET_COMPLEX(alpha2, alphar, alphai);
1970
1971   return GSL_SUCCESS;
1972 } /* gen_compute_eigenvals() */
1973
1974 /*
1975 gen_store_eigval1()
1976   Store eigenvalue of a 1-by-1 block into the alpha and beta
1977 output vectors. This routine ensures that eigenvalues are stored
1978 in the same order as they appear in the Schur form and updates
1979 various internal workspace quantities.
1980 */
1981
1982 static void
1983 gen_store_eigval1(const gsl_matrix *H, const double a, const double b,
1984                   gsl_vector_complex *alpha,
1985                   gsl_vector *beta, gsl_eigen_gen_workspace *w)
1986 {
1987   size_t top = gen_get_submatrix(w->H, H);
1988   gsl_complex z;
1989
1990   GSL_SET_COMPLEX(&z, a, 0.0);
1991
1992   gsl_vector_complex_set(alpha, top, z);
1993   gsl_vector_set(beta, top, b);
1994
1995   w->n_evals += 1;
1996   w->n_iter = 0;
1997   w->eshift = 0.0;
1998 } /* gen_store_eigval1() */
1999
2000 /*
2001 gen_store_eigval2()
2002   Store eigenvalues of a 2-by-2 block into the alpha and beta
2003 output vectors. This routine ensures that eigenvalues are stored
2004 in the same order as they appear in the Schur form and updates
2005 various internal workspace quantities.
2006 */
2007
2008 static void
2009 gen_store_eigval2(const gsl_matrix *H, const gsl_complex *alpha1,
2010                   const double beta1, const gsl_complex *alpha2,
2011                   const double beta2, gsl_vector_complex *alpha,
2012                   gsl_vector *beta, gsl_eigen_gen_workspace *w)
2013 {
2014   size_t top = gen_get_submatrix(w->H, H);
2015
2016   gsl_vector_complex_set(alpha, top, *alpha1);
2017   gsl_vector_set(beta, top, beta1);
2018
2019   gsl_vector_complex_set(alpha, top + 1, *alpha2);
2020   gsl_vector_set(beta, top + 1, beta2);
2021
2022   w->n_evals += 2;
2023   w->n_iter = 0;
2024   w->eshift = 0.0;
2025 } /* gen_store_eigval2() */
2026
2027 /*
2028 gen_get_submatrix()
2029   B is a submatrix of A. The goal of this function is to
2030 compute the indices in A of where the matrix B resides
2031 */
2032
2033 static inline size_t
2034 gen_get_submatrix(const gsl_matrix *A, const gsl_matrix *B)
2035 {
2036   size_t diff;
2037   double ratio;
2038   size_t top;
2039
2040   diff = (size_t) (B->data - A->data);
2041
2042   /* B is on the diagonal of A, so measure distance in units of
2043      tda+1 */
2044
2045   ratio = (double)diff / ((double) (A->tda + 1));
2046
2047   top = (size_t) floor(ratio);
2048
2049   return top;
2050 } /* gen_get_submatrix() */
2051
2052 /* Frobenius norm */
2053 inline static double
2054 normF (gsl_matrix * A)
2055 {
2056   size_t i, j, M = A->size1, N = A->size2;
2057   double sum = 0.0, scale = 0.0, ssq = 1.0;
2058
2059   for (i = 0; i < M; i++)
2060     {
2061       for (j = 0; j < N; j++)
2062         {
2063           double Aij = gsl_matrix_get (A, i, j);
2064
2065           if (Aij != 0.0)
2066             {
2067               double ax = fabs (Aij);
2068
2069               if (scale < ax)
2070                 {
2071                   ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2072                   scale = ax;
2073                 }
2074               else
2075                 {
2076                   ssq += (ax / scale) * (ax / scale);
2077                 }
2078             }
2079
2080         }
2081     }
2082
2083   sum = scale * sqrt (ssq);
2084
2085   return sum;
2086 }
2087
2088 /* Generate a Givens rotation (cos,sin) which takes v=(x,y) to (|v|,0)
2089
2090    From Golub and Van Loan, "Matrix Computations", Section 5.1.8 */
2091
2092 inline static void
2093 create_givens (const double a, const double b, double *c, double *s)
2094 {
2095   if (b == 0)
2096     {
2097       *c = 1;
2098       *s = 0;
2099     }
2100   else if (fabs (b) > fabs (a))
2101     {
2102       double t = -a / b;
2103       double s1 = 1.0 / sqrt (1 + t * t);
2104       *s = s1;
2105       *c = s1 * t;
2106     }
2107   else
2108     {
2109       double t = -b / a;
2110       double c1 = 1.0 / sqrt (1 + t * t);
2111       *c = c1;
2112       *s = c1 * t;
2113     }
2114 }