* unsymm eigenvector patch
@ 2006-09-21 18:43 Patrick Alken
2006-09-22 13:49 ` Brian Gough
0 siblings, 1 reply; 2+ messages in thread
From: Patrick Alken @ 2006-09-21 18:43 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 693 bytes --]
Hello,
I've attached a patch which implements eigenvector support for
real nonsymmetric matrices. The algorithm used is backsubstitution
to find eigenvectors of the Schur form, followed by backtransformation
by the Schur vectors. I followed the lapack code closely (dtrevc),
so this implementation handles degeneracies very well.
The patch includes:
* unsymmv.c - all the eigenvector stuff
* updated eigen.texi to give documentation on all new functions plus
an example
* updated eigen/sort.c to add gsl_eigen_unsymmv_sort()
I also attached a new test program which no longer needs lapack
to verify the results. I haven't added anything to eigen/test.c yet.
Patrick Alken
[-- Attachment #2: gsl.unsymmv.patch --]
[-- Type: text/plain, Size: 63744 bytes --]
? doc/examples/a.out
Index: doc/eigen.texi
===================================================================
RCS file: /cvs/gsl/gsl/doc/eigen.texi,v
retrieving revision 1.25
diff -u -r1.25 eigen.texi
--- doc/eigen.texi 18 Sep 2006 21:26:59 -0000 1.25
+++ doc/eigen.texi 21 Sep 2006 18:26:27 -0000
@@ -235,6 +235,31 @@
computes the Schur vectors and stores them into @var{Z}.
@end deftypefun
+@deftypefun {gsl_eigen_unsymmv_workspace *} gsl_eigen_unsymmv_alloc (const size_t @var{n})
+This function allocates a workspace for computing eigenvalues and
+eigenvectors of @var{n}-by-@var{n} real unsymmetric matrices. The
+size of the workspace is @math{O(5n)}.
+@end deftypefun
+
+@deftypefun void gsl_eigen_unsymmv_free (gsl_eigen_unsymmv_workspace * @var{w})
+This function frees the memory associated with the workspace @var{w}.
+@end deftypefun
+
+@deftypefun int gsl_eigen_unsymmv (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_unsymmv_workspace * @var{w})
+This function computes eigenvalues and right eigenvectors of the
+@var{n}-by-@var{n} real unsymmetric matrix @var{A}. It first calls
+@code{gsl_eigen_unsymm} to compute the eigenvalues, Schur form @math{T}, and
+Schur vectors. Then it finds eigenvectors of @math{T} and backtransforms
+them using the Schur vectors. The Schur vectors are destroyed in the
+process, but can be saved by using @code{gsl_eigen_unsymmv_Z}. The
+computed eigenvectors are normalized to have Euclidean norm 1.
+@end deftypefun
+
+@deftypefun int gsl_eigen_unsymmv_Z (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_matrix * @var{Z}, gsl_eigen_unsymmv_workspace * @var{w})
+This function is identical to @code{gsl_eigen_unsymmv} except it also saves
+the Schur vectors into @var{Z}.
+@end deftypefun
+
@node Sorting Eigenvalues and Eigenvectors
@section Sorting Eigenvalues and Eigenvectors
@cindex sorting eigenvalues and eigenvectors
@@ -265,6 +290,15 @@
according to the value of the parameter @var{sort_type} as shown above.
@end deftypefun
+@deftypefun int gsl_eigen_unsymmv_sort (gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_sort_t @var{sort_type})
+This function simultaneously sorts the eigenvalues stored in the vector
+@var{eval} and the corresponding complex eigenvectors stored in the
+columns of the matrix @var{evec} into ascending or descending order
+according to the value of the parameter @var{sort_type} as shown above.
+Only GSL_EIGEN_SORT_ABS_ASC and GSL_EIGEN_SORT_ABS_DESC are supported
+due to the eigenvalues being complex.
+@end deftypefun
+
@comment @deftypefun int gsl_eigen_jacobi (gsl_matrix * @var{matrix}, gsl_vector * @var{eval}, gsl_matrix * @var{evec}, unsigned int @var{max_rot}, unsigned int * @var{nrot})
@comment This function finds the eigenvectors and eigenvalues of a real symmetric
@@ -332,6 +366,51 @@
Note that the eigenvectors can differ by a change of sign, since the
sign of an eigenvector is arbitrary.
+The following program illustrates the use of the unsymmetric
+eigensolver, by computing the eigenvalues and eigenvectors of
+the Vandermonde matrix @math{V(x;i,j) = x_i^{n - j}} with
+@math{x = (1,2,3,4)}.
+
+@example
+@verbatiminclude examples/eigen_unsymm.c
+@end example
+
+@noindent
+Here is the beginning of the output from the program,
+
+@example
+$ ./a.out
+eigenvalue = 15.4897 + 0i
+eigenvector =
+0.106363 + 0i
+0.228569 + 0i
+0.462688 + 0i
+0.849919 + 0i
+...
+@end example
+
+@noindent
+This can be compared with the corresponding output from @sc{gnu octave},
+
+@example
+octave> [v,d] = eig(vander([1 2 3 4]));
+octave> diag(d)
+ans =
+
+ 15.489740
+ -7.706289
+ 1.294224
+ -0.077675
+
+octave> v
+v =
+
+ 0.106363 0.137820 -0.129196 0.052254
+ 0.228569 0.036136 0.407586 -0.375568
+ 0.462688 -0.289537 0.376326 0.794157
+ 0.849919 -0.946503 -0.821925 -0.474902
+@end example
+
@node Eigenvalue and Eigenvector References
@section References and Further Reading
Index: doc/linalg.texi
===================================================================
RCS file: /cvs/gsl/gsl/doc/linalg.texi,v
retrieving revision 1.52
diff -u -r1.52 linalg.texi
--- doc/linalg.texi 18 Sep 2006 21:26:59 -0000 1.52
+++ doc/linalg.texi 21 Sep 2006 18:26:29 -0000
@@ -966,7 +966,7 @@
where @math{D} is a diagonal matrix whose entries are powers
of the floating point radix in order to reduce roundoff error.
-@deftypefun int gsl_linalg_balance (gsl_matrix * @var{A}, gsl_vector * @var{D})
+@deftypefun int gsl_linalg_balance_matrix (gsl_matrix * @var{A}, gsl_vector * @var{D})
This function replaces the matrix @var{A} with its balanced counterpart
and stores the diagonal elements of the similarity transformation
into the vector @var{D}.
Index: doc/examples/eigen_unsymm.c
===================================================================
RCS file: doc/examples/eigen_unsymm.c
diff -N doc/examples/eigen_unsymm.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ doc/examples/eigen_unsymm.c 21 Sep 2006 18:26:29 -0000
@@ -0,0 +1,51 @@
+#include <stdio.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_eigen.h>
+
+int
+main (void)
+{
+ double data[] = { 1.0, 1.0, 1.0, 1.0,
+ 8.0, 4.0, 2.0, 1.0,
+ 27.0, 9.0, 3.0, 1.0,
+ 64.0, 16.0, 4.0, 1.0 };
+
+ gsl_matrix_view m
+ = gsl_matrix_view_array (data, 4, 4);
+
+ gsl_vector_complex *eval = gsl_vector_complex_alloc (4);
+ gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);
+
+ gsl_eigen_unsymmv_workspace * w =
+ gsl_eigen_unsymmv_alloc (4);
+
+ gsl_eigen_unsymmv (&m.matrix, eval, evec, w);
+
+ gsl_eigen_unsymmv_free (w);
+
+ gsl_eigen_unsymmv_sort (eval, evec,
+ GSL_EIGEN_SORT_ABS_ASC);
+
+ {
+ int i, j;
+
+ for (i = 0; i < 4; i++)
+ {
+ gsl_complex eval_i
+ = gsl_vector_complex_get (eval, i);
+ gsl_vector_complex_view evec_i
+ = gsl_matrix_complex_column (evec, i);
+
+ printf ("eigenvalue = %g + %gi\n",
+ GSL_REAL(eval_i), GSL_IMAG(eval_i));
+ printf ("eigenvector = \n");
+ for (j = 0; j < 4; ++j)
+ {
+ gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
+ printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
+ }
+ }
+ }
+
+ return 0;
+}
Index: eigen/Makefile.am
===================================================================
RCS file: /cvs/gsl/gsl/eigen/Makefile.am,v
retrieving revision 1.19
diff -u -r1.19 Makefile.am
--- eigen/Makefile.am 18 Sep 2006 21:26:59 -0000 1.19
+++ eigen/Makefile.am 21 Sep 2006 18:26:29 -0000
@@ -3,7 +3,7 @@
check_PROGRAMS = test
pkginclude_HEADERS = gsl_eigen.h
-libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c unsymm.c herm.c hermv.c sort.c francis.c schur.c
+libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c unsymm.c unsymmv.c herm.c hermv.c sort.c francis.c schur.c
INCLUDES= -I$(top_builddir)
noinst_HEADERS = qrstep.c
Index: eigen/gsl_eigen.h
===================================================================
RCS file: /cvs/gsl/gsl/eigen/gsl_eigen.h,v
retrieving revision 1.15
diff -u -r1.15 gsl_eigen.h
--- eigen/gsl_eigen.h 18 Sep 2006 21:26:59 -0000 1.15
+++ eigen/gsl_eigen.h 21 Sep 2006 18:26:29 -0000
@@ -22,6 +22,7 @@
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_permutation.h>
#undef __BEGIN_DECLS
#undef __END_DECLS
@@ -103,6 +104,28 @@
gsl_matrix * Z, gsl_eigen_unsymm_workspace * w);
typedef struct {
+ size_t size; /* size of matrices */
+ gsl_vector *work; /* scratch workspace */
+ gsl_vector *work2; /* scratch workspace */
+ gsl_vector *work3; /* scratch workspace */
+
+ gsl_permutation *p2; /* 2 element permutation */
+
+ gsl_matrix *Z; /* pointer to Schur vectors */
+
+ gsl_eigen_unsymm_workspace *unsymm_workspace_p;
+} gsl_eigen_unsymmv_workspace;
+
+gsl_eigen_unsymmv_workspace * gsl_eigen_unsymmv_alloc (const size_t n);
+void gsl_eigen_unsymmv_free (gsl_eigen_unsymmv_workspace * w);
+int gsl_eigen_unsymmv (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_matrix_complex * evec,
+ gsl_eigen_unsymmv_workspace * w);
+int gsl_eigen_unsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_matrix_complex * evec, gsl_matrix * Z,
+ gsl_eigen_unsymmv_workspace * w);
+
+typedef struct {
size_t size;
double * d;
double * sd;
@@ -152,6 +175,10 @@
int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
gsl_eigen_sort_t sort_type);
+int gsl_eigen_unsymmv_sort(gsl_vector_complex * eval,
+ gsl_matrix_complex * evec,
+ gsl_eigen_sort_t sort_type);
+
/* The following functions are obsolete: */
Index: eigen/sort.c
===================================================================
RCS file: /cvs/gsl/gsl/eigen/sort.c,v
retrieving revision 1.4
diff -u -r1.4 sort.c
--- eigen/sort.c 26 Jun 2005 13:27:01 -0000 1.4
+++ eigen/sort.c 21 Sep 2006 18:26:29 -0000
@@ -165,3 +165,69 @@
return GSL_SUCCESS;
}
}
+
+int
+gsl_eigen_unsymmv_sort (gsl_vector_complex * eval,
+ gsl_matrix_complex * evec,
+ gsl_eigen_sort_t sort_type)
+{
+ if (evec->size1 != evec->size2)
+ {
+ GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
+ }
+ else if (eval->size != evec->size1)
+ {
+ GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
+ }
+ else
+ {
+ const size_t N = eval->size;
+ size_t i;
+
+ for (i = 0; i < N - 1; i++)
+ {
+ size_t j;
+ size_t k = i;
+
+ gsl_complex ek = gsl_vector_complex_get (eval, i);
+
+ /* search for something to swap */
+ for (j = i + 1; j < N; j++)
+ {
+ int test;
+ const gsl_complex ej = gsl_vector_complex_get (eval, j);
+
+ switch (sort_type)
+ {
+ case GSL_EIGEN_SORT_ABS_ASC:
+ test = (gsl_complex_abs (ej) < gsl_complex_abs (ek));
+ break;
+ case GSL_EIGEN_SORT_ABS_DESC:
+ test = (gsl_complex_abs (ej) > gsl_complex_abs (ek));
+ break;
+ case GSL_EIGEN_SORT_VAL_ASC:
+ case GSL_EIGEN_SORT_VAL_DESC:
+ default:
+ GSL_ERROR ("invalid sort type", GSL_EINVAL);
+ }
+
+ if (test)
+ {
+ k = j;
+ ek = ej;
+ }
+ }
+
+ if (k != i)
+ {
+ /* swap eigenvalues */
+ gsl_vector_complex_swap_elements (eval, i, k);
+
+ /* swap eigenvectors */
+ gsl_matrix_complex_swap_columns (evec, i, k);
+ }
+ }
+
+ return GSL_SUCCESS;
+ }
+}
Index: eigen/unsymmv.c
===================================================================
RCS file: eigen/unsymmv.c
diff -N eigen/unsymmv.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ eigen/unsymmv.c 21 Sep 2006 18:26:30 -0000
@@ -0,0 +1,1480 @@
+/* eigen/unsymmv.c
+ *
+ * Copyright (C) 2006 Patrick Alken
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <stdlib.h>
+#include <math.h>
+
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_vector_complex.h>
+#include <gsl/gsl_matrix.h>
+
+/*
+ * This module computes the eigenvalues and eigenvectors of a real
+ * unsymmetric matrix.
+ */
+
+#define GSL_UNSYMMV_SMLNUM (2.0 * GSL_DBL_MIN)
+#define GSL_UNSYMMV_BIGNUM ((1.0 - GSL_DBL_EPSILON) / GSL_UNSYMMV_SMLNUM)
+
+static void unsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
+ gsl_vector_complex *eval,
+ gsl_matrix_complex *evec,
+ gsl_eigen_unsymmv_workspace *w);
+static inline void unsymmv_solve_equation(gsl_matrix *A, double z,
+ gsl_vector *b, gsl_vector *x,
+ double *s, double *xnorm,
+ double smin,
+ gsl_eigen_unsymmv_workspace *w);
+static inline void unsymmv_solve_equation_z(gsl_matrix *A, gsl_complex *z,
+ gsl_vector_complex *b,
+ gsl_vector_complex *x,
+ double *s, double *xnorm,
+ double smin,
+ gsl_eigen_unsymmv_workspace *w);
+static void unsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
+ gsl_matrix_complex *evec);
+
+/*
+gsl_eigen_unsymmv_alloc()
+
+Allocate a workspace for solving the unsymmetric eigenvalue problem.
+The size of this workspace is O(5n).
+
+Inputs: n - size of matrices
+
+Return: pointer to workspace
+*/
+
+gsl_eigen_unsymmv_workspace *
+gsl_eigen_unsymmv_alloc(const size_t n)
+{
+ gsl_eigen_unsymmv_workspace *w;
+
+ if (n == 0)
+ {
+ GSL_ERROR_NULL ("matrix dimension must be positive integer",
+ GSL_EINVAL);
+ }
+
+ w = (gsl_eigen_unsymmv_workspace *)
+ malloc (sizeof (gsl_eigen_unsymmv_workspace));
+
+ if (w == 0)
+ {
+ GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
+ }
+
+ w->size = n;
+ w->Z = NULL;
+ w->unsymm_workspace_p = gsl_eigen_unsymm_alloc(n);
+
+ if (w->unsymm_workspace_p == 0)
+ {
+ GSL_ERROR_NULL ("failed to allocate space for unsymm workspace", GSL_ENOMEM);
+ }
+
+ w->p2 = gsl_permutation_alloc(2);
+ if (w->p2 == 0)
+ {
+ GSL_ERROR_NULL ("failed to allocate space for unsymmv permutation", GSL_ENOMEM);
+ }
+
+ /*
+ * set parameters to compute the full Schur form T and balance
+ * the matrices
+ */
+ gsl_eigen_unsymm_params(1, 1, w->unsymm_workspace_p);
+
+ w->work = gsl_vector_alloc(n);
+ w->work2 = gsl_vector_alloc(n);
+ w->work3 = gsl_vector_alloc(n);
+ if (w->work == 0 || w->work2 == 0 || w->work3 == 0)
+ {
+ GSL_ERROR_NULL ("failed to allocate space for unsymmv additional workspace", GSL_ENOMEM);
+ }
+
+ return (w);
+} /* gsl_eigen_unsymmv_alloc() */
+
+/*
+gsl_eigen_unsymmv_free()
+ Free workspace w
+*/
+
+void
+gsl_eigen_unsymmv_free (gsl_eigen_unsymmv_workspace * w)
+{
+ gsl_eigen_unsymm_free(w->unsymm_workspace_p);
+ gsl_vector_free(w->work);
+ gsl_vector_free(w->work2);
+ gsl_vector_free(w->work3);
+ gsl_permutation_free(w->p2);
+
+ free(w);
+} /* gsl_eigen_unsymmv_free() */
+
+/*
+gsl_eigen_unsymmv()
+
+Solve the unsymmetric eigensystem problem
+
+A x = \lambda x
+
+for the eigenvalues \lambda and right eigenvectors x
+
+Inputs: A - general real matrix
+ eval - where to store eigenvalues
+ evec - where to store eigenvectors
+ w - workspace
+
+Return: number of eigenvalues found; if this number is less than
+ N, the QR algorithm failed and no eigenvectors are computed
+*/
+
+int
+gsl_eigen_unsymmv (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_matrix_complex * evec,
+ gsl_eigen_unsymmv_workspace * w)
+{
+ const size_t N = A->size1;
+
+ /* check matrix and vector sizes */
+
+ if (N != A->size2)
+ {
+ GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
+ }
+ else if (eval->size != N)
+ {
+ GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
+ }
+ else if (evec->size1 != evec->size2)
+ {
+ GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
+ }
+ else if (evec->size1 != N)
+ {
+ GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
+ }
+ else
+ {
+ int s;
+ gsl_matrix Z;
+
+ /*
+ * We need a place to store the Schur vectors, so we will
+ * treat evec as a real matrix and store them in the left
+ * half - the factor of 2 in the tda corresponds to the
+ * complex multiplicity
+ */
+ Z.size1 = N;
+ Z.size2 = N;
+ Z.tda = 2 * N;
+ Z.data = evec->data;
+ Z.block = 0;
+ Z.owner = 0;
+
+ /* compute eigenvalues, Schur form, and Schur vectors */
+ s = gsl_eigen_unsymm_Z(A, eval, &Z, w->unsymm_workspace_p);
+
+ if (w->Z)
+ {
+ /*
+ * save the Schur vectors in user supplied matrix, since
+ * they will be destroyed when computing eigenvectors
+ */
+ gsl_matrix_memcpy(w->Z, &Z);
+ }
+
+ /* only compute eigenvectors if we found all eigenvalues */
+ if ((size_t) s == w->size)
+ {
+ /* compute eigenvectors */
+ unsymmv_get_right_eigenvectors(A, &Z, eval, evec, w);
+
+ /* normalize so that Euclidean norm is 1 */
+ unsymmv_normalize_eigenvectors(eval, evec);
+ }
+
+ return s;
+ }
+} /* gsl_eigen_unsymmv() */
+
+/*
+gsl_eigen_unsymmv_Z()
+ Compute eigenvalues and eigenvectors of a real unsymmetric matrix
+and also save the Schur vectors. See comments in gsl_eigen_unsymm_Z
+for more information.
+
+Inputs: A - real unsymmetric matrix
+ eval - where to store eigenvalues
+ evec - where to store eigenvectors
+ Z - where to store Schur vectors
+ w - unsymmv workspace
+*/
+
+int
+gsl_eigen_unsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_matrix_complex * evec, gsl_matrix * Z,
+ gsl_eigen_unsymmv_workspace * w)
+{
+ /* check matrix and vector sizes */
+
+ if (A->size1 != A->size2)
+ {
+ GSL_ERROR ("matrix must be square to compute eigenvalues/eigenvectors", GSL_ENOTSQR);
+ }
+ else if (eval->size != A->size1)
+ {
+ GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
+ }
+ else if (evec->size1 != evec->size2)
+ {
+ GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
+ }
+ else if (evec->size1 != A->size1)
+ {
+ GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
+ }
+ else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
+ {
+ GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
+ }
+ else
+ {
+ int s;
+
+ w->Z = Z;
+
+ s = gsl_eigen_unsymmv(A, eval, evec, w);
+
+ w->Z = NULL;
+
+ return s;
+ }
+} /* gsl_eigen_unsymmv_Z() */
+
+/********************************************
+ * INTERNAL ROUTINES *
+ ********************************************/
+
+/*
+unsymmv_get_right_eigenvectors()
+ Compute the right eigenvectors of the Schur form T and then
+backtransform them using the Schur vectors to get right eigenvectors of
+the original matrix.
+
+Inputs: T - Schur form
+ Z - Schur vectors
+ eval - where to store eigenvalues (to ensure that the
+ correct eigenvalue is stored in the same position
+ as the eigenvectors)
+ evec - where to store eigenvectors
+ w - unsymmv workspace
+
+Return: none
+
+Notes: 1) based on LAPACK routine DTREVC - the algorithm used is
+ backsubstitution on the upper quasi triangular system T
+ followed by backtransformation by Z to get vectors of the
+ original matrix.
+
+ 2) The Schur vectors in Z are destroyed and replaced with
+ eigenvectors stored with the same storage scheme as DTREVC.
+ The eigenvectors are also stored in 'evec'
+
+ 3) The matrix T is unchanged on output
+
+ 4) Each eigenvector is normalized so that the element of
+ largest magnitude has magnitude 1; here the magnitude of
+ a complex number (x,y) is taken to be |x| + |y|
+*/
+
+static void
+unsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
+ gsl_vector_complex *eval,
+ gsl_matrix_complex *evec,
+ gsl_eigen_unsymmv_workspace *w)
+{
+ const size_t N = T->size1;
+ const double smlnum = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
+ const double bignum = (1.0 - GSL_DBL_EPSILON) / smlnum;
+ int i; /* looping */
+ size_t iu, /* looping */
+ ju,
+ ii;
+ gsl_complex lambda; /* current eigenvalue */
+ double lambda_re, /* Re(lambda) */
+ lambda_im; /* Im(lambda) */
+ gsl_matrix_view Tv, /* temporary views */
+ Zv;
+ gsl_vector_view y, /* temporary views */
+ y2,
+ ev,
+ ev2;
+ double dat[4], /* scratch arrays */
+ dat_X[4];
+ double scale; /* scale factor */
+ double xnorm; /* |X| */
+ gsl_vector_complex_view ecol, /* column of evec */
+ ecol2;
+ int complex_pair; /* complex eigenvalue pair? */
+ double smin;
+
+ /*
+ * Compute 1-norm of each column of upper triangular part of T
+ * to control overflow in triangular solver
+ */
+
+ gsl_vector_set(w->work3, 0, 0.0);
+ for (ju = 1; ju < N; ++ju)
+ {
+ gsl_vector_set(w->work3, ju, 0.0);
+ for (iu = 0; iu < ju; ++iu)
+ {
+ gsl_vector_set(w->work3, ju,
+ gsl_vector_get(w->work3, ju) +
+ fabs(gsl_matrix_get(T, iu, ju)));
+ }
+ }
+
+ for (i = (int) N - 1; i >= 0; --i)
+ {
+ iu = (size_t) i;
+
+ /* get current eigenvalue and store it in lambda */
+ lambda_re = gsl_matrix_get(T, iu, iu);
+
+ if (iu != 0 && gsl_matrix_get(T, iu, iu - 1) != 0.0)
+ {
+ lambda_im = sqrt(fabs(gsl_matrix_get(T, iu, iu - 1))) *
+ sqrt(fabs(gsl_matrix_get(T, iu - 1, iu)));
+ }
+ else
+ {
+ lambda_im = 0.0;
+ }
+
+ GSL_SET_COMPLEX(&lambda, lambda_re, lambda_im);
+
+ smin = GSL_MAX(GSL_DBL_EPSILON * (fabs(lambda_re) + fabs(lambda_im)),
+ smlnum);
+ smin = GSL_MAX(smin, GSL_UNSYMMV_SMLNUM);
+
+ if (lambda_im == 0.0)
+ {
+ int k, l;
+ gsl_vector_view bv, xv;
+
+ /* real eigenvector */
+
+ /*
+ * The ordering of eigenvalues in 'eval' is arbitrary and
+ * does not necessarily follow the Schur form T, so store
+ * lambda in the right slot in eval to ensure it corresponds
+ * to the eigenvector we are about to compute
+ */
+ gsl_vector_complex_set(eval, iu, lambda);
+
+ /*
+ * We need to solve the system:
+ *
+ * (T(1:iu-1, 1:iu-1) - lambda*I)*X = -T(1:iu-1,iu)
+ */
+
+ /* construct right hand side */
+ for (k = 0; k < i; ++k)
+ {
+ gsl_vector_set(w->work,
+ (size_t) k,
+ -gsl_matrix_get(T, (size_t) k, iu));
+ }
+
+ gsl_vector_set(w->work, iu, 1.0);
+
+ for (l = i - 1; l >= 0; --l)
+ {
+ size_t lu = (size_t) l;
+
+ if (lu == 0)
+ complex_pair = 0;
+ else
+ complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;
+
+ if (!complex_pair)
+ {
+ double x;
+
+ /*
+ * 1-by-1 diagonal block - solve the system:
+ *
+ * (T_{ll} - lambda)*x = -T_{l(iu)}
+ */
+
+ Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
+ bv = gsl_vector_view_array(dat, 1);
+ gsl_vector_set(&bv.vector, 0,
+ gsl_vector_get(w->work, lu));
+ xv = gsl_vector_view_array(dat_X, 1);
+
+ unsymmv_solve_equation(&Tv.matrix,
+ lambda_re,
+ &bv.vector,
+ &xv.vector,
+ &scale,
+ &xnorm,
+ smin,
+ w);
+
+ /* scale x to avoid overflow */
+ x = gsl_vector_get(&xv.vector, 0);
+ if (xnorm > 1.0)
+ {
+ if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
+ {
+ x /= xnorm;
+ scale /= xnorm;
+ }
+ }
+
+ if (scale != 1.0)
+ {
+ gsl_vector_view wv;
+
+ wv = gsl_vector_subvector(w->work, 0, iu + 1);
+ gsl_blas_dscal(scale, &wv.vector);
+ }
+
+ gsl_vector_set(w->work, lu, x);
+
+ if (lu > 0)
+ {
+ gsl_vector_view v1, v2;
+
+ /* update right hand side */
+
+ v1 = gsl_matrix_column(T, lu);
+ v1 = gsl_vector_subvector(&v1.vector, 0, lu);
+
+ v2 = gsl_vector_subvector(w->work, 0, lu);
+
+ gsl_blas_daxpy(-x, &v1.vector, &v2.vector);
+ } /* if (l > 0) */
+ } /* if (!complex_pair) */
+ else
+ {
+ double x11, x21;
+
+ /*
+ * 2-by-2 diagonal block
+ */
+
+ Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
+ bv = gsl_vector_view_array(dat, 2);
+ gsl_vector_set(&bv.vector, 0,
+ gsl_vector_get(w->work, lu - 1));
+ gsl_vector_set(&bv.vector, 1,
+ gsl_vector_get(w->work, lu));
+ xv = gsl_vector_view_array(dat_X, 2);
+
+ unsymmv_solve_equation(&Tv.matrix,
+ lambda_re,
+ &bv.vector,
+ &xv.vector,
+ &scale,
+ &xnorm,
+ smin,
+ w);
+
+ /* scale X(1,1) and X(2,1) to avoid overflow */
+ x11 = gsl_vector_get(&xv.vector, 0);
+ x21 = gsl_vector_get(&xv.vector, 1);
+
+ if (xnorm > 1.0)
+ {
+ double beta;
+
+ beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
+ gsl_vector_get(w->work3, lu));
+ if (beta > bignum / xnorm)
+ {
+ x11 /= xnorm;
+ x21 /= xnorm;
+ scale /= xnorm;
+ }
+ }
+
+ /* scale if necessary */
+ if (scale != 1.0)
+ {
+ gsl_vector_view wv;
+
+ wv = gsl_vector_subvector(w->work, 0, iu + 1);
+ gsl_blas_dscal(scale, &wv.vector);
+ }
+
+ gsl_vector_set(w->work, lu - 1, x11);
+ gsl_vector_set(w->work, lu, x21);
+
+ /* update right hand side */
+ if (lu > 1)
+ {
+ gsl_vector_view v1, v2;
+
+ v1 = gsl_matrix_column(T, lu - 1);
+ v1 = gsl_vector_subvector(&v1.vector, 0, lu - 1);
+ v2 = gsl_vector_subvector(w->work, 0, lu - 1);
+ gsl_blas_daxpy(-x11, &v1.vector, &v2.vector);
+
+ v1 = gsl_matrix_column(T, lu);
+ v1 = gsl_vector_subvector(&v1.vector, 0, lu - 1);
+ gsl_blas_daxpy(-x21, &v1.vector, &v2.vector);
+ }
+
+ --l;
+ } /* if (complex_pair) */
+ } /* for (l = i - 1; l >= 0; --l) */
+
+ /*
+ * At this point, w->work is an eigenvector of the
+ * Schur form T. To get an eigenvector of the original
+ * matrix, we multiply on the left by Z, the matrix of
+ * Schur vectors
+ */
+
+ ecol = gsl_matrix_complex_column(evec, iu);
+ y = gsl_matrix_column(Z, iu);
+
+ if (iu > 0)
+ {
+ gsl_vector_view x;
+
+ Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu);
+
+ x = gsl_vector_subvector(w->work, 0, iu);
+
+ /* compute Z * w->work and store it in Z(:,iu) */
+ gsl_blas_dgemv(CblasNoTrans,
+ 1.0,
+ &Zv.matrix,
+ &x.vector,
+ gsl_vector_get(w->work, iu),
+ &y.vector);
+ } /* if (iu > 0) */
+
+ /* store eigenvector into evec */
+
+ ev = gsl_vector_complex_real(&ecol.vector);
+ ev2 = gsl_vector_complex_imag(&ecol.vector);
+
+ scale = 0.0;
+ for (ii = 0; ii < N; ++ii)
+ {
+ double a = gsl_vector_get(&y.vector, ii);
+
+ /* store real part of eigenvector */
+ gsl_vector_set(&ev.vector, ii, a);
+
+ /* set imaginary part to 0 */
+ gsl_vector_set(&ev2.vector, ii, 0.0);
+
+ if (fabs(a) > scale)
+ scale = fabs(a);
+ }
+
+ if (scale != 0.0)
+ scale = 1.0 / scale;
+
+ /* scale by magnitude of largest element */
+ gsl_blas_dscal(scale, &ev.vector);
+ } /* if (GSL_IMAG(lambda) == 0.0) */
+ else
+ {
+ gsl_vector_complex_view bv, xv;
+ size_t k;
+ int l;
+ gsl_complex lambda2;
+
+ /* complex eigenvector */
+
+ /*
+ * Store the complex conjugate eigenvalues in the right
+ * slots in eval
+ */
+ GSL_SET_REAL(&lambda2, GSL_REAL(lambda));
+ GSL_SET_IMAG(&lambda2, -GSL_IMAG(lambda));
+ gsl_vector_complex_set(eval, iu - 1, lambda);
+ gsl_vector_complex_set(eval, iu, lambda2);
+
+ /*
+ * First solve:
+ *
+ * [ T(i:i+1,i:i+1) - lambda*I ] * X = 0
+ */
+
+ if (fabs(gsl_matrix_get(T, iu - 1, iu)) >=
+ fabs(gsl_matrix_get(T, iu, iu - 1)))
+ {
+ gsl_vector_set(w->work, iu - 1, 1.0);
+ gsl_vector_set(w->work2, iu,
+ lambda_im / gsl_matrix_get(T, iu - 1, iu));
+ }
+ else
+ {
+ gsl_vector_set(w->work, iu - 1,
+ -lambda_im / gsl_matrix_get(T, iu, iu - 1));
+ gsl_vector_set(w->work2, iu, 1.0);
+ }
+ gsl_vector_set(w->work, iu, 0.0);
+ gsl_vector_set(w->work2, iu - 1, 0.0);
+
+ /* construct right hand side */
+ for (k = 0; k < iu - 1; ++k)
+ {
+ gsl_vector_set(w->work, k,
+ -gsl_vector_get(w->work, iu - 1) *
+ gsl_matrix_get(T, k, iu - 1));
+ gsl_vector_set(w->work2, k,
+ -gsl_vector_get(w->work2, iu) *
+ gsl_matrix_get(T, k, iu));
+ }
+
+ /*
+ * We must solve the upper quasi-triangular system:
+ *
+ * [ T(1:i-2,1:i-2) - lambda*I ] * X = s*(work + i*work2)
+ */
+
+ for (l = i - 2; l >= 0; --l)
+ {
+ size_t lu = (size_t) l;
+
+ if (lu == 0)
+ complex_pair = 0;
+ else
+ complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;
+
+ if (!complex_pair)
+ {
+ gsl_complex bval;
+ gsl_complex x;
+
+ /*
+ * 1-by-1 diagonal block - solve the system:
+ *
+ * (T_{ll} - lambda)*x = work + i*work2
+ */
+
+ Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
+ bv = gsl_vector_complex_view_array(dat, 1);
+ xv = gsl_vector_complex_view_array(dat_X, 1);
+
+ GSL_SET_COMPLEX(&bval,
+ gsl_vector_get(w->work, lu),
+ gsl_vector_get(w->work2, lu));
+ gsl_vector_complex_set(&bv.vector, 0, bval);
+
+ unsymmv_solve_equation_z(&Tv.matrix,
+ &lambda,
+ &bv.vector,
+ &xv.vector,
+ &scale,
+ &xnorm,
+ smin,
+ w);
+
+ if (xnorm > 1.0)
+ {
+ if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
+ {
+ gsl_blas_zdscal(1.0/xnorm, &xv.vector);
+ scale /= xnorm;
+ }
+ }
+
+ /* scale if necessary */
+ if (scale != 1.0)
+ {
+ gsl_vector_view wv;
+
+ wv = gsl_vector_subvector(w->work, 0, iu + 1);
+ gsl_blas_dscal(scale, &wv.vector);
+ wv = gsl_vector_subvector(w->work2, 0, iu + 1);
+ gsl_blas_dscal(scale, &wv.vector);
+ }
+
+ x = gsl_vector_complex_get(&xv.vector, 0);
+ gsl_vector_set(w->work, lu, GSL_REAL(x));
+ gsl_vector_set(w->work2, lu, GSL_IMAG(x));
+
+ /* update the right hand side */
+ if (lu > 0)
+ {
+ gsl_vector_view v1, v2;
+
+ v1 = gsl_matrix_column(T, lu);
+ v1 = gsl_vector_subvector(&v1.vector, 0, lu);
+ v2 = gsl_vector_subvector(w->work, 0, lu);
+ gsl_blas_daxpy(-GSL_REAL(x), &v1.vector, &v2.vector);
+
+ v2 = gsl_vector_subvector(w->work2, 0, lu);
+ gsl_blas_daxpy(-GSL_IMAG(x), &v1.vector, &v2.vector);
+ } /* if (lu > 0) */
+ } /* if (!complex_pair) */
+ else
+ {
+ gsl_complex b1, b2, x1, x2;
+
+ /*
+ * 2-by-2 diagonal block - solve the system
+ */
+
+ Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
+ bv = gsl_vector_complex_view_array(dat, 2);
+ xv = gsl_vector_complex_view_array(dat_X, 2);
+
+ GSL_SET_COMPLEX(&b1,
+ gsl_vector_get(w->work, lu - 1),
+ gsl_vector_get(w->work2, lu - 1));
+ GSL_SET_COMPLEX(&b2,
+ gsl_vector_get(w->work, lu),
+ gsl_vector_get(w->work2, lu));
+ gsl_vector_complex_set(&bv.vector, 0, b1);
+ gsl_vector_complex_set(&bv.vector, 1, b2);
+
+ unsymmv_solve_equation_z(&Tv.matrix,
+ &lambda,
+ &bv.vector,
+ &xv.vector,
+ &scale,
+ &xnorm,
+ smin,
+ w);
+
+ x1 = gsl_vector_complex_get(&xv.vector, 0);
+ x2 = gsl_vector_complex_get(&xv.vector, 1);
+
+ if (xnorm > 1.0)
+ {
+ double beta;
+
+ beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
+ gsl_vector_get(w->work3, lu));
+ if (beta > bignum / xnorm)
+ {
+ gsl_blas_zdscal(1.0/xnorm, &xv.vector);
+ scale /= xnorm;
+ }
+ }
+
+ /* scale if necessary */
+ if (scale != 1.0)
+ {
+ gsl_vector_view wv;
+
+ wv = gsl_vector_subvector(w->work, 0, iu + 1);
+ gsl_blas_dscal(scale, &wv.vector);
+ wv = gsl_vector_subvector(w->work2, 0, iu + 1);
+ gsl_blas_dscal(scale, &wv.vector);
+ }
+ gsl_vector_set(w->work, lu - 1, GSL_REAL(x1));
+ gsl_vector_set(w->work, lu, GSL_REAL(x2));
+ gsl_vector_set(w->work2, lu - 1, GSL_IMAG(x1));
+ gsl_vector_set(w->work2, lu, GSL_IMAG(x2));
+
+ /* update right hand side */
+ if (lu > 1)
+ {
+ gsl_vector_view v1, v2, v3, v4;
+
+ v1 = gsl_matrix_column(T, lu - 1);
+ v1 = gsl_vector_subvector(&v1.vector, 0, lu - 1);
+ v4 = gsl_matrix_column(T, lu);
+ v4 = gsl_vector_subvector(&v4.vector, 0, lu - 1);
+ v2 = gsl_vector_subvector(w->work, 0, lu - 1);
+ v3 = gsl_vector_subvector(w->work2, 0, lu - 1);
+
+ gsl_blas_daxpy(-GSL_REAL(x1), &v1.vector, &v2.vector);
+ gsl_blas_daxpy(-GSL_REAL(x2), &v4.vector, &v2.vector);
+ gsl_blas_daxpy(-GSL_IMAG(x1), &v1.vector, &v3.vector);
+ gsl_blas_daxpy(-GSL_IMAG(x2), &v4.vector, &v3.vector);
+ } /* if (lu > 1) */
+
+ --l;
+ } /* if (complex_pair) */
+ } /* for (l = i - 2; l >= 0; --l) */
+
+ /*
+ * At this point, work + i*work2 is an eigenvector
+ * of T - backtransform to get an eigenvector of the
+ * original matrix
+ */
+
+ y = gsl_matrix_column(Z, iu - 1);
+ y2 = gsl_matrix_column(Z, iu);
+
+ if (iu > 1)
+ {
+ gsl_vector_view x;
+
+ /* compute real part of eigenvectors */
+
+ Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu - 1);
+ x = gsl_vector_subvector(w->work, 0, iu - 1);
+
+ gsl_blas_dgemv(CblasNoTrans,
+ 1.0,
+ &Zv.matrix,
+ &x.vector,
+ gsl_vector_get(w->work, iu - 1),
+ &y.vector);
+
+
+ /* now compute the imaginary part */
+ x = gsl_vector_subvector(w->work2, 0, iu - 1);
+
+ gsl_blas_dgemv(CblasNoTrans,
+ 1.0,
+ &Zv.matrix,
+ &x.vector,
+ gsl_vector_get(w->work2, iu),
+ &y2.vector);
+ }
+ else
+ {
+ gsl_blas_dscal(gsl_vector_get(w->work, iu - 1), &y.vector);
+ gsl_blas_dscal(gsl_vector_get(w->work2, iu), &y2.vector);
+ }
+
+ /*
+ * Now store the eigenvectors into evec - the real parts
+ * are Z(:,iu - 1) and the imaginary parts are
+ * +/- Z(:,iu)
+ */
+
+ /* get views of the two eigenvector slots */
+ ecol = gsl_matrix_complex_column(evec, iu - 1);
+ ecol2 = gsl_matrix_complex_column(evec, iu);
+
+ /*
+ * save imaginary part first as it may get overwritten
+ * when copying the real part due to our storage scheme
+ * in Z/evec
+ */
+ ev = gsl_vector_complex_imag(&ecol.vector);
+ ev2 = gsl_vector_complex_imag(&ecol2.vector);
+ scale = 0.0;
+ for (ii = 0; ii < N; ++ii)
+ {
+ double a = gsl_vector_get(&y2.vector, ii);
+
+ scale = GSL_MAX(scale,
+ fabs(a) + fabs(gsl_vector_get(&y.vector, ii)));
+
+ gsl_vector_set(&ev.vector, ii, a);
+ gsl_vector_set(&ev2.vector, ii, -a);
+ }
+
+ /* now save the real part */
+ ev = gsl_vector_complex_real(&ecol.vector);
+ ev2 = gsl_vector_complex_real(&ecol2.vector);
+ for (ii = 0; ii < N; ++ii)
+ {
+ double a = gsl_vector_get(&y.vector, ii);
+
+ gsl_vector_set(&ev.vector, ii, a);
+ gsl_vector_set(&ev2.vector, ii, a);
+ }
+
+ if (scale != 0.0)
+ scale = 1.0 / scale;
+
+ /* scale by largest element magnitude */
+
+ gsl_blas_zdscal(scale, &ecol.vector);
+ gsl_blas_zdscal(scale, &ecol2.vector);
+
+ /*
+ * decrement i since we took care of two eigenvalues at
+ * the same time
+ */
+ --i;
+ } /* if (GSL_IMAG(lambda) != 0.0) */
+ } /* for (i = (int) N - 1; i >= 0; --i) */
+} /* unsymmv_get_right_eigenvectors() */
+
+/*
+unsymmv_solve_equation()
+
+ Solve the equation which comes up in the back substitution
+when computing eigenvectors corresponding to real eigenvalues.
+The equation that is solved is:
+
+(A - z*I)*x = s*b
+
+where
+
+A is n-by-n with n = 1 or 2
+b and x are n-by-1 real vectors
+s is a scaling factor set by this function to prevent overflow in x
+
+Inputs: A - square matrix (n-by-n)
+ z - real scalar (eigenvalue)
+ b - right hand side vector
+ x - (output) where to store solution
+ s - (output) scale factor
+ xnorm - (output) infinity norm of X
+ smin - lower bound on singular values of A - if A - z*I
+ is less than this value, we'll use smin*I instead.
+ This value should be a safe distance above underflow.
+ w - unsymmv workspace
+
+Notes: 1) A and b are not changed on output
+ 2) Based on lapack routine DLALN2
+*/
+
+static inline void
+unsymmv_solve_equation(gsl_matrix *A, double z, gsl_vector *b,
+ gsl_vector *x, double *s, double *xnorm,
+ double smin, gsl_eigen_unsymmv_workspace *w)
+{
+ size_t N = A->size1;
+ double bnorm;
+ double scale = 1.0;
+
+ if (N == 1)
+ {
+ double c, /* denominator */
+ cnorm; /* |c| */
+
+ /*
+ * we have a 1-by-1 (real) scalar system to solve:
+ *
+ * (a - z)*x = b
+ * with z real
+ */
+
+ /* c = a - z */
+ c = gsl_matrix_get(A, 0, 0) - z;
+ cnorm = fabs(c);
+
+ if (cnorm < smin)
+ {
+ /* set c = smin*I */
+ c = smin;
+ cnorm = smin;
+ }
+
+ /* check scaling for x = b / c */
+ bnorm = fabs(gsl_vector_get(b, 0));
+ if (cnorm < 1.0 && bnorm > 1.0)
+ {
+ if (bnorm > GSL_UNSYMMV_BIGNUM*cnorm)
+ scale = 1.0 / bnorm;
+ }
+
+ /* compute x */
+ gsl_vector_set(x, 0, gsl_vector_get(b, 0) * scale / c);
+ *xnorm = fabs(gsl_vector_get(x, 0));
+ } /* if (N == 1) */
+ else
+ {
+ double cr[2][2];
+ double *crv;
+ double cmax;
+ size_t icmax, j;
+ double bval1, bval2;
+ double ur11, ur12, ur22, ur11r;
+ double cr21, cr22;
+ double lr21;
+ double b1, b2, bbnd;
+ double x1, x2;
+ double temp;
+ size_t ipivot[4][4] = { { 0, 1, 2, 3 },
+ { 1, 0, 3, 2 },
+ { 2, 3, 0, 1 },
+ { 3, 2, 1, 0 } };
+ int rswap[4] = { 0, 1, 0, 1 };
+ int zswap[4] = { 0, 0, 1, 1 };
+
+ /*
+ * we have a 2-by-2 real system to solve:
+ *
+ * [ A11 - z A12 ] [ x1 ] = [ b1 ]
+ * [ A21 A22 - z ] [ x2 ] [ b2 ]
+ *
+ * (z real)
+ */
+
+ crv = (double *) cr;
+
+ /*
+ * compute the real part of C = A - z*I - use column ordering
+ * here since porting from lapack
+ */
+ cr[0][0] = gsl_matrix_get(A, 0, 0) - z;
+ cr[1][1] = gsl_matrix_get(A, 1, 1) - z;
+ cr[0][1] = gsl_matrix_get(A, 1, 0);
+ cr[1][0] = gsl_matrix_get(A, 0, 1);
+
+ /* find the largest element in C */
+ cmax = 0.0;
+ icmax = 0;
+ for (j = 0; j < 4; ++j)
+ {
+ if (fabs(crv[j]) > cmax)
+ {
+ cmax = fabs(crv[j]);
+ icmax = j;
+ }
+ }
+
+ bval1 = gsl_vector_get(b, 0);
+ bval2 = gsl_vector_get(b, 1);
+
+ /* if norm(C) < smin, use smin*I */
+
+ if (cmax < smin)
+ {
+ bnorm = GSL_MAX(fabs(bval1), fabs(bval2));
+ if (smin < 1.0 && bnorm > 1.0)
+ {
+ if (bnorm > GSL_UNSYMMV_BIGNUM*smin)
+ scale = 1.0 / bnorm;
+ }
+ temp = scale / smin;
+ gsl_vector_set(x, 0, temp * bval1);
+ gsl_vector_set(x, 1, temp * bval2);
+ *xnorm = temp * bnorm;
+ *s = scale;
+ return;
+ }
+
+ /* gaussian elimination with complete pivoting */
+ ur11 = crv[icmax];
+ cr21 = crv[ipivot[1][icmax]];
+ ur12 = crv[ipivot[2][icmax]];
+ cr22 = crv[ipivot[3][icmax]];
+ ur11r = 1.0 / ur11;
+ lr21 = ur11r * cr21;
+ ur22 = cr22 - ur12 * lr21;
+
+ /* if smaller pivot < smin, use smin */
+ if (fabs(ur22) < smin)
+ ur22 = smin;
+
+ if (rswap[icmax])
+ {
+ b1 = bval2;
+ b2 = bval1;
+ }
+ else
+ {
+ b1 = bval1;
+ b2 = bval2;
+ }
+
+ b2 -= lr21 * b1;
+ bbnd = GSL_MAX(fabs(b1 * (ur22 * ur11r)), fabs(b2));
+ if (bbnd > 1.0 && fabs(ur22) < 1.0)
+ {
+ if (bbnd >= GSL_UNSYMMV_BIGNUM * fabs(ur22))
+ scale = 1.0 / bbnd;
+ }
+
+ x2 = (b2 * scale) / ur22;
+ x1 = (scale * b1) * ur11r - x2 * (ur11r * ur12);
+ if (zswap[icmax])
+ {
+ gsl_vector_set(x, 0, x2);
+ gsl_vector_set(x, 1, x1);
+ }
+ else
+ {
+ gsl_vector_set(x, 0, x1);
+ gsl_vector_set(x, 1, x2);
+ }
+
+ *xnorm = GSL_MAX(fabs(x1), fabs(x2));
+
+ /* further scaling if norm(A) norm(X) > overflow */
+ if (*xnorm > 1.0 && cmax > 1.0)
+ {
+ if (*xnorm > GSL_UNSYMMV_BIGNUM / cmax)
+ {
+ temp = cmax / GSL_UNSYMMV_BIGNUM;
+ gsl_blas_dscal(temp, x);
+ *xnorm *= temp;
+ scale *= temp;
+ }
+ }
+ } /* if (N == 2) */
+
+ *s = scale;
+} /* unsymmv_solve_equation() */
+
+/*
+unsymmv_solve_equation_z()
+
+ Solve the equation which comes up in the back substitution
+when computing eigenvectors corresponding to complex eigenvalues.
+The equation that is solved is:
+
+(A - z*I)*x = s*b
+
+where
+
+A is n-by-n with n = 1 or 2
+b and x are n-by-1 complex vectors
+s is a scaling factor set by this function to prevent overflow in x
+
+Inputs: A - square matrix (n-by-n)
+ z - complex scalar (eigenvalue)
+ b - right hand side vector
+ x - (output) where to store solution
+ s - (output) scale factor
+ xnorm - (output) infinity norm of X
+ smin - lower bound on singular values of A - if A - z*I
+ is less than this value, we'll use smin*I instead.
+ This value should be a safe distance above underflow.
+ w - unsymmv workspace
+
+Notes: 1) A and b are not changed on output
+ 2) Based on lapack routine DLALN2
+*/
+
+static inline void
+unsymmv_solve_equation_z(gsl_matrix *A, gsl_complex *z,
+ gsl_vector_complex *b, gsl_vector_complex *x,
+ double *s, double *xnorm, double smin,
+ gsl_eigen_unsymmv_workspace *w)
+{
+ size_t N = A->size1;
+ double scale = 1.0;
+ double bnorm;
+
+ if (N == 1)
+ {
+ double cr, /* denominator */
+ ci,
+ cnorm; /* |c| */
+ gsl_complex bval, c, xval, tmp;
+
+ /*
+ * we have a 1-by-1 (complex) scalar system to solve:
+ *
+ * (a - z)*x = b
+ * (z is complex, a is real)
+ */
+
+ /* c = a - z */
+ cr = gsl_matrix_get(A, 0, 0) - GSL_REAL(*z);
+ ci = -GSL_IMAG(*z);
+ cnorm = fabs(cr) + fabs(ci);
+
+ if (cnorm < smin)
+ {
+ /* set c = smin*I */
+ cr = smin;
+ ci = 0.0;
+ cnorm = smin;
+ }
+
+ /* check scaling for x = b / c */
+ bval = gsl_vector_complex_get(b, 0);
+ bnorm = fabs(GSL_REAL(bval)) + fabs(GSL_IMAG(bval));
+ if (cnorm < 1.0 && bnorm > 1.0)
+ {
+ if (bnorm > GSL_UNSYMMV_BIGNUM*cnorm)
+ scale = 1.0 / bnorm;
+ }
+
+ /* compute x */
+ GSL_SET_COMPLEX(&tmp, scale*GSL_REAL(bval), scale*GSL_IMAG(bval));
+ GSL_SET_COMPLEX(&c, cr, ci);
+ xval = gsl_complex_div(tmp, c);
+
+ gsl_vector_complex_set(x, 0, xval);
+
+ *xnorm = fabs(GSL_REAL(xval)) + fabs(GSL_IMAG(xval));
+ } /* if (N == 1) */
+ else
+ {
+ double cr[2][2], ci[2][2];
+ double *civ, *crv;
+ double cmax;
+ gsl_complex bval1, bval2;
+ gsl_complex xval1, xval2;
+ double xr1, xi1;
+ size_t icmax;
+ size_t j;
+ double temp;
+ double ur11, ur12, ur22, ui11, ui12, ui22, ur11r, ui11r;
+ double ur12s, ui12s;
+ double u22abs;
+ double lr21, li21;
+ double cr21, cr22, ci21, ci22;
+ double br1, bi1, br2, bi2, bbnd;
+ gsl_complex b1, b2;
+ size_t ipivot[4][4] = { { 0, 1, 2, 3 },
+ { 1, 0, 3, 2 },
+ { 2, 3, 0, 1 },
+ { 3, 2, 1, 0 } };
+ int rswap[4] = { 0, 1, 0, 1 };
+ int zswap[4] = { 0, 0, 1, 1 };
+
+ /*
+ * complex 2-by-2 system:
+ *
+ * [ A11 - z A12 ] [ X1 ] = [ B1 ]
+ * [ A21 A22 - z ] [ X2 ] [ B2 ]
+ *
+ * (z complex)
+ *
+ * where the X and B values are complex.
+ */
+
+ civ = (double *) ci;
+ crv = (double *) cr;
+
+ /*
+ * compute the real part of C = A - z*I - use column ordering
+ * here since porting from lapack
+ */
+ cr[0][0] = gsl_matrix_get(A, 0, 0) - GSL_REAL(*z);
+ cr[1][1] = gsl_matrix_get(A, 1, 1) - GSL_REAL(*z);
+ cr[0][1] = gsl_matrix_get(A, 1, 0);
+ cr[1][0] = gsl_matrix_get(A, 0, 1);
+
+ /* compute the imaginary part */
+ ci[0][0] = -GSL_IMAG(*z);
+ ci[0][1] = 0.0;
+ ci[1][0] = 0.0;
+ ci[1][1] = -GSL_IMAG(*z);
+
+ cmax = 0.0;
+ icmax = 0;
+
+ for (j = 0; j < 4; ++j)
+ {
+ if (fabs(crv[j]) + fabs(civ[j]) > cmax)
+ {
+ cmax = fabs(crv[j]) + fabs(civ[j]);
+ icmax = j;
+ }
+ }
+
+ bval1 = gsl_vector_complex_get(b, 0);
+ bval2 = gsl_vector_complex_get(b, 1);
+
+ /* if norm(C) < smin, use smin*I */
+ if (cmax < smin)
+ {
+ bnorm = GSL_MAX(fabs(GSL_REAL(bval1)) + fabs(GSL_IMAG(bval1)),
+ fabs(GSL_REAL(bval2)) + fabs(GSL_IMAG(bval2)));
+ if (smin < 1.0 && bnorm > 1.0)
+ {
+ if (bnorm > GSL_UNSYMMV_BIGNUM*smin)
+ scale = 1.0 / bnorm;
+ }
+
+ temp = scale / smin;
+ xval1 = gsl_complex_mul_real(bval1, temp);
+ xval2 = gsl_complex_mul_real(bval2, temp);
+ gsl_vector_complex_set(x, 0, xval1);
+ gsl_vector_complex_set(x, 1, xval2);
+ *xnorm = temp * bnorm;
+ *s = scale;
+ return;
+ }
+
+ /* gaussian elimination with complete pivoting */
+ ur11 = crv[icmax];
+ ui11 = civ[icmax];
+ cr21 = crv[ipivot[1][icmax]];
+ ci21 = civ[ipivot[1][icmax]];
+ ur12 = crv[ipivot[2][icmax]];
+ ui12 = civ[ipivot[2][icmax]];
+ cr22 = crv[ipivot[3][icmax]];
+ ci22 = civ[ipivot[3][icmax]];
+
+ if (icmax == 0 || icmax == 3)
+ {
+ /* off diagonals of pivoted C are real */
+ if (fabs(ur11) > fabs(ui11))
+ {
+ temp = ui11 / ur11;
+ ur11r = 1.0 / (ur11 * (1.0 + temp*temp));
+ ui11r = -temp * ur11r;
+ }
+ else
+ {
+ temp = ur11 / ui11;
+ ui11r = -1.0 / (ui11 * (1.0 + temp*temp));
+ ur11r = -temp*ui11r;
+ }
+ lr21 = cr21 * ur11r;
+ li21 = cr21 * ui11r;
+ ur12s = ur12 * ur11r;
+ ui12s = ur12 * ui11r;
+ ur22 = cr22 - ur12 * lr21;
+ ui22 = ci22 - ur12 * li21;
+ }
+ else
+ {
+ /* diagonals of pivoted C are real */
+ ur11r = 1.0 / ur11;
+ ui11r = 0.0;
+ lr21 = cr21 * ur11r;
+ li21 = ci21 * ur11r;
+ ur12s = ur12 * ur11r;
+ ui12s = ui12 * ur11r;
+ ur22 = cr22 - ur12 * lr21 + ui12 * li21;
+ ui22 = -ur12 * li21 - ui12 * lr21;
+ }
+
+ u22abs = fabs(ur22) + fabs(ui22);
+
+ /* if smaller pivot < smin, use smin */
+ if (u22abs < smin)
+ {
+ ur22 = smin;
+ ui22 = 0.0;
+ }
+
+ if (rswap[icmax])
+ {
+ br2 = GSL_REAL(bval1);
+ bi2 = GSL_IMAG(bval1);
+ br1 = GSL_REAL(bval2);
+ bi1 = GSL_IMAG(bval2);
+ }
+ else
+ {
+ br1 = GSL_REAL(bval1);
+ bi1 = GSL_IMAG(bval1);
+ br2 = GSL_REAL(bval2);
+ bi2 = GSL_IMAG(bval2);
+ }
+
+ br2 += li21*bi1 - lr21*br1;
+ bi2 -= li21*br1 + lr21*bi1;
+ bbnd = GSL_MAX((fabs(br1) + fabs(bi1)) *
+ (u22abs * (fabs(ur11r) + fabs(ui11r))),
+ fabs(br2) + fabs(bi2));
+ if (bbnd > 1.0 && u22abs < 1.0)
+ {
+ if (bbnd >= GSL_UNSYMMV_BIGNUM*u22abs)
+ {
+ scale = 1.0 / bbnd;
+ br1 *= scale;
+ bi1 *= scale;
+ br2 *= scale;
+ bi2 *= scale;
+ }
+ }
+
+ GSL_SET_COMPLEX(&b1, br2, bi2);
+ GSL_SET_COMPLEX(&b2, ur22, ui22);
+ xval2 = gsl_complex_div(b1, b2);
+
+ xr1 = ur11r*br1 - ui11r*bi1 - ur12s*GSL_REAL(xval2) + ui12s*GSL_IMAG(xval2);
+ xi1 = ui11r*br1 + ur11r*bi1 - ui12s*GSL_REAL(xval2) - ur12s*GSL_IMAG(xval2);
+ GSL_SET_COMPLEX(&xval1, xr1, xi1);
+
+ if (zswap[icmax])
+ {
+ gsl_vector_complex_set(x, 0, xval2);
+ gsl_vector_complex_set(x, 1, xval1);
+ }
+ else
+ {
+ gsl_vector_complex_set(x, 0, xval1);
+ gsl_vector_complex_set(x, 1, xval2);
+ }
+
+ *xnorm = GSL_MAX(fabs(GSL_REAL(xval1)) + fabs(GSL_IMAG(xval1)),
+ fabs(GSL_REAL(xval2)) + fabs(GSL_IMAG(xval2)));
+
+ /* further scaling if norm(A) norm(X) > overflow */
+ if (*xnorm > 1.0 && cmax > 1.0)
+ {
+ if (*xnorm > GSL_UNSYMMV_BIGNUM / cmax)
+ {
+ temp = cmax / GSL_UNSYMMV_BIGNUM;
+ gsl_blas_zdscal(temp, x);
+ *xnorm *= temp;
+ scale *= temp;
+ }
+ }
+ } /* if (N == 2) */
+
+ *s = scale;
+} /* unsymmv_solve_equation_z() */
+
+/*
+unsymmv_normalize_eigenvectors()
+ Normalize eigenvectors so that their Euclidean norm is 1
+
+Inputs: eval - eigenvalues
+ evec - eigenvectors
+*/
+
+static void
+unsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
+ gsl_matrix_complex *evec)
+{
+ const size_t N = evec->size1;
+ size_t i; /* looping */
+ gsl_complex ei;
+ gsl_vector_complex_view vi;
+ gsl_vector_view re, im;
+ double scale; /* scaling factor */
+
+ for (i = 0; i < N; ++i)
+ {
+ ei = gsl_vector_complex_get(eval, i);
+ vi = gsl_matrix_complex_column(evec, i);
+
+ re = gsl_vector_complex_real(&vi.vector);
+
+ if (GSL_IMAG(ei) == 0.0)
+ {
+ scale = 1.0 / gsl_blas_dnrm2(&re.vector);
+ gsl_blas_dscal(scale, &re.vector);
+ }
+ else if (GSL_IMAG(ei) > 0.0)
+ {
+ im = gsl_vector_complex_imag(&vi.vector);
+
+ scale = 1.0 / gsl_hypot(gsl_blas_dnrm2(&re.vector),
+ gsl_blas_dnrm2(&im.vector));
+ gsl_blas_zdscal(scale, &vi.vector);
+
+ vi = gsl_matrix_complex_column(evec, i + 1);
+ gsl_blas_zdscal(scale, &vi.vector);
+ }
+ }
+} /* unsymmv_normalize_eigenvectors() */
Index: linalg/balancemat.c
===================================================================
RCS file: /cvs/gsl/gsl/linalg/balancemat.c,v
retrieving revision 1.2
diff -u -r1.2 balancemat.c
--- linalg/balancemat.c 18 Sep 2006 21:26:59 -0000 1.2
+++ linalg/balancemat.c 21 Sep 2006 18:26:30 -0000
@@ -41,101 +41,106 @@
#include <gsl/gsl_linalg.h>
+#define FLOAT_RADIX 2.0
+#define FLOAT_RADIX_SQ (FLOAT_RADIX * FLOAT_RADIX)
+
int
-gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D)
+gsl_linalg_balance_matrix(gsl_matrix * A, gsl_vector * D)
{
- double row_norm, col_norm;
- int not_converged;
- const size_t M = A->size1;
- const size_t N = A->size2;
+ const size_t N = A->size1;
- if (M != N)
+ if (N != D->size)
{
- GSL_ERROR ("matrix for balancing must be square", GSL_EINVAL);
+ GSL_ERROR ("vector must match matrix size", GSL_EBADLEN);
}
-
- if (D->size != N)
+ else
{
- GSL_ERROR ("length of D must match dimension of A", GSL_EINVAL);
- }
+ double row_norm,
+ col_norm;
+ int not_converged;
+ gsl_vector_view v;
- gsl_vector_set_all (D, 1.0);
+ /* initialize D to the identity matrix */
+ gsl_vector_set_all(D, 1.0);
- not_converged = 1;
-
- while (not_converged)
- {
- size_t i, j;
- double g, f, s;
+ not_converged = 1;
- not_converged = 0;
-
- for (i = 0; i < N; i++)
+ while (not_converged)
{
- row_norm = 0.0;
- col_norm = 0.0;
+ size_t i, j;
+ double g, f, s;
+
+ not_converged = 0;
- for (j = 0; j < N; ++j)
+ for (i = 0; i < N; ++i)
{
- if (j != i)
+ row_norm = 0.0;
+ col_norm = 0.0;
+
+ for (j = 0; j < N; ++j)
{
- col_norm += fabs (gsl_matrix_get (A, j, i));
- row_norm += fabs (gsl_matrix_get (A, i, j));
+ if (j != i)
+ {
+ col_norm += fabs(gsl_matrix_get(A, j, i));
+ row_norm += fabs(gsl_matrix_get(A, i, j));
+ }
}
- }
- if ((col_norm == 0.0) || (row_norm == 0.0)
- || !gsl_finite (col_norm) || !gsl_finite (row_norm))
- {
- continue;
- }
+ if ((col_norm == 0.0) || (row_norm == 0.0))
+ {
+ continue;
+ }
- g = row_norm / 2.0;
- f = 1.0;
- s = col_norm + row_norm;
-
- /* FIXME: we could use frexp() here */
-
- /*
- * find the integer power of the machine radix which
- * comes closest to balancing the matrix
- */
- while (col_norm < g)
- {
- f *= 2.0;
- col_norm *= 4.0;
- }
+ g = row_norm / FLOAT_RADIX;
+ f = 1.0;
+ s = col_norm + row_norm;
+
+ /*
+ * find the integer power of the machine radix which
+ * comes closest to balancing the matrix
+ */
+ while (col_norm < g)
+ {
+ f *= FLOAT_RADIX;
+ col_norm *= FLOAT_RADIX_SQ;
+ }
- g = row_norm * 2.0;
+ g = row_norm * FLOAT_RADIX;
- while (col_norm > g)
- {
- f /= 2.0;
- col_norm /= 4.0;
- }
+ while (col_norm > g)
+ {
+ f /= FLOAT_RADIX;
+ col_norm /= FLOAT_RADIX_SQ;
+ }
- if ((row_norm + col_norm) < 0.95 * s * f)
- {
- not_converged = 1;
+ if ((row_norm + col_norm) < 0.95 * s * f)
+ {
+ not_converged = 1;
- gsl_vector_set (D, i, f);
+ g = 1.0 / f;
- /* apply similarity transformation */
+ /*
+ * apply similarity transformation D, where
+ * D_{ij} = f_i * delta_{ij}
+ */
+
+ /* multiply by D^{-1} on the left */
+ v = gsl_matrix_row(A, i);
+ gsl_blas_dscal(g, &v.vector);
+
+ /* multiply by D on the right */
+ v = gsl_matrix_column(A, i);
+ gsl_blas_dscal(f, &v.vector);
- if (f != 1.0)
- {
- gsl_vector_view A_row_i = gsl_matrix_row (A, i);
- gsl_vector_view A_col_i = gsl_matrix_column (A, i);
- g = 1.0 / f;
- gsl_blas_dscal (g, &A_row_i.vector);
- gsl_blas_dscal (f, &A_col_i.vector);
+ /* keep track of transformation */
+ gsl_vector_set(D, i, gsl_vector_get(D, i) * f);
}
}
}
- }
- return GSL_SUCCESS;
-}
+ return GSL_SUCCESS;
+ }
+} /* gsl_linalg_balance_matrix() */
/*
gsl_linalg_balance_accum()
[-- Attachment #3: testunsymm.c --]
[-- Type: text/plain, Size: 17837 bytes --]
/*
* testunsymm.c
* Patrick Alken
*
* Compile: gcc -g -O2 -Wall -lm -lgsl -llapack -lf77blas -lcblas -latlas -lg2c
*
* Usage: testunsymm [options]
*
* -i : incremental matrices
* -b : balance the matrices
* -z : compute Schur vectors and test them
* -n size : size of matrices
* -l lower-bound : lower bound for matrix elements
* -u upper-bound : upper bound for matrix elements
* -c num : number of matrices to solve
*/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <getopt.h>
#include <sys/times.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_test.h>
typedef struct
{
gsl_eigen_unsymm_workspace *unsymm_p;
gsl_eigen_unsymmv_workspace *unsymmv_p;
gsl_matrix *A;
gsl_matrix *A2;
gsl_vector_complex *eval;
gsl_vector_complex *evalv;
gsl_matrix_complex *evec;
int compute_z;
gsl_matrix *Z;
gsl_matrix *Zv;
} unsymm_workspace;
unsymm_workspace *unsymm_alloc(size_t n, int compute_z, int do_balance);
void unsymm_free(unsymm_workspace *w);
int unsymm_proc(unsymm_workspace *w);
/*
* Global variables
*/
unsigned long count = 0;
/*
* Prototypes
*/
void make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower,
int upper);
void make_start_matrix(gsl_matrix *m, int lower);
int inc_matrix (gsl_matrix *m, int lower, int upper);
void output_matrix(gsl_matrix *m);
void print_octave(gsl_matrix *m, const char *str);
void print_matrix(gsl_matrix *m, const char *str);
void print_hess(gsl_matrix *H, const char *str);
void print_vector(gsl_vector_complex *eval, const char *str);
int cmp(double a, double b);
int compare(const void *a, const void *b);
void sort_complex_vector(gsl_vector_complex *v);
int test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
gsl_matrix *A, const char *obsname, const char *expname);
int test_Z(gsl_matrix *A, gsl_matrix *Z, gsl_matrix *T);
int test_eigenvectors(gsl_matrix *A, gsl_vector_complex *eval,
gsl_matrix_complex *evec);
void my_error_handler(const char *reason, const char *file, int line,
int err);
/**********************************************
* GSL routines
**********************************************/
unsymm_workspace *
unsymm_alloc(size_t n, int compute_z, int do_balance)
{
unsymm_workspace *w;
w = (unsymm_workspace *) malloc(sizeof(unsymm_workspace));
memset(w, '\0', sizeof(unsymm_workspace));
w->unsymm_p = gsl_eigen_unsymm_alloc(n);
w->unsymmv_p = gsl_eigen_unsymmv_alloc(n);
w->A = gsl_matrix_alloc(n, n);
w->A2 = gsl_matrix_alloc(n, n);
if (compute_z)
{
w->Z = gsl_matrix_alloc(n, n);
w->Zv = gsl_matrix_alloc(n, n);
w->compute_z = 1;
gsl_eigen_unsymm_params(1, do_balance, w->unsymm_p);
}
else
gsl_eigen_unsymm_params(0, do_balance, w->unsymm_p);
w->eval = gsl_vector_complex_alloc(n);
w->evalv = gsl_vector_complex_alloc(n);
w->evec = gsl_matrix_complex_alloc(n, n);
return (w);
} /* unsymm_alloc() */
void
unsymm_free(unsymm_workspace *w)
{
if (!w)
return;
if (w->unsymm_p)
gsl_eigen_unsymm_free(w->unsymm_p);
if (w->unsymmv_p)
gsl_eigen_unsymmv_free(w->unsymmv_p);
if (w->A)
gsl_matrix_free(w->A);
if (w->A2)
gsl_matrix_free(w->A2);
if (w->Z)
gsl_matrix_free(w->Z);
if (w->Zv)
gsl_matrix_free(w->Zv);
if (w->eval)
gsl_vector_complex_free(w->eval);
if (w->evalv)
gsl_vector_complex_free(w->evalv);
if (w->evec)
gsl_matrix_complex_free(w->evec);
free(w);
} /* unsymm_free() */
int
unsymm_proc(unsymm_workspace *w)
{
int s1, s2, s;
if (w->compute_z)
{
s1 = gsl_eigen_unsymm_Z(w->A, w->eval, w->Z, w->unsymm_p);
s2 = gsl_eigen_unsymmv_Z(w->A2, w->evalv, w->evec, w->Zv, w->unsymmv_p);
}
else
{
s1 = gsl_eigen_unsymm(w->A, w->eval, w->unsymm_p);
s2 = gsl_eigen_unsymmv(w->A2, w->evalv, w->evec, w->unsymmv_p);
}
s = s1;
if (s > s2)
s = s2;
return s;
} /* unsymm_proc() */
/**********************************************
* General routines
**********************************************/
void
make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)
{
size_t i, j;
size_t N = m->size1;
for (i = 0; i < N; ++i)
{
for (j = 0; j < N; ++j)
{
gsl_matrix_set(m,
i,
j,
gsl_rng_uniform(r) * (upper - lower) + lower);
}
}
} /* make_random_matrix() */
void
make_start_matrix(gsl_matrix *m, int lower)
{
size_t i, j;
size_t N = m->size1;
for (i = 0; i < N; ++i)
for (j = 0; j < N; ++j)
gsl_matrix_set(m, i, j, (double)lower);
} /* make_start_matrix() */
int
inc_matrix (gsl_matrix *m, int lower, int upper)
{
size_t i = 0;
size_t N = m->size1 * m->size2;
int carry = 1;
for (i = 0; carry > 0 && i < N; i++)
{
double v = m->data[i] + carry;
carry = (v > upper) ? 1 : 0;
m->data[i] = (v > upper) ? lower : v;
}
return carry;
} /* inc_matrix() */
void
output_matrix(gsl_matrix *m)
{
size_t i, j;
size_t N = m->size1;
size_t M = m->size2;
for (i = 0; i < N; ++i)
{
for (j = 0; j < M; ++j)
{
/*printf("%10.9f%s",*/
printf("%10.18e%s",
gsl_matrix_get(m, i, j),
(j < M - 1) ? "," : ";\n");
}
}
}
void
print_octave(gsl_matrix *m, const char *str)
{
FILE *fp;
size_t i, j;
const size_t N = m->size1;
const size_t M = m->size2;
fp = fopen(str, "w");
if (!fp)
return;
fprintf(fp, "# Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
fprintf(fp, "# name: %s\n", str);
fprintf(fp, "# type: matrix\n");
fprintf(fp, "# rows: %u\n", N);
fprintf(fp, "# columns: %u\n", N);
for (i = 0; i < N; ++i)
{
for (j = 0; j < M; ++j)
{
fprintf(fp,
"%10.9f%s",
gsl_matrix_get(m, i, j),
(j < M - 1) ? " " : "\n");
}
}
fclose(fp);
}
void
print_matrix(gsl_matrix *m, const char *str)
{
size_t i, j;
size_t N = m->size1;
size_t M = m->size2;
gsl_matrix_view v;
size_t rows, cols;
size_t r, c;
char buf[100];
/*print_octave(m, str);
return;*/
/*rows = GSL_MIN(15, N);*/
rows = N;
cols = GSL_MIN(15, N);
/*cols = N;*/
for (i = 0; i < N; i += rows)
{
for (j = 0; j < M; j += cols)
{
r = GSL_MIN(rows, N - i);
c = GSL_MIN(cols, N - j);
v = gsl_matrix_submatrix(m,
i,
j,
r,
c);
sprintf(buf, "%s(%u:%u,%u:%u)",
str,
i + 1,
i + r,
j + 1,
j + c);
printf("%s = [\n", buf);
output_matrix(&v.matrix);
printf("]\n");
}
}
} /* print_matrix() */
void
print_vector(gsl_vector_complex *eval, const char *str)
{
size_t N = eval->size;
size_t i;
gsl_complex z;
printf("%s = [\n", str);
for (i = 0; i < N; ++i)
{
z = gsl_vector_complex_get(eval, i);
printf("%.18e %.18e;\n", GSL_REAL(z), GSL_IMAG(z));
}
printf("\n]\n");
} /* print_vector() */
int
cmp(double a, double b)
{
return ((a > b) ? 1 : ((a < b) ? -1 : 0));
} /* cmp() */
int
compare(const void *a, const void *b)
{
const double *x = a;
const double *y = b;
int r1 = cmp(y[0], x[0]);
int r2 = cmp(y[1], x[1]);
if (fabs(x[0] - y[0]) < 1.0e-8)
{
/* real parts are very close to each other */
return r2;
}
else
{
return r1 ? r1 : r2;
}
} /* compare() */
void
sort_complex_vector(gsl_vector_complex *v)
{
qsort(v->data, v->size, 2 * sizeof(double), &compare);
} /* sort_complex_vector() */
int
test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
gsl_matrix *A, const char *obsname, const char *expname)
{
size_t N = expected->size;
size_t i, k;
double max, max_abserr, max_relerr;
max = 0.0;
max_abserr = 0.0;
max_relerr = 0.0;
k = 0;
for (i = 0; i < N; ++i)
{
gsl_complex z = gsl_vector_complex_get(expected, i);
max = GSL_MAX_DBL(max, gsl_complex_abs(z));
}
for (i = 0; i < N; ++i)
{
gsl_complex z_obs = gsl_vector_complex_get(obs, i);
gsl_complex z_exp = gsl_vector_complex_get(expected, i);
double x_obs = GSL_REAL(z_obs);
double y_obs = GSL_IMAG(z_obs);
double x_exp = GSL_REAL(z_exp);
double y_exp = GSL_IMAG(z_exp);
double abserr_x = fabs(x_obs - x_exp);
double abserr_y = fabs(y_obs - y_exp);
double noise = max * GSL_DBL_EPSILON * N * N;
max_abserr = GSL_MAX_DBL(max_abserr, abserr_x + abserr_y);
if (abserr_x < noise && abserr_y < noise)
continue;
if (abserr_x > 1.0e-6 || abserr_y > 1.0e-6)
++k;
}
if (k)
{
printf("==== CASE %lu ===========================\n\n", count);
print_matrix(A, "A");
printf("=== eval - %s ===\n", expname);
print_vector(expected, expname);
printf("=== eval - %s ===\n", obsname);
print_vector(obs, obsname);
printf("max abserr = %g max relerr = %g\n", max_abserr, max_relerr);
printf("=========================================\n\n");
}
return k;
} /* test_evals() */
/* test if A = ZTZ^t (or AZ = ZT) */
int
test_Z(gsl_matrix *A, gsl_matrix *Z, gsl_matrix *T)
{
size_t N = A->size1;
size_t i, j, k;
double rhs, lhs;
double abserr;
gsl_matrix *T1, *T2;
T1 = gsl_matrix_alloc(N, N);
T2 = gsl_matrix_alloc(N, N);
/* zero the lower triangle of T */
if (N > 3)
{
for (i = 2; i < N; ++i)
{
for (j = 0; j < (i - 1); ++j)
{
gsl_matrix_set(T, i, j, 0.0);
}
}
}
else if (N == 3)
{
gsl_matrix_set(T, 2, 0, 0.0);
}
/* compute T1 = A Z */
gsl_blas_dgemm(CblasNoTrans,
CblasNoTrans,
1.0,
A,
Z,
0.0,
T1);
/* compute T2 = Z T */
gsl_blas_dgemm(CblasNoTrans,
CblasNoTrans,
1.0,
Z,
T,
0.0,
T2);
k = 0;
for (i = 0; i < N; ++i)
{
for (j = 0; j < N; ++j)
{
lhs = gsl_matrix_get(T1, i, j);
rhs = gsl_matrix_get(T2, i, j);
abserr = fabs(lhs - rhs);
if (abserr > 1.0e-6)
++k;
}
}
if (k)
{
printf("==== CASE %lu ===========================\n\n", count);
print_matrix(A, "A");
printf("=== Schur Form matrix ===\n");
print_matrix(T, "T");
printf("=== Similarity matrix ===\n");
print_matrix(Z, "Z");
printf("=== A Z ===\n");
print_matrix(T1, "A Z");
printf("=== Z T ===\n");
print_matrix(T2, "Z T");
printf("=== A Z - Z T ===\n");
gsl_matrix_sub(T1, T2);
print_matrix(T1, "A Z - Z T");
printf("=========================================\n\n");
}
gsl_matrix_free(T1);
gsl_matrix_free(T2);
return k;
} /* test_Z() */
int
test_eigenvectors(gsl_matrix *A, gsl_vector_complex *eval,
gsl_matrix_complex *evec)
{
const size_t N = A->size1;
size_t i, j;
int k, s;
gsl_matrix_complex *m;
gsl_vector_complex *x, *y;
gsl_complex z_one, z_zero;
m = gsl_matrix_complex_alloc(N, N);
y = gsl_vector_complex_alloc(N);
x = gsl_vector_complex_alloc(N);
/* m <- A */
for (i = 0; i < N; ++i)
{
for (j = 0; j < N; ++j)
{
gsl_complex z;
GSL_SET_COMPLEX(&z, gsl_matrix_get(A, i, j), 0.0);
gsl_matrix_complex_set(m, i, j, z);
}
}
GSL_SET_COMPLEX(&z_one, 1.0, 0.0);
GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);
s = 0;
/* check eigenvalues */
for (i = 0; i < N; ++i)
{
gsl_vector_complex_view vi = gsl_matrix_complex_column(evec, i);
gsl_complex ei = gsl_vector_complex_get(eval, i);
double norm = gsl_blas_dznrm2(&vi.vector);
/* check that eigenvector is normalized */
gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON, "case %u, normalized",
count);
/* compute x = lambda * v */
gsl_vector_complex_memcpy(x, &vi.vector);
gsl_blas_zscal(ei, x);
/* compute y = A v */
gsl_blas_zgemv(CblasNoTrans, z_one, m, &vi.vector, z_zero, y);
k = 0;
/* now test if y = x */
for (j = 0; j < N; ++j)
{
gsl_complex z;
double lhs_r, lhs_i;
double rhs_r, rhs_i;
z = gsl_vector_complex_get(y, j);
lhs_r = GSL_REAL(z);
lhs_i = GSL_IMAG(z);
z = gsl_vector_complex_get(x, j);
rhs_r = GSL_REAL(z);
rhs_i = GSL_IMAG(z);
if (fabs(lhs_r - rhs_r) > 1e8 * GSL_DBL_EPSILON)
++k;
if (fabs(lhs_i - rhs_i) > 1e8 * GSL_DBL_EPSILON)
++k;
}
if (k)
{
s++;
printf("==== CASE %lu ===========================\n\n", count);
print_matrix(A, "A");
printf("eval = %f + %fi\n", GSL_REAL(ei), GSL_IMAG(ei));
print_vector(&vi.vector, "v");
print_vector(y, "Av");
print_vector(x, "ev_v");
printf("=========================================\n\n");
}
}
gsl_matrix_complex_free(m);
gsl_vector_complex_free(y);
gsl_vector_complex_free(x);
return s;
} /* test_eigenvectors() */
void
my_error_handler(const char *reason, const char *file, int line,
int err)
{
printf("[caught: %s:%d: errno=%d %s]\n", file, line, err, reason);
} /* my_error_handler() */
int
main(int argc, char *argv[])
{
unsymm_workspace *unsymm_workspace_p;
size_t N;
int c;
gsl_matrix *A;
gsl_rng *r;
int incremental; /* incremental/random matrices */
int lower; /* lower bound */
int upper; /* upper bound */
unsigned int nmat; /* number of matrices to solve */
int status;
int compute_z;
int do_balance;
gsl_ieee_env_setup();
gsl_rng_env_setup();
/*gsl_set_error_handler(&my_error_handler);*/
N = 30;
incremental = 0;
compute_z = 0;
do_balance = 0;
lower = -10;
upper = 10;
nmat = 0;
while ((c = getopt(argc, argv, "izbc:n:l:u:t:")) != (-1))
{
switch (c)
{
case 'i':
incremental = 1;
break;
case 'n':
N = strtol(optarg, NULL, 0);
break;
case 'l':
lower = strtol(optarg, NULL, 0);
break;
case 'u':
upper = strtol(optarg, NULL, 0);
break;
case 'c':
nmat = strtoul(optarg, NULL, 0);
break;
case 'b':
do_balance = 1;
break;
case 'z':
compute_z = 1;
break;
case '?':
default:
printf("usage: %s [-i] [-b] [-z] [-n size] [-l lower-bound] [-u upper-bound] [-c num]\n", argv[0]);
exit(1);
break;
} /* switch (c) */
}
A = gsl_matrix_alloc(N, N);
unsymm_workspace_p = unsymm_alloc(N, compute_z, do_balance);
if (!incremental)
r = gsl_rng_alloc(gsl_rng_default);
else
{
r = 0;
make_start_matrix(A, lower);
}
fprintf(stderr, "testing N = %d", N);
if (incremental)
fprintf(stderr, " incrementally");
else
fprintf(stderr, " randomly");
fprintf(stderr, " on element range [%d, %d]", lower, upper);
if (compute_z)
fprintf(stderr, ", with Schur vectors");
if (do_balance)
fprintf(stderr, ", with balancing");
fprintf(stderr, "\n");
while (1)
{
if (nmat && (count >= nmat))
break;
++count;
if (!incremental)
make_random_matrix(A, r, lower, upper);
else
{
status = inc_matrix(A, lower, upper);
if (status)
break; /* all done */
}
/*if (count != 186979)
continue;*/
/* make copies of the matrix */
gsl_matrix_memcpy(unsymm_workspace_p->A, A);
gsl_matrix_memcpy(unsymm_workspace_p->A2, A);
/* compute eigenvalues with GSL code */
status = unsymm_proc(unsymm_workspace_p);
if (status != (int) N)
{
printf("=========== CASE %lu ============\n", count);
printf("Failed to converge: found %d eigenvalues\n", status);
print_matrix(A, "A");
}
/*gsl_eigen_unsymmv_sort(unsymm_workspace_p->evalv,
unsymm_workspace_p->evec,
GSL_EIGEN_SORT_ABS_ASC);*/
status = test_eigenvectors(A,
unsymm_workspace_p->evalv,
unsymm_workspace_p->evec);
sort_complex_vector(unsymm_workspace_p->eval);
sort_complex_vector(unsymm_workspace_p->evalv);
status = test_evals(unsymm_workspace_p->eval,
unsymm_workspace_p->evalv,
A,
"unsymm",
"unsymmv");
if (compute_z)
{
test_Z(A, unsymm_workspace_p->Z, unsymm_workspace_p->A);
test_Z(A, unsymm_workspace_p->Zv, unsymm_workspace_p->A2);
}
}
gsl_matrix_free(A);
unsymm_free(unsymm_workspace_p);
if (r)
gsl_rng_free(r);
return 0;
} /* main() */
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: unsymm eigenvector patch
2006-09-21 18:43 unsymm eigenvector patch Patrick Alken
@ 2006-09-22 13:49 ` Brian Gough
0 siblings, 0 replies; 2+ messages in thread
From: Brian Gough @ 2006-09-22 13:49 UTC (permalink / raw)
To: Patrick Alken; +Cc: gsl-discuss
Patrick Alken wrote:
> I've attached a patch which implements eigenvector support for
> real nonsymmetric matrices. The algorithm used is backsubstitution
> to find eigenvectors of the Schur form, followed by backtransformation
> by the Schur vectors. I followed the lapack code closely (dtrevc),
> so this implementation handles degeneracies very well.
Thanks, I have committed the patches to the sourceware cvs. Good work.
--
Brian Gough
Network Theory Ltd,
Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2006-09-22 13:49 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2006-09-21 18:43 unsymm eigenvector patch Patrick Alken
2006-09-22 13:49 ` Brian Gough
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).