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