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

   1 /* eigen/genv.c
   2  *
   3  * Copyright (C) 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 #include <gsl/gsl_errno.h>
  32
  33 /*
  34  * This module computes the eigenvalues and eigenvectors of a
  35  * real generalized eigensystem A x = \lambda B x. Left and right
  36  * Schur vectors are optionally computed as well.
  37  *
  38  * This file contains routines based on original code from LAPACK
  39  * which is distributed under the modified BSD license.
  40  */
  41
  42 static int genv_get_right_eigenvectors(const gsl_matrix *S,
  43                                        const gsl_matrix *T,
  44                                        gsl_matrix *Z,
  45                                        gsl_matrix_complex *evec,
  46                                        gsl_eigen_genv_workspace *w);
  47 static void genv_normalize_eigenvectors(gsl_vector_complex *alpha,
  48                                         gsl_matrix_complex *evec);
  49
  50 /*
  51 gsl_eigen_genv_alloc()
  52   Allocate a workspace for solving the generalized eigenvalue problem.
  53 The size of this workspace is O(7n).
  54
  55 Inputs: n - size of matrices
  56
  57 Return: pointer to workspace
  58 */
  59
  60 gsl_eigen_genv_workspace *
  61 gsl_eigen_genv_alloc(const size_t n)
  62 {
  63   gsl_eigen_genv_workspace *w;
  64
  65   if (n == 0)
  66     {
  67       GSL_ERROR_NULL ("matrix dimension must be positive integer",
  68                       GSL_EINVAL);
  69     }
  70
  71   w = calloc (1, sizeof (gsl_eigen_genv_workspace));
  72
  73   if (w == 0)
  74     {
  75       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
  76     }
  77
  78   w->size = n;
  79   w->Q = NULL;
  80   w->Z = NULL;
  81
  82   w->gen_workspace_p = gsl_eigen_gen_alloc(n);
  83
  84   if (w->gen_workspace_p == 0)
  85     {
  86       gsl_eigen_genv_free(w);
  87       GSL_ERROR_NULL ("failed to allocate space for gen workspace", GSL_ENOMEM);
  88     }
  89
  90   /* compute the full Schur forms */
  91   gsl_eigen_gen_params(1, 1, 1, w->gen_workspace_p);
  92
  93   w->work1 = gsl_vector_alloc(n);
  94   w->work2 = gsl_vector_alloc(n);
  95   w->work3 = gsl_vector_alloc(n);
  96   w->work4 = gsl_vector_alloc(n);
  97   w->work5 = gsl_vector_alloc(n);
  98   w->work6 = gsl_vector_alloc(n);
  99
 100   if (w->work1 == 0 || w->work2 == 0 || w->work3 == 0 ||
 101       w->work4 == 0 || w->work5 == 0 || w->work6 == 0)
 102     {
 103       gsl_eigen_genv_free(w);
 104       GSL_ERROR_NULL ("failed to allocate space for additional workspace", GSL_ENOMEM);
 105     }
 106
 107   return (w);
 108 } /* gsl_eigen_genv_alloc() */
 109
 110 /*
 111 gsl_eigen_genv_free()
 112   Free workspace w
 113 */
 114
 115 void
 116 gsl_eigen_genv_free(gsl_eigen_genv_workspace *w)
 117 {
 118   if (w->gen_workspace_p)
 119     gsl_eigen_gen_free(w->gen_workspace_p);
 120
 121   if (w->work1)
 122     gsl_vector_free(w->work1);
 123
 124   if (w->work2)
 125     gsl_vector_free(w->work2);
 126
 127   if (w->work3)
 128     gsl_vector_free(w->work3);
 129
 130   if (w->work4)
 131     gsl_vector_free(w->work4);
 132
 133   if (w->work5)
 134     gsl_vector_free(w->work5);
 135
 136   if (w->work6)
 137     gsl_vector_free(w->work6);
 138
 139   free(w);
 140 } /* gsl_eigen_genv_free() */
 141
 142 /*
 143 gsl_eigen_genv()
 144
 145 Solve the generalized eigenvalue problem
 146
 147 A x = \lambda B x
 148
 149 for the eigenvalues \lambda and right eigenvectors x.
 150
 151 Inputs: A     - general real matrix
 152         B     - general real matrix
 153         alpha - (output) where to store eigenvalue numerators
 154         beta  - (output) where to store eigenvalue denominators
 155         evec  - (output) where to store eigenvectors
 156         w     - workspace
 157
 158 Return: success or error
 159 */
 160
 161 int
 162 gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B, gsl_vector_complex * alpha,
 163                 gsl_vector * beta, gsl_matrix_complex *evec,
 164                 gsl_eigen_genv_workspace * w)
 165 {
 166   const size_t N = A->size1;
 167
 168   /* check matrix and vector sizes */
 169
 170   if (N != A->size2)
 171     {
 172       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 173     }
 174   else if ((N != B->size1) || (N != B->size2))
 175     {
 176       GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
 177     }
 178   else if (alpha->size != N || beta->size != N)
 179     {
 180       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 181     }
 182   else if (w->size != N)
 183     {
 184       GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
 185     }
 186   else if (evec->size1 != N)
 187     {
 188       GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
 189     }
 190   else
 191     {
 192       int s;
 193       gsl_matrix Z;
 194
 195       /*
 196        * We need a place to store the right Schur vectors, so we will
 197        * treat evec as a real matrix and store them in the left
 198        * half - the factor of 2 in the tda corresponds to the
 199        * complex multiplicity
 200        */
 201       Z.size1 = N;
 202       Z.size2 = N;
 203       Z.tda = 2 * N;
 204       Z.data = evec->data;
 205       Z.block = 0;
 206       Z.owner = 0;
 207
 208       s = gsl_eigen_gen_QZ(A, B, alpha, beta, w->Q, &Z, w->gen_workspace_p);
 209
 210       if (w->Z)
 211         {
 212           /* save right Schur vectors */
 213           gsl_matrix_memcpy(w->Z, &Z);
 214         }
 215
 216       /* only compute eigenvectors if we found all eigenvalues */
 217       if (s == GSL_SUCCESS)
 218         {
 219           /* compute eigenvectors */
 220           s = genv_get_right_eigenvectors(A, B, &Z, evec, w);
 221
 222           if (s == GSL_SUCCESS)
 223             genv_normalize_eigenvectors(alpha, evec);
 224         }
 225
 226       return s;
 227     }
 228 } /* gsl_eigen_genv() */
 229
 230 /*
 231 gsl_eigen_genv_QZ()
 232
 233 Solve the generalized eigenvalue problem
 234
 235 A x = \lambda B x
 236
 237 for the eigenvalues \lambda and right eigenvectors x. Optionally
 238 compute left and/or right Schur vectors Q and Z which satisfy:
 239
 240 A = Q S Z^t
 241 B = Q T Z^t
 242
 243 where (S, T) is the generalized Schur form of (A, B)
 244
 245 Inputs: A     - general real matrix
 246         B     - general real matrix
 247         alpha - (output) where to store eigenvalue numerators
 248         beta  - (output) where to store eigenvalue denominators
 249         evec  - (output) where to store eigenvectors
 250         Q     - (output) if non-null, where to store left Schur vectors
 251         Z     - (output) if non-null, where to store right Schur vectors
 252         w     - workspace
 253
 254 Return: success or error
 255 */
 256
 257 int
 258 gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
 259                    gsl_vector_complex * alpha, gsl_vector * beta,
 260                    gsl_matrix_complex * evec,
 261                    gsl_matrix * Q, gsl_matrix * Z,
 262                    gsl_eigen_genv_workspace * w)
 263 {
 264   if (Q && (A->size1 != Q->size1 || A->size1 != Q->size2))
 265     {
 266       GSL_ERROR("Q matrix has wrong dimensions", GSL_EBADLEN);
 267     }
 268   else if (Z && (A->size1 != Z->size1 || A->size1 != Z->size2))
 269     {
 270       GSL_ERROR("Z matrix has wrong dimensions", GSL_EBADLEN);
 271     }
 272   else
 273     {
 274       int s;
 275
 276       w->Q = Q;
 277       w->Z = Z;
 278
 279       s = gsl_eigen_genv(A, B, alpha, beta, evec, w);
 280
 281       w->Q = NULL;
 282       w->Z = NULL;
 283
 284       return s;
 285     }
 286 } /* gsl_eigen_genv_QZ() */
 287
 288 /********************************************
 289  *           INTERNAL ROUTINES              *
 290  ********************************************/
 291
 292 /*
 293 genv_get_right_eigenvectors()
 294   Compute right eigenvectors of the Schur form (S, T) and then
 295 backtransform them using the right Schur vectors to get right
 296 eigenvectors of the original system.
 297
 298 Inputs: S     - upper quasi-triangular Schur form of A
 299         T     - upper triangular Schur form of B
 300         Z     - right Schur vectors
 301         evec  - (output) where to store eigenvectors
 302         w     - workspace
 303
 304 Return: success or error
 305
 306 Notes: 1) based on LAPACK routine DTGEVC
 307        2) eigenvectors are stored in the order that their
 308           eigenvalues appear in the Schur form
 309 */
 310
 311 static int
 312 genv_get_right_eigenvectors(const gsl_matrix *S, const gsl_matrix *T,
 313                             gsl_matrix *Z,
 314                             gsl_matrix_complex *evec,
 315                             gsl_eigen_genv_workspace *w)
 316 {
 317   const size_t N = w->size;
 318   const double small = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
 319   const double big = 1.0 / small;
 320   const double bignum = 1.0 / (GSL_DBL_MIN * N);
 321   size_t i, j, k, end;
 322   int is;
 323   double anorm, bnorm;
 324   double temp, temp2, temp2r, temp2i;
 325   double ascale, bscale;
 326   double salfar, sbeta;
 327   double acoef, bcoefr, bcoefi, acoefa, bcoefa;
 328   double creala, cimaga, crealb, cimagb, cre2a, cim2a, cre2b, cim2b;
 329   double dmin, xmax;
 330   double scale;
 331   size_t nw, na;
 332   int lsa, lsb;
 333   int complex_pair;
 334   gsl_complex z_zero, z_one;
 335   double bdiag[2] = { 0.0, 0.0 };
 336   double sum[4];
 337   int il2by2;
 338   size_t jr, jc, ja;
 339   double xscale;
 340   gsl_vector_complex_view ecol;
 341   gsl_vector_view re, im, re2, im2;
 342
 343   GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);
 344   GSL_SET_COMPLEX(&z_one, 1.0, 0.0);
 345
 346   /*
 347    * Compute the 1-norm of each column of (S, T) excluding elements
 348    * belonging to the diagonal blocks to check for possible overflow
 349    * in the triangular solver
 350    */
 351
 352   anorm = fabs(gsl_matrix_get(S, 0, 0));
 353   if (N > 1)
 354     anorm += fabs(gsl_matrix_get(S, 1, 0));
 355   bnorm = fabs(gsl_matrix_get(T, 0, 0));
 356
 357   gsl_vector_set(w->work1, 0, 0.0);
 358   gsl_vector_set(w->work2, 0, 0.0);
 359
 360   for (j = 1; j < N; ++j)
 361     {
 362       temp = temp2 = 0.0;
 363       if (gsl_matrix_get(S, j, j - 1) == 0.0)
 364         end = j;
 365       else
 366         end = j - 1;
 367
 368       for (i = 0; i < end; ++i)
 369         {
 370           temp += fabs(gsl_matrix_get(S, i, j));
 371           temp2 += fabs(gsl_matrix_get(T, i, j));
 372         }
 373
 374       gsl_vector_set(w->work1, j, temp);
 375       gsl_vector_set(w->work2, j, temp2);
 376
 377       for (i = end; i < GSL_MIN(j + 2, N); ++i)
 378         {
 379           temp += fabs(gsl_matrix_get(S, i, j));
 380           temp2 += fabs(gsl_matrix_get(T, i, j));
 381         }
 382
 383       anorm = GSL_MAX(anorm, temp);
 384       bnorm = GSL_MAX(bnorm, temp2);
 385     }
 386
 387   ascale = 1.0 / GSL_MAX(anorm, GSL_DBL_MIN);
 388   bscale = 1.0 / GSL_MAX(bnorm, GSL_DBL_MIN);
 389
 390   complex_pair = 0;
 391   for (k = 0; k < N; ++k)
 392     {
 393       size_t je = N - 1 - k;
 394
 395       if (complex_pair)
 396         {
 397           complex_pair = 0;
 398           continue;
 399         }
 400
 401       nw = 1;
 402       if (je > 0)
 403         {
 404           if (gsl_matrix_get(S, je, je - 1) != 0.0)
 405             {
 406               complex_pair = 1;
 407               nw = 2;
 408             }
 409         }
 410
 411       if (!complex_pair)
 412         {
 413           if (fabs(gsl_matrix_get(S, je, je)) <= GSL_DBL_MIN &&
 414               fabs(gsl_matrix_get(T, je, je)) <= GSL_DBL_MIN)
 415             {
 416               /* singular matrix pencil - unit eigenvector */
 417               for (i = 0; i < N; ++i)
 418                 gsl_matrix_complex_set(evec, i, je, z_zero);
 419
 420               gsl_matrix_complex_set(evec, je, je, z_one);
 421
 422               continue;
 423             }
 424
 425           /* clear vector */
 426           for (i = 0; i < N; ++i)
 427             gsl_vector_set(w->work3, i, 0.0);
 428         }
 429       else
 430         {
 431           /* clear vectors */
 432           for (i = 0; i < N; ++i)
 433             {
 434               gsl_vector_set(w->work3, i, 0.0);
 435               gsl_vector_set(w->work4, i, 0.0);
 436             }
 437         }
 438
 439       if (!complex_pair)
 440         {
 441           /* real eigenvalue */
 442
 443           temp = 1.0 / GSL_MAX(GSL_DBL_MIN,
 444                                GSL_MAX(fabs(gsl_matrix_get(S, je, je)) * ascale,
 445                                        fabs(gsl_matrix_get(T, je, je)) * bscale));
 446           salfar = (temp * gsl_matrix_get(S, je, je)) * ascale;
 447           sbeta = (temp * gsl_matrix_get(T, je, je)) * bscale;
 448           acoef = sbeta * ascale;
 449           bcoefr = salfar * bscale;
 450           bcoefi = 0.0;
 451
 452           /* scale to avoid underflow */
 453           scale = 1.0;
 454           lsa = fabs(sbeta) >= GSL_DBL_MIN && fabs(acoef) < small;
 455           lsb = fabs(salfar) >= GSL_DBL_MIN && fabs(bcoefr) < small;
 456           if (lsa)
 457             scale = (small / fabs(sbeta)) * GSL_MIN(anorm, big);
 458           if (lsb)
 459             scale = GSL_MAX(scale, (small / fabs(salfar)) * GSL_MIN(bnorm, big));
 460
 461           if (lsa || lsb)
 462             {
 463               scale = GSL_MIN(scale,
 464                         1.0 / (GSL_DBL_MIN *
 465                                GSL_MAX(1.0,
 466                                  GSL_MAX(fabs(acoef), fabs(bcoefr)))));
 467               if (lsa)
 468                 acoef = ascale * (scale * sbeta);
 469               else
 470                 acoef *= scale;
 471
 472               if (lsb)
 473                 bcoefr = bscale * (scale * salfar);
 474               else
 475                 bcoefr *= scale;
 476             }
 477
 478           acoefa = fabs(acoef);
 479           bcoefa = fabs(bcoefr);
 480
 481           /* first component is 1 */
 482           gsl_vector_set(w->work3, je, 1.0);
 483           xmax = 1.0;
 484
 485           /* compute contribution from column je of A and B to sum */
 486
 487           for (i = 0; i < je; ++i)
 488             {
 489               gsl_vector_set(w->work3, i,
 490                 bcoefr*gsl_matrix_get(T, i, je) -
 491                 acoef * gsl_matrix_get(S, i, je));
 492             }
 493         }
 494       else
 495         {
 496           gsl_matrix_const_view vs =
 497             gsl_matrix_const_submatrix(S, je - 1, je - 1, 2, 2);
 498           gsl_matrix_const_view vt =
 499             gsl_matrix_const_submatrix(T, je - 1, je - 1, 2, 2);
 500
 501           /* complex eigenvalue */
 502
 503           gsl_schur_gen_eigvals(&vs.matrix,
 504                                 &vt.matrix,
 505                                 &bcoefr,
 506                                 &temp2,
 507                                 &bcoefi,
 508                                 &acoef,
 509                                 &temp);
 510           if (bcoefi == 0.0)
 511             {
 512               GSL_ERROR("gsl_schur_gen_eigvals failed on complex block", GSL_FAILURE);
 513             }
 514
 515           /* scale to avoid over/underflow */
 516           acoefa = fabs(acoef);
 517           bcoefa = fabs(bcoefr) + fabs(bcoefi);
 518           scale = 1.0;
 519
 520           if (acoefa*GSL_DBL_EPSILON < GSL_DBL_MIN && acoefa >= GSL_DBL_MIN)
 521             scale = (GSL_DBL_MIN / GSL_DBL_EPSILON) / acoefa;
 522           if (bcoefa*GSL_DBL_EPSILON < GSL_DBL_MIN && bcoefa >= GSL_DBL_MIN)
 523             scale = GSL_MAX(scale, (GSL_DBL_MIN/GSL_DBL_EPSILON) / bcoefa);
 524           if (GSL_DBL_MIN*acoefa > ascale)
 525             scale = ascale / (GSL_DBL_MIN * acoefa);
 526           if (GSL_DBL_MIN*bcoefa > bscale)
 527             scale = GSL_MIN(scale, bscale / (GSL_DBL_MIN*bcoefa));
 528           if (scale != 1.0)
 529             {
 530               acoef *= scale;
 531               acoefa = fabs(acoef);
 532               bcoefr *= scale;
 533               bcoefi *= scale;
 534               bcoefa = fabs(bcoefr) + fabs(bcoefi);
 535             }
 536
 537           /* compute first two components of eigenvector */
 538
 539           temp = acoef * gsl_matrix_get(S, je, je - 1);
 540           temp2r = acoef * gsl_matrix_get(S, je, je) -
 541                    bcoefr * gsl_matrix_get(T, je, je);
 542           temp2i = -bcoefi * gsl_matrix_get(T, je, je);
 543
 544           if (fabs(temp) >= fabs(temp2r) + fabs(temp2i))
 545             {
 546               gsl_vector_set(w->work3, je, 1.0);
 547               gsl_vector_set(w->work4, je, 0.0);
 548               gsl_vector_set(w->work3, je - 1, -temp2r / temp);
 549               gsl_vector_set(w->work4, je - 1, -temp2i / temp);
 550             }
 551           else
 552             {
 553               gsl_vector_set(w->work3, je - 1, 1.0);
 554               gsl_vector_set(w->work4, je - 1, 0.0);
 555               temp = acoef * gsl_matrix_get(S, je - 1, je);
 556               gsl_vector_set(w->work3, je,
 557                 (bcoefr*gsl_matrix_get(T, je - 1, je - 1) -
 558                  acoef*gsl_matrix_get(S, je - 1, je - 1)) / temp);
 559               gsl_vector_set(w->work4, je,
 560                 bcoefi*gsl_matrix_get(T, je - 1, je - 1) / temp);
 561             }
 562
 563           xmax = GSL_MAX(fabs(gsl_vector_get(w->work3, je)) +
 564                          fabs(gsl_vector_get(w->work4, je)),
 565                          fabs(gsl_vector_get(w->work3, je - 1)) +
 566                          fabs(gsl_vector_get(w->work4, je - 1)));
 567
 568           /* compute contribution from column je and je - 1 */
 569
 570           creala = acoef * gsl_vector_get(w->work3, je - 1);
 571           cimaga = acoef * gsl_vector_get(w->work4, je - 1);
 572           crealb = bcoefr * gsl_vector_get(w->work3, je - 1) -
 573                    bcoefi * gsl_vector_get(w->work4, je - 1);
 574           cimagb = bcoefi * gsl_vector_get(w->work3, je - 1) +
 575                    bcoefr * gsl_vector_get(w->work4, je - 1);
 576           cre2a = acoef * gsl_vector_get(w->work3, je);
 577           cim2a = acoef * gsl_vector_get(w->work4, je);
 578           cre2b = bcoefr * gsl_vector_get(w->work3, je) -
 579                   bcoefi * gsl_vector_get(w->work4, je);
 580           cim2b = bcoefi * gsl_vector_get(w->work3, je) +
 581                   bcoefr * gsl_vector_get(w->work4, je);
 582
 583           for (i = 0; i < je - 1; ++i)
 584             {
 585               gsl_vector_set(w->work3, i,
 586                 -creala * gsl_matrix_get(S, i, je - 1) +
 587                 crealb * gsl_matrix_get(T, i, je - 1) -
 588                 cre2a * gsl_matrix_get(S, i, je) +
 589                 cre2b * gsl_matrix_get(T, i, je));
 590               gsl_vector_set(w->work4, i,
 591                 -cimaga * gsl_matrix_get(S, i, je - 1) +
 592                 cimagb * gsl_matrix_get(T, i, je - 1) -
 593                 cim2a * gsl_matrix_get(S, i, je) +
 594                 cim2b * gsl_matrix_get(T, i, je));
 595             }
 596         }
 597
 598       dmin = GSL_MAX(GSL_DBL_MIN,
 599                GSL_MAX(GSL_DBL_EPSILON*acoefa*anorm,
 600                        GSL_DBL_EPSILON*bcoefa*bnorm));
 601
 602       /* triangular solve of (a A - b B) x = 0 */
 603
 604       il2by2 = 0;
 605       for (is = (int) je - (int) nw; is >= 0; --is)
 606         {
 607           j = (size_t) is;
 608
 609           if (!il2by2 && j > 0)
 610             {
 611               if (gsl_matrix_get(S, j, j - 1) != 0.0)
 612                 {
 613                   il2by2 = 1;
 614                   continue;
 615                 }
 616             }
 617
 618           bdiag[0] = gsl_matrix_get(T, j, j);
 619           if (il2by2)
 620             {
 621               na = 2;
 622               bdiag[1] = gsl_matrix_get(T, j + 1, j + 1);
 623             }
 624           else
 625             na = 1;
 626
 627
 628           if (nw == 1)
 629             {
 630               gsl_matrix_const_view sv =
 631                 gsl_matrix_const_submatrix(S, j, j, na, na);
 632               gsl_vector_view xv, bv;
 633
 634               bv = gsl_vector_subvector(w->work3, j, na);
 635
 636               /*
 637                * the loop below expects the solution in the first column
 638                * of sum, so set stride to 2
 639                */
 640               xv = gsl_vector_view_array_with_stride(sum, 2, na);
 641
 642               gsl_schur_solve_equation(acoef,
 643                                        &sv.matrix,
 644                                        bcoefr,
 645                                        bdiag[0],
 646                                        bdiag[1],
 647                                        &bv.vector,
 648                                        &xv.vector,
 649                                        &scale,
 650                                        &temp,
 651                                        dmin);
 652             }
 653           else
 654             {
 655               double bdat[4];
 656               gsl_matrix_const_view sv =
 657                 gsl_matrix_const_submatrix(S, j, j, na, na);
 658               gsl_vector_complex_view xv =
 659                 gsl_vector_complex_view_array(sum, na);
 660               gsl_vector_complex_view bv =
 661                 gsl_vector_complex_view_array(bdat, na);
 662               gsl_complex z;
 663
 664               bdat[0] = gsl_vector_get(w->work3, j);
 665               bdat[1] = gsl_vector_get(w->work4, j);
 666               if (na == 2)
 667                 {
 668                   bdat[2] = gsl_vector_get(w->work3, j + 1);
 669                   bdat[3] = gsl_vector_get(w->work4, j + 1);
 670                 }
 671
 672               GSL_SET_COMPLEX(&z, bcoefr, bcoefi);
 673
 674               gsl_schur_solve_equation_z(acoef,
 675                                          &sv.matrix,
 676                                          &z,
 677                                          bdiag[0],
 678                                          bdiag[1],
 679                                          &bv.vector,
 680                                          &xv.vector,
 681                                          &scale,
 682                                          &temp,
 683                                          dmin);
 684             }
 685
 686           if (scale < 1.0)
 687             {
 688               for (jr = 0; jr <= je; ++jr)
 689                 {
 690                   gsl_vector_set(w->work3, jr,
 691                     scale * gsl_vector_get(w->work3, jr));
 692                   if (nw == 2)
 693                     {
 694                       gsl_vector_set(w->work4, jr,
 695                         scale * gsl_vector_get(w->work4, jr));
 696                     }
 697                 }
 698             }
 699
 700           xmax = GSL_MAX(scale * xmax, temp);
 701
 702           for (jr = 0; jr < na; ++jr)
 703             {
 704               gsl_vector_set(w->work3, j + jr, sum[jr*na]);
 705               if (nw == 2)
 706                 gsl_vector_set(w->work4, j + jr, sum[jr*na + 1]);
 707             }
 708
 709           if (j > 0)
 710             {
 711               xscale = 1.0 / GSL_MAX(1.0, xmax);
 712               temp = acoefa * gsl_vector_get(w->work1, j) +
 713                      bcoefa * gsl_vector_get(w->work2, j);
 714               if (il2by2)
 715                 {
 716                   temp = GSL_MAX(temp,
 717                            acoefa * gsl_vector_get(w->work1, j + 1) +
 718                            bcoefa * gsl_vector_get(w->work2, j + 1));
 719                 }
 720
 721               temp = GSL_MAX(temp, GSL_MAX(acoefa, bcoefa));
 722               if (temp > bignum * xscale)
 723                 {
 724                   for (jr = 0; jr <= je; ++jr)
 725                     {
 726                       gsl_vector_set(w->work3, jr,
 727                         xscale * gsl_vector_get(w->work3, jr));
 728                       if (nw == 2)
 729                         {
 730                           gsl_vector_set(w->work4, jr,
 731                             xscale * gsl_vector_get(w->work4, jr));
 732                         }
 733                     }
 734                   xmax *= xscale;
 735                 }
 736
 737               for (ja = 0; ja < na; ++ja)
 738                 {
 739                   if (complex_pair)
 740                     {
 741                       creala = acoef * gsl_vector_get(w->work3, j + ja);
 742                       cimaga = acoef * gsl_vector_get(w->work4, j + ja);
 743                       crealb = bcoefr * gsl_vector_get(w->work3, j + ja) -
 744                                bcoefi * gsl_vector_get(w->work4, j + ja);
 745                       cimagb = bcoefi * gsl_vector_get(w->work3, j + ja) +
 746                                bcoefr * gsl_vector_get(w->work4, j + ja);
 747                       for (jr = 0; jr <= j - 1; ++jr)
 748                         {
 749                           gsl_vector_set(w->work3, jr,
 750                             gsl_vector_get(w->work3, jr) -
 751                             creala * gsl_matrix_get(S, jr, j + ja) +
 752                             crealb * gsl_matrix_get(T, jr, j + ja));
 753                           gsl_vector_set(w->work4, jr,
 754                             gsl_vector_get(w->work4, jr) -
 755                             cimaga * gsl_matrix_get(S, jr, j + ja) +
 756                             cimagb * gsl_matrix_get(T, jr, j + ja));
 757                         }
 758                     }
 759                   else
 760                     {
 761                       creala = acoef * gsl_vector_get(w->work3, j + ja);
 762                       crealb = bcoefr * gsl_vector_get(w->work3, j + ja);
 763                       for (jr = 0; jr <= j - 1; ++jr)
 764                         {
 765                           gsl_vector_set(w->work3, jr,
 766                             gsl_vector_get(w->work3, jr) -
 767                             creala * gsl_matrix_get(S, jr, j + ja) +
 768                             crealb * gsl_matrix_get(T, jr, j + ja));
 769                         }
 770                     } /* if (!complex_pair) */
 771                 } /* for (ja = 0; ja < na; ++ja) */
 772             } /* if (j > 0) */
 773
 774           il2by2 = 0;
 775         } /* for (i = 0; i < je - nw; ++i) */
 776
 777       for (jr = 0; jr < N; ++jr)
 778         {
 779           gsl_vector_set(w->work5, jr,
 780             gsl_vector_get(w->work3, 0) * gsl_matrix_get(Z, jr, 0));
 781           if (nw == 2)
 782             {
 783               gsl_vector_set(w->work6, jr,
 784                 gsl_vector_get(w->work4, 0) * gsl_matrix_get(Z, jr, 0));
 785             }
 786         }
 787
 788       for (jc = 1; jc <= je; ++jc)
 789         {
 790           for (jr = 0; jr < N; ++jr)
 791             {
 792               gsl_vector_set(w->work5, jr,
 793                 gsl_vector_get(w->work5, jr) +
 794                 gsl_vector_get(w->work3, jc) * gsl_matrix_get(Z, jr, jc));
 795               if (nw == 2)
 796                 {
 797                   gsl_vector_set(w->work6, jr,
 798                     gsl_vector_get(w->work6, jr) +
 799                     gsl_vector_get(w->work4, jc) * gsl_matrix_get(Z, jr, jc));
 800                 }
 801             }
 802         }
 803
 804       /* store the eigenvector */
 805
 806       if (complex_pair)
 807         {
 808           ecol = gsl_matrix_complex_column(evec, je - 1);
 809           re = gsl_vector_complex_real(&ecol.vector);
 810           im = gsl_vector_complex_imag(&ecol.vector);
 811
 812           ecol = gsl_matrix_complex_column(evec, je);
 813           re2 = gsl_vector_complex_real(&ecol.vector);
 814           im2 = gsl_vector_complex_imag(&ecol.vector);
 815         }
 816       else
 817         {
 818           ecol = gsl_matrix_complex_column(evec, je);
 819           re = gsl_vector_complex_real(&ecol.vector);
 820           im = gsl_vector_complex_imag(&ecol.vector);
 821         }
 822
 823       for (jr = 0; jr < N; ++jr)
 824         {
 825           gsl_vector_set(&re.vector, jr, gsl_vector_get(w->work5, jr));
 826           if (complex_pair)
 827             {
 828               gsl_vector_set(&im.vector, jr, gsl_vector_get(w->work6, jr));
 829               gsl_vector_set(&re2.vector, jr, gsl_vector_get(w->work5, jr));
 830               gsl_vector_set(&im2.vector, jr, -gsl_vector_get(w->work6, jr));
 831             }
 832           else
 833             {
 834               gsl_vector_set(&im.vector, jr, 0.0);
 835             }
 836         }
 837
 838       /* scale eigenvector */
 839       xmax = 0.0;
 840       if (complex_pair)
 841         {
 842           for (j = 0; j < N; ++j)
 843             {
 844               xmax = GSL_MAX(xmax,
 845                              fabs(gsl_vector_get(&re.vector, j)) +
 846                              fabs(gsl_vector_get(&im.vector, j)));
 847             }
 848         }
 849       else
 850         {
 851           for (j = 0; j < N; ++j)
 852             {
 853               xmax = GSL_MAX(xmax, fabs(gsl_vector_get(&re.vector, j)));
 854             }
 855         }
 856
 857       if (xmax > GSL_DBL_MIN)
 858         {
 859           xscale = 1.0 / xmax;
 860           for (j = 0; j < N; ++j)
 861             {
 862               gsl_vector_set(&re.vector, j,
 863                              gsl_vector_get(&re.vector, j) * xscale);
 864               if (complex_pair)
 865                 {
 866                   gsl_vector_set(&im.vector, j,
 867                                  gsl_vector_get(&im.vector, j) * xscale);
 868                   gsl_vector_set(&re2.vector, j,
 869                                  gsl_vector_get(&re2.vector, j) * xscale);
 870                   gsl_vector_set(&im2.vector, j,
 871                                  gsl_vector_get(&im2.vector, j) * xscale);
 872                 }
 873             }
 874         }
 875     } /* for (k = 0; k < N; ++k) */
 876
 877   return GSL_SUCCESS;
 878 } /* genv_get_right_eigenvectors() */
 879
 880 /*
 881 genv_normalize_eigenvectors()
 882   Normalize eigenvectors so that their Euclidean norm is 1
 883
 884 Inputs: alpha - eigenvalue numerators
 885         evec  - eigenvectors
 886 */
 887
 888 static void
 889 genv_normalize_eigenvectors(gsl_vector_complex *alpha,
 890                             gsl_matrix_complex *evec)
 891 {
 892   const size_t N = evec->size1;
 893   size_t i;     /* looping */
 894   gsl_complex ai;
 895   gsl_vector_complex_view vi;
 896   gsl_vector_view re, im;
 897   double scale; /* scaling factor */
 898
 899   for (i = 0; i < N; ++i)
 900     {
 901       ai = gsl_vector_complex_get(alpha, i);
 902       vi = gsl_matrix_complex_column(evec, i);
 903
 904       re = gsl_vector_complex_real(&vi.vector);
 905
 906       if (GSL_IMAG(ai) == 0.0)
 907         {
 908           scale = 1.0 / gsl_blas_dnrm2(&re.vector);
 909           gsl_blas_dscal(scale, &re.vector);
 910         }
 911       else if (GSL_IMAG(ai) > 0.0)
 912         {
 913           im = gsl_vector_complex_imag(&vi.vector);
 914
 915           scale = 1.0 / gsl_hypot(gsl_blas_dnrm2(&re.vector),
 916                                   gsl_blas_dnrm2(&im.vector));
 917           gsl_blas_zdscal(scale, &vi.vector);
 918
 919           vi = gsl_matrix_complex_column(evec, i + 1);
 920           gsl_blas_zdscal(scale, &vi.vector);
 921         }
 922     }
 923 } /* genv_normalize_eigenvectors() */