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