Annotated Source Code

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