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

   1 /* eigen/nonsymmv.c
   2  *
   3  * Copyright (C) 2006 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 <config.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23
  24 #include <gsl/gsl_complex.h>
  25 #include <gsl/gsl_complex_math.h>
  26 #include <gsl/gsl_eigen.h>
  27 #include <gsl/gsl_linalg.h>
  28 #include <gsl/gsl_math.h>
  29 #include <gsl/gsl_blas.h>
  30 #include <gsl/gsl_cblas.h>
  31 #include <gsl/gsl_vector.h>
  32 #include <gsl/gsl_vector_complex.h>
  33 #include <gsl/gsl_matrix.h>
  34
  35 /*
  36  * This module computes the eigenvalues and eigenvectors of a real
  37  * nonsymmetric matrix.
  38  *
  39  * This file contains routines based on original code from LAPACK
  40  * which is distributed under the modified BSD license. The LAPACK
  41  * routines used are DTREVC and DLALN2.
  42  */
  43
  44 #define GSL_NONSYMMV_SMLNUM (2.0 * GSL_DBL_MIN)
  45 #define GSL_NONSYMMV_BIGNUM ((1.0 - GSL_DBL_EPSILON) / GSL_NONSYMMV_SMLNUM)
  46
  47 static void nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
  48                                             gsl_vector_complex *eval,
  49                                             gsl_matrix_complex *evec,
  50                                             gsl_eigen_nonsymmv_workspace *w);
  51 static void nonsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
  52                                             gsl_matrix_complex *evec);
  53
  54 /*
  55 gsl_eigen_nonsymmv_alloc()
  56
  57 Allocate a workspace for solving the nonsymmetric eigenvalue problem.
  58 The size of this workspace is O(5n).
  59
  60 Inputs: n - size of matrices
  61
  62 Return: pointer to workspace
  63 */
  64
  65 gsl_eigen_nonsymmv_workspace *
  66 gsl_eigen_nonsymmv_alloc(const size_t n)
  67 {
  68   gsl_eigen_nonsymmv_workspace *w;
  69
  70   if (n == 0)
  71     {
  72       GSL_ERROR_NULL ("matrix dimension must be positive integer",
  73                       GSL_EINVAL);
  74     }
  75
  76   w = (gsl_eigen_nonsymmv_workspace *)
  77       calloc (1, sizeof (gsl_eigen_nonsymmv_workspace));
  78
  79   if (w == 0)
  80     {
  81       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
  82     }
  83
  84   w->size = n;
  85   w->Z = NULL;
  86   w->nonsymm_workspace_p = gsl_eigen_nonsymm_alloc(n);
  87
  88   if (w->nonsymm_workspace_p == 0)
  89     {
  90       gsl_eigen_nonsymmv_free(w);
  91       GSL_ERROR_NULL ("failed to allocate space for nonsymm workspace", GSL_ENOMEM);
  92     }
  93
  94   /*
  95    * set parameters to compute the full Schur form T and balance
  96    * the matrices
  97    */
  98   gsl_eigen_nonsymm_params(1, 1, w->nonsymm_workspace_p);
  99
 100   w->work = gsl_vector_alloc(n);
 101   w->work2 = gsl_vector_alloc(n);
 102   w->work3 = gsl_vector_alloc(n);
 103   if (w->work == 0 || w->work2 == 0 || w->work3 == 0)
 104     {
 105       gsl_eigen_nonsymmv_free(w);
 106       GSL_ERROR_NULL ("failed to allocate space for nonsymmv additional workspace", GSL_ENOMEM);
 107     }
 108
 109   return (w);
 110 } /* gsl_eigen_nonsymmv_alloc() */
 111
 112 /*
 113 gsl_eigen_nonsymmv_free()
 114   Free workspace w
 115 */
 116
 117 void
 118 gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w)
 119 {
 120   if (w->nonsymm_workspace_p)
 121     gsl_eigen_nonsymm_free(w->nonsymm_workspace_p);
 122
 123   if (w->work)
 124     gsl_vector_free(w->work);
 125
 126   if (w->work2)
 127     gsl_vector_free(w->work2);
 128
 129   if (w->work3)
 130     gsl_vector_free(w->work3);
 131
 132   free(w);
 133 } /* gsl_eigen_nonsymmv_free() */
 134
 135 /*
 136 gsl_eigen_nonsymmv()
 137
 138 Solve the nonsymmetric eigensystem problem
 139
 140 A x = \lambda x
 141
 142 for the eigenvalues \lambda and right eigenvectors x
 143
 144 Inputs: A    - general real matrix
 145         eval - where to store eigenvalues
 146         evec - where to store eigenvectors
 147         w    - workspace
 148
 149 Return: success or error
 150 */
 151
 152 int
 153 gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
 154                     gsl_matrix_complex * evec,
 155                     gsl_eigen_nonsymmv_workspace * w)
 156 {
 157   const size_t N = A->size1;
 158
 159   /* check matrix and vector sizes */
 160
 161   if (N != A->size2)
 162     {
 163       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 164     }
 165   else if (eval->size != N)
 166     {
 167       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 168     }
 169   else if (evec->size1 != evec->size2)
 170     {
 171       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
 172     }
 173   else if (evec->size1 != N)
 174     {
 175       GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
 176     }
 177   else
 178     {
 179       int s;
 180       gsl_matrix Z;
 181
 182       /*
 183        * We need a place to store the Schur vectors, so we will
 184        * treat evec as a real matrix and store them in the left
 185        * half - the factor of 2 in the tda corresponds to the
 186        * complex multiplicity
 187        */
 188       Z.size1 = N;
 189       Z.size2 = N;
 190       Z.tda = 2 * N;
 191       Z.data = evec->data;
 192       Z.block = 0;
 193       Z.owner = 0;
 194
 195       /* compute eigenvalues, Schur form, and Schur vectors */
 196       s = gsl_eigen_nonsymm_Z(A, eval, &Z, w->nonsymm_workspace_p);
 197
 198       if (w->Z)
 199         {
 200           /*
 201            * save the Schur vectors in user supplied matrix, since
 202            * they will be destroyed when computing eigenvectors
 203            */
 204           gsl_matrix_memcpy(w->Z, &Z);
 205         }
 206
 207       /* only compute eigenvectors if we found all eigenvalues */
 208       if (s == GSL_SUCCESS)
 209         {
 210           /* compute eigenvectors */
 211           nonsymmv_get_right_eigenvectors(A, &Z, eval, evec, w);
 212
 213           /* normalize so that Euclidean norm is 1 */
 214           nonsymmv_normalize_eigenvectors(eval, evec);
 215         }
 216
 217       return s;
 218     }
 219 } /* gsl_eigen_nonsymmv() */
 220
 221 /*
 222 gsl_eigen_nonsymmv_Z()
 223   Compute eigenvalues and eigenvectors of a real nonsymmetric matrix
 224 and also save the Schur vectors. See comments in gsl_eigen_nonsymm_Z
 225 for more information.
 226
 227 Inputs: A    - real nonsymmetric matrix
 228         eval - where to store eigenvalues
 229         evec - where to store eigenvectors
 230         Z    - where to store Schur vectors
 231         w    - nonsymmv workspace
 232
 233 Return: success or error
 234 */
 235
 236 int
 237 gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
 238                       gsl_matrix_complex * evec, gsl_matrix * Z,
 239                       gsl_eigen_nonsymmv_workspace * w)
 240 {
 241   /* check matrix and vector sizes */
 242
 243   if (A->size1 != A->size2)
 244     {
 245       GSL_ERROR ("matrix must be square to compute eigenvalues/eigenvectors", GSL_ENOTSQR);
 246     }
 247   else if (eval->size != A->size1)
 248     {
 249       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 250     }
 251   else if (evec->size1 != evec->size2)
 252     {
 253       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
 254     }
 255   else if (evec->size1 != A->size1)
 256     {
 257       GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
 258     }
 259   else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
 260     {
 261       GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
 262     }
 263   else
 264     {
 265       int s;
 266
 267       w->Z = Z;
 268
 269       s = gsl_eigen_nonsymmv(A, eval, evec, w);
 270
 271       w->Z = NULL;
 272
 273       return s;
 274     }
 275 } /* gsl_eigen_nonsymmv_Z() */
 276
 277 /********************************************
 278  *           INTERNAL ROUTINES              *
 279  ********************************************/
 280
 281 /*
 282 nonsymmv_get_right_eigenvectors()
 283   Compute the right eigenvectors of the Schur form T and then
 284 backtransform them using the Schur vectors to get right eigenvectors of
 285 the original matrix.
 286
 287 Inputs: T    - Schur form
 288         Z    - Schur vectors
 289         eval - where to store eigenvalues (to ensure that the
 290                correct eigenvalue is stored in the same position
 291                as the eigenvectors)
 292         evec - where to store eigenvectors
 293         w    - nonsymmv workspace
 294
 295 Return: none
 296
 297 Notes: 1) based on LAPACK routine DTREVC - the algorithm used is
 298           backsubstitution on the upper quasi triangular system T
 299           followed by backtransformation by Z to get vectors of the
 300           original matrix.
 301
 302        2) The Schur vectors in Z are destroyed and replaced with
 303           eigenvectors stored with the same storage scheme as DTREVC.
 304           The eigenvectors are also stored in 'evec'
 305
 306        3) The matrix T is unchanged on output
 307
 308        4) Each eigenvector is normalized so that the element of
 309           largest magnitude has magnitude 1; here the magnitude of
 310           a complex number (x,y) is taken to be |x| + |y|
 311 */
 312
 313 static void
 314 nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
 315                                 gsl_vector_complex *eval,
 316                                 gsl_matrix_complex *evec,
 317                                 gsl_eigen_nonsymmv_workspace *w)
 318 {
 319   const size_t N = T->size1;
 320   const double smlnum = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
 321   const double bignum = (1.0 - GSL_DBL_EPSILON) / smlnum;
 322   int i;              /* looping */
 323   size_t iu,          /* looping */
 324          ju,
 325          ii;
 326   gsl_complex lambda; /* current eigenvalue */
 327   double lambda_re,   /* Re(lambda) */
 328          lambda_im;   /* Im(lambda) */
 329   gsl_matrix_view Tv, /* temporary views */
 330                   Zv;
 331   gsl_vector_view y,  /* temporary views */
 332                   y2,
 333                   ev,
 334                   ev2;
 335   double dat[4],      /* scratch arrays */
 336          dat_X[4];
 337   double scale;       /* scale factor */
 338   double xnorm;       /* |X| */
 339   gsl_vector_complex_view ecol, /* column of evec */
 340                           ecol2;
 341   int complex_pair;   /* complex eigenvalue pair? */
 342   double smin;
 343
 344   /*
 345    * Compute 1-norm of each column of upper triangular part of T
 346    * to control overflow in triangular solver
 347    */
 348
 349   gsl_vector_set(w->work3, 0, 0.0);
 350   for (ju = 1; ju < N; ++ju)
 351     {
 352       gsl_vector_set(w->work3, ju, 0.0);
 353       for (iu = 0; iu < ju; ++iu)
 354         {
 355           gsl_vector_set(w->work3, ju,
 356                          gsl_vector_get(w->work3, ju) +
 357                          fabs(gsl_matrix_get(T, iu, ju)));
 358         }
 359     }
 360
 361   for (i = (int) N - 1; i >= 0; --i)
 362     {
 363       iu = (size_t) i;
 364
 365       /* get current eigenvalue and store it in lambda */
 366       lambda_re = gsl_matrix_get(T, iu, iu);
 367
 368       if (iu != 0 && gsl_matrix_get(T, iu, iu - 1) != 0.0)
 369         {
 370           lambda_im = sqrt(fabs(gsl_matrix_get(T, iu, iu - 1))) *
 371                       sqrt(fabs(gsl_matrix_get(T, iu - 1, iu)));
 372         }
 373       else
 374         {
 375           lambda_im = 0.0;
 376         }
 377
 378       GSL_SET_COMPLEX(&lambda, lambda_re, lambda_im);
 379
 380       smin = GSL_MAX(GSL_DBL_EPSILON * (fabs(lambda_re) + fabs(lambda_im)),
 381                      smlnum);
 382       smin = GSL_MAX(smin, GSL_NONSYMMV_SMLNUM);
 383
 384       if (lambda_im == 0.0)
 385         {
 386           int k, l;
 387           gsl_vector_view bv, xv;
 388
 389           /* real eigenvector */
 390
 391           /*
 392            * The ordering of eigenvalues in 'eval' is arbitrary and
 393            * does not necessarily follow the Schur form T, so store
 394            * lambda in the right slot in eval to ensure it corresponds
 395            * to the eigenvector we are about to compute
 396            */
 397           gsl_vector_complex_set(eval, iu, lambda);
 398
 399           /*
 400            * We need to solve the system:
 401            *
 402            * (T(1:iu-1, 1:iu-1) - lambda*I)*X = -T(1:iu-1,iu)
 403            */
 404
 405           /* construct right hand side */
 406           for (k = 0; k < i; ++k)
 407             {
 408               gsl_vector_set(w->work,
 409                              (size_t) k,
 410                              -gsl_matrix_get(T, (size_t) k, iu));
 411             }
 412
 413           gsl_vector_set(w->work, iu, 1.0);
 414
 415           for (l = i - 1; l >= 0; --l)
 416             {
 417               size_t lu = (size_t) l;
 418
 419               if (lu == 0)
 420                 complex_pair = 0;
 421               else
 422                 complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;
 423
 424               if (!complex_pair)
 425                 {
 426                   double x;
 427
 428                   /*
 429                    * 1-by-1 diagonal block - solve the system:
 430                    *
 431                    * (T_{ll} - lambda)*x = -T_{l(iu)}
 432                    */
 433
 434                   Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
 435                   bv = gsl_vector_view_array(dat, 1);
 436                   gsl_vector_set(&bv.vector, 0,
 437                                  gsl_vector_get(w->work, lu));
 438                   xv = gsl_vector_view_array(dat_X, 1);
 439
 440                   gsl_schur_solve_equation(1.0,
 441                                            &Tv.matrix,
 442                                            lambda_re,
 443                                            1.0,
 444                                            1.0,
 445                                            &bv.vector,
 446                                            &xv.vector,
 447                                            &scale,
 448                                            &xnorm,
 449                                            smin);
 450
 451                   /* scale x to avoid overflow */
 452                   x = gsl_vector_get(&xv.vector, 0);
 453                   if (xnorm > 1.0)
 454                     {
 455                       if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
 456                         {
 457                           x /= xnorm;
 458                           scale /= xnorm;
 459                         }
 460                     }
 461
 462                   if (scale != 1.0)
 463                     {
 464                       gsl_vector_view wv;
 465
 466                       wv = gsl_vector_subvector(w->work, 0, iu + 1);
 467                       gsl_blas_dscal(scale, &wv.vector);
 468                     }
 469
 470                   gsl_vector_set(w->work, lu, x);
 471
 472                   if (lu > 0)
 473                     {
 474                       gsl_vector_view v1, v2;
 475
 476                       /* update right hand side */
 477
 478                       v1 = gsl_matrix_subcolumn(T, lu, 0, lu);
 479                       v2 = gsl_vector_subvector(w->work, 0, lu);
 480                       gsl_blas_daxpy(-x, &v1.vector, &v2.vector);
 481                     } /* if (l > 0) */
 482                 } /* if (!complex_pair) */
 483               else
 484                 {
 485                   double x11, x21;
 486
 487                   /*
 488                    * 2-by-2 diagonal block
 489                    */
 490
 491                   Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
 492                   bv = gsl_vector_view_array(dat, 2);
 493                   gsl_vector_set(&bv.vector, 0,
 494                                  gsl_vector_get(w->work, lu - 1));
 495                   gsl_vector_set(&bv.vector, 1,
 496                                  gsl_vector_get(w->work, lu));
 497                   xv = gsl_vector_view_array(dat_X, 2);
 498
 499                   gsl_schur_solve_equation(1.0,
 500                                            &Tv.matrix,
 501                                            lambda_re,
 502                                            1.0,
 503                                            1.0,
 504                                            &bv.vector,
 505                                            &xv.vector,
 506                                            &scale,
 507                                            &xnorm,
 508                                            smin);
 509
 510                   /* scale X(1,1) and X(2,1) to avoid overflow */
 511                   x11 = gsl_vector_get(&xv.vector, 0);
 512                   x21 = gsl_vector_get(&xv.vector, 1);
 513
 514                   if (xnorm > 1.0)
 515                     {
 516                       double beta;
 517
 518                       beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
 519                                      gsl_vector_get(w->work3, lu));
 520                       if (beta > bignum / xnorm)
 521                         {
 522                           x11 /= xnorm;
 523                           x21 /= xnorm;
 524                           scale /= xnorm;
 525                         }
 526                     }
 527
 528                   /* scale if necessary */
 529                   if (scale != 1.0)
 530                     {
 531                       gsl_vector_view wv;
 532
 533                       wv = gsl_vector_subvector(w->work, 0, iu + 1);
 534                       gsl_blas_dscal(scale, &wv.vector);
 535                     }
 536
 537                   gsl_vector_set(w->work, lu - 1, x11);
 538                   gsl_vector_set(w->work, lu, x21);
 539
 540                   /* update right hand side */
 541                   if (lu > 1)
 542                     {
 543                       gsl_vector_view v1, v2;
 544
 545                       v1 = gsl_matrix_subcolumn(T, lu - 1, 0, lu - 1);
 546                       v2 = gsl_vector_subvector(w->work, 0, lu - 1);
 547                       gsl_blas_daxpy(-x11, &v1.vector, &v2.vector);
 548
 549                       v1 = gsl_matrix_subcolumn(T, lu, 0, lu - 1);
 550                       gsl_blas_daxpy(-x21, &v1.vector, &v2.vector);
 551                     }
 552
 553                   --l;
 554                 } /* if (complex_pair) */
 555             } /* for (l = i - 1; l >= 0; --l) */
 556
 557           /*
 558            * At this point, w->work is an eigenvector of the
 559            * Schur form T. To get an eigenvector of the original
 560            * matrix, we multiply on the left by Z, the matrix of
 561            * Schur vectors
 562            */
 563
 564           ecol = gsl_matrix_complex_column(evec, iu);
 565           y = gsl_matrix_column(Z, iu);
 566
 567           if (iu > 0)
 568             {
 569               gsl_vector_view x;
 570
 571               Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu);
 572
 573               x = gsl_vector_subvector(w->work, 0, iu);
 574
 575               /* compute Z * w->work and store it in Z(:,iu) */
 576               gsl_blas_dgemv(CblasNoTrans,
 577                              1.0,
 578                              &Zv.matrix,
 579                              &x.vector,
 580                              gsl_vector_get(w->work, iu),
 581                              &y.vector);
 582             } /* if (iu > 0) */
 583
 584           /* store eigenvector into evec */
 585
 586           ev = gsl_vector_complex_real(&ecol.vector);
 587           ev2 = gsl_vector_complex_imag(&ecol.vector);
 588
 589           scale = 0.0;
 590           for (ii = 0; ii < N; ++ii)
 591             {
 592               double a = gsl_vector_get(&y.vector, ii);
 593
 594               /* store real part of eigenvector */
 595               gsl_vector_set(&ev.vector, ii, a);
 596
 597               /* set imaginary part to 0 */
 598               gsl_vector_set(&ev2.vector, ii, 0.0);
 599
 600               if (fabs(a) > scale)
 601                 scale = fabs(a);
 602             }
 603
 604           if (scale != 0.0)
 605             scale = 1.0 / scale;
 606
 607           /* scale by magnitude of largest element */
 608           gsl_blas_dscal(scale, &ev.vector);
 609         } /* if (GSL_IMAG(lambda) == 0.0) */
 610       else
 611         {
 612           gsl_vector_complex_view bv, xv;
 613           size_t k;
 614           int l;
 615           gsl_complex lambda2;
 616
 617           /* complex eigenvector */
 618
 619           /*
 620            * Store the complex conjugate eigenvalues in the right
 621            * slots in eval
 622            */
 623           GSL_SET_REAL(&lambda2, GSL_REAL(lambda));
 624           GSL_SET_IMAG(&lambda2, -GSL_IMAG(lambda));
 625           gsl_vector_complex_set(eval, iu - 1, lambda);
 626           gsl_vector_complex_set(eval, iu, lambda2);
 627
 628           /*
 629            * First solve:
 630            *
 631            * [ T(i:i+1,i:i+1) - lambda*I ] * X = 0
 632            */
 633
 634           if (fabs(gsl_matrix_get(T, iu - 1, iu)) >=
 635               fabs(gsl_matrix_get(T, iu, iu - 1)))
 636             {
 637               gsl_vector_set(w->work, iu - 1, 1.0);
 638               gsl_vector_set(w->work2, iu,
 639                              lambda_im / gsl_matrix_get(T, iu - 1, iu));
 640             }
 641           else
 642             {
 643               gsl_vector_set(w->work, iu - 1,
 644                              -lambda_im / gsl_matrix_get(T, iu, iu - 1));
 645               gsl_vector_set(w->work2, iu, 1.0);
 646             }
 647           gsl_vector_set(w->work, iu, 0.0);
 648           gsl_vector_set(w->work2, iu - 1, 0.0);
 649
 650           /* construct right hand side */
 651           for (k = 0; k < iu - 1; ++k)
 652             {
 653               gsl_vector_set(w->work, k,
 654                              -gsl_vector_get(w->work, iu - 1) *
 655                              gsl_matrix_get(T, k, iu - 1));
 656               gsl_vector_set(w->work2, k,
 657                              -gsl_vector_get(w->work2, iu) *
 658                              gsl_matrix_get(T, k, iu));
 659             }
 660
 661           /*
 662            * We must solve the upper quasi-triangular system:
 663            *
 664            * [ T(1:i-2,1:i-2) - lambda*I ] * X = s*(work + i*work2)
 665            */
 666
 667           for (l = i - 2; l >= 0; --l)
 668             {
 669               size_t lu = (size_t) l;
 670
 671               if (lu == 0)
 672                 complex_pair = 0;
 673               else
 674                 complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;
 675
 676               if (!complex_pair)
 677                 {
 678                   gsl_complex bval;
 679                   gsl_complex x;
 680
 681                   /*
 682                    * 1-by-1 diagonal block - solve the system:
 683                    *
 684                    * (T_{ll} - lambda)*x = work + i*work2
 685                    */
 686
 687                   Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
 688                   bv = gsl_vector_complex_view_array(dat, 1);
 689                   xv = gsl_vector_complex_view_array(dat_X, 1);
 690
 691                   GSL_SET_COMPLEX(&bval,
 692                                   gsl_vector_get(w->work, lu),
 693                                   gsl_vector_get(w->work2, lu));
 694                   gsl_vector_complex_set(&bv.vector, 0, bval);
 695
 696                   gsl_schur_solve_equation_z(1.0,
 697                                              &Tv.matrix,
 698                                              &lambda,
 699                                              1.0,
 700                                              1.0,
 701                                              &bv.vector,
 702                                              &xv.vector,
 703                                              &scale,
 704                                              &xnorm,
 705                                              smin);
 706
 707                   if (xnorm > 1.0)
 708                     {
 709                       if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
 710                         {
 711                           gsl_blas_zdscal(1.0/xnorm, &xv.vector);
 712                           scale /= xnorm;
 713                         }
 714                     }
 715
 716                   /* scale if necessary */
 717                   if (scale != 1.0)
 718                     {
 719                       gsl_vector_view wv;
 720
 721                       wv = gsl_vector_subvector(w->work, 0, iu + 1);
 722                       gsl_blas_dscal(scale, &wv.vector);
 723                       wv = gsl_vector_subvector(w->work2, 0, iu + 1);
 724                       gsl_blas_dscal(scale, &wv.vector);
 725                     }
 726
 727                   x = gsl_vector_complex_get(&xv.vector, 0);
 728                   gsl_vector_set(w->work, lu, GSL_REAL(x));
 729                   gsl_vector_set(w->work2, lu, GSL_IMAG(x));
 730
 731                   /* update the right hand side */
 732                   if (lu > 0)
 733                     {
 734                       gsl_vector_view v1, v2;
 735
 736                       v1 = gsl_matrix_subcolumn(T, lu, 0, lu);
 737                       v2 = gsl_vector_subvector(w->work, 0, lu);
 738                       gsl_blas_daxpy(-GSL_REAL(x), &v1.vector, &v2.vector);
 739
 740                       v2 = gsl_vector_subvector(w->work2, 0, lu);
 741                       gsl_blas_daxpy(-GSL_IMAG(x), &v1.vector, &v2.vector);
 742                     } /* if (lu > 0) */
 743                 } /* if (!complex_pair) */
 744               else
 745                 {
 746                   gsl_complex b1, b2, x1, x2;
 747
 748                   /*
 749                    * 2-by-2 diagonal block - solve the system
 750                    */
 751
 752                   Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
 753                   bv = gsl_vector_complex_view_array(dat, 2);
 754                   xv = gsl_vector_complex_view_array(dat_X, 2);
 755
 756                   GSL_SET_COMPLEX(&b1,
 757                                   gsl_vector_get(w->work, lu - 1),
 758                                   gsl_vector_get(w->work2, lu - 1));
 759                   GSL_SET_COMPLEX(&b2,
 760                                   gsl_vector_get(w->work, lu),
 761                                   gsl_vector_get(w->work2, lu));
 762                   gsl_vector_complex_set(&bv.vector, 0, b1);
 763                   gsl_vector_complex_set(&bv.vector, 1, b2);
 764
 765                   gsl_schur_solve_equation_z(1.0,
 766                                              &Tv.matrix,
 767                                              &lambda,
 768                                              1.0,
 769                                              1.0,
 770                                              &bv.vector,
 771                                              &xv.vector,
 772                                              &scale,
 773                                              &xnorm,
 774                                              smin);
 775
 776                   x1 = gsl_vector_complex_get(&xv.vector, 0);
 777                   x2 = gsl_vector_complex_get(&xv.vector, 1);
 778
 779                   if (xnorm > 1.0)
 780                     {
 781                       double beta;
 782
 783                       beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
 784                                      gsl_vector_get(w->work3, lu));
 785                       if (beta > bignum / xnorm)
 786                         {
 787                           gsl_blas_zdscal(1.0/xnorm, &xv.vector);
 788                           scale /= xnorm;
 789                         }
 790                     }
 791
 792                   /* scale if necessary */
 793                   if (scale != 1.0)
 794                     {
 795                       gsl_vector_view wv;
 796
 797                       wv = gsl_vector_subvector(w->work, 0, iu + 1);
 798                       gsl_blas_dscal(scale, &wv.vector);
 799                       wv = gsl_vector_subvector(w->work2, 0, iu + 1);
 800                       gsl_blas_dscal(scale, &wv.vector);
 801                     }
 802                   gsl_vector_set(w->work, lu - 1, GSL_REAL(x1));
 803                   gsl_vector_set(w->work, lu, GSL_REAL(x2));
 804                   gsl_vector_set(w->work2, lu - 1, GSL_IMAG(x1));
 805                   gsl_vector_set(w->work2, lu, GSL_IMAG(x2));
 806
 807                   /* update right hand side */
 808                   if (lu > 1)
 809                     {
 810                       gsl_vector_view v1, v2, v3, v4;
 811
 812                       v1 = gsl_matrix_subcolumn(T, lu - 1, 0, lu - 1);
 813                       v4 = gsl_matrix_subcolumn(T, lu, 0, lu - 1);
 814                       v2 = gsl_vector_subvector(w->work, 0, lu - 1);
 815                       v3 = gsl_vector_subvector(w->work2, 0, lu - 1);
 816
 817                       gsl_blas_daxpy(-GSL_REAL(x1), &v1.vector, &v2.vector);
 818                       gsl_blas_daxpy(-GSL_REAL(x2), &v4.vector, &v2.vector);
 819                       gsl_blas_daxpy(-GSL_IMAG(x1), &v1.vector, &v3.vector);
 820                       gsl_blas_daxpy(-GSL_IMAG(x2), &v4.vector, &v3.vector);
 821                     } /* if (lu > 1) */
 822
 823                   --l;
 824                 } /* if (complex_pair) */
 825             } /* for (l = i - 2; l >= 0; --l) */
 826
 827           /*
 828            * At this point, work + i*work2 is an eigenvector
 829            * of T - backtransform to get an eigenvector of the
 830            * original matrix
 831            */
 832
 833           y = gsl_matrix_column(Z, iu - 1);
 834           y2 = gsl_matrix_column(Z, iu);
 835
 836           if (iu > 1)
 837             {
 838               gsl_vector_view x;
 839
 840               /* compute real part of eigenvectors */
 841
 842               Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu - 1);
 843               x = gsl_vector_subvector(w->work, 0, iu - 1);
 844
 845               gsl_blas_dgemv(CblasNoTrans,
 846                              1.0,
 847                              &Zv.matrix,
 848                              &x.vector,
 849                              gsl_vector_get(w->work, iu - 1),
 850                              &y.vector);
 851
 852
 853               /* now compute the imaginary part */
 854               x = gsl_vector_subvector(w->work2, 0, iu - 1);
 855
 856               gsl_blas_dgemv(CblasNoTrans,
 857                              1.0,
 858                              &Zv.matrix,
 859                              &x.vector,
 860                              gsl_vector_get(w->work2, iu),
 861                              &y2.vector);
 862             }
 863           else
 864             {
 865               gsl_blas_dscal(gsl_vector_get(w->work, iu - 1), &y.vector);
 866               gsl_blas_dscal(gsl_vector_get(w->work2, iu), &y2.vector);
 867             }
 868
 869           /*
 870            * Now store the eigenvectors into evec - the real parts
 871            * are Z(:,iu - 1) and the imaginary parts are
 872            * +/- Z(:,iu)
 873            */
 874
 875           /* get views of the two eigenvector slots */
 876           ecol = gsl_matrix_complex_column(evec, iu - 1);
 877           ecol2 = gsl_matrix_complex_column(evec, iu);
 878
 879           /*
 880            * save imaginary part first as it may get overwritten
 881            * when copying the real part due to our storage scheme
 882            * in Z/evec
 883            */
 884           ev = gsl_vector_complex_imag(&ecol.vector);
 885           ev2 = gsl_vector_complex_imag(&ecol2.vector);
 886           scale = 0.0;
 887           for (ii = 0; ii < N; ++ii)
 888             {
 889               double a = gsl_vector_get(&y2.vector, ii);
 890
 891               scale = GSL_MAX(scale,
 892                               fabs(a) + fabs(gsl_vector_get(&y.vector, ii)));
 893
 894               gsl_vector_set(&ev.vector, ii, a);
 895               gsl_vector_set(&ev2.vector, ii, -a);
 896             }
 897
 898           /* now save the real part */
 899           ev = gsl_vector_complex_real(&ecol.vector);
 900           ev2 = gsl_vector_complex_real(&ecol2.vector);
 901           for (ii = 0; ii < N; ++ii)
 902             {
 903               double a = gsl_vector_get(&y.vector, ii);
 904
 905               gsl_vector_set(&ev.vector, ii, a);
 906               gsl_vector_set(&ev2.vector, ii, a);
 907             }
 908
 909           if (scale != 0.0)
 910             scale = 1.0 / scale;
 911
 912           /* scale by largest element magnitude */
 913
 914           gsl_blas_zdscal(scale, &ecol.vector);
 915           gsl_blas_zdscal(scale, &ecol2.vector);
 916
 917           /*
 918            * decrement i since we took care of two eigenvalues at
 919            * the same time
 920            */
 921           --i;
 922         } /* if (GSL_IMAG(lambda) != 0.0) */
 923     } /* for (i = (int) N - 1; i >= 0; --i) */
 924 } /* nonsymmv_get_right_eigenvectors() */
 925
 926 /*
 927 nonsymmv_normalize_eigenvectors()
 928   Normalize eigenvectors so that their Euclidean norm is 1
 929
 930 Inputs: eval - eigenvalues
 931         evec - eigenvectors
 932 */
 933
 934 static void
 935 nonsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
 936                                 gsl_matrix_complex *evec)
 937 {
 938   const size_t N = evec->size1;
 939   size_t i;     /* looping */
 940   gsl_complex ei;
 941   gsl_vector_complex_view vi;
 942   gsl_vector_view re, im;
 943   double scale; /* scaling factor */
 944
 945   for (i = 0; i < N; ++i)
 946     {
 947       ei = gsl_vector_complex_get(eval, i);
 948       vi = gsl_matrix_complex_column(evec, i);
 949
 950       re = gsl_vector_complex_real(&vi.vector);
 951
 952       if (GSL_IMAG(ei) == 0.0)
 953         {
 954           scale = 1.0 / gsl_blas_dnrm2(&re.vector);
 955           gsl_blas_dscal(scale, &re.vector);
 956         }
 957       else if (GSL_IMAG(ei) > 0.0)
 958         {
 959           im = gsl_vector_complex_imag(&vi.vector);
 960
 961           scale = 1.0 / gsl_hypot(gsl_blas_dnrm2(&re.vector),
 962                                   gsl_blas_dnrm2(&im.vector));
 963           gsl_blas_zdscal(scale, &vi.vector);
 964
 965           vi = gsl_matrix_complex_column(evec, i + 1);
 966           gsl_blas_zdscal(scale, &vi.vector);
 967         }
 968     }
 969 } /* nonsymmv_normalize_eigenvectors() */