This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
Re: unsymmetric eigenvalue problem
- From: Patrick dot Alken at colorado dot edu
- To: gsl-discuss at sourceware dot org
- Date: Thu, 18 May 2006 15:36:23 -0600
- Subject: Re: unsymmetric eigenvalue problem
Hi again,
There was a nasty underflow problem in the last patch I submitted, which was a
result of the gsl_linalg_householder_transform routine. In some steps of the
Francis QR algorithm it is necessary to compute a householder vector with
small elements, and the householder_transform wasn't set up to do this. I have
attached a new patch which adds some code to gsl_linalg_householder_transform
to basically scale the input vector v -> v / |v|. This seems to correct the
underflow problem.
On another note, I haven't gotten the convergence criteria completely right
yet for the eigenvalue solver, so I need to study Golub & Van Loan a bit more -
they aren't completely clear about the convergence test. But the patch I
submitted should work a lot better than the previous one.
Thanks,
Patrick Alken
diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.c ./eigen/balance.c
--- /home/palken/tmp2/gsl-1.8/eigen/balance.c 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.c 2006-05-17 11:57:52.000000000 -0600
@@ -0,0 +1,117 @@
+/* eigen/balance.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 "balance.h"
+
+/*
+balance_matrix()
+ Balance a given matrix by applying a diagonal matrix
+similarity transformation so that the rows and columns
+of the new matrix have norms which are the same order of
+magnitude. This is necessary for the unsymmetric eigenvalue
+problem since the calculation can become numerically unstable
+for unbalanced matrices.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 7.5.7
+
+and
+
+Numerical Recipes in C, Section 11.5
+*/
+
+void
+balance_matrix(gsl_matrix * A)
+{
+ double row_norm,
+ col_norm;
+ int not_converged;
+ const size_t N = A->size1;
+
+ not_converged = 1;
+
+ while (not_converged)
+ {
+ size_t i, j;
+ double g, f, s;
+
+ not_converged = 0;
+
+ for (i = 0; i < N; ++i)
+ {
+ row_norm = 0.0;
+ col_norm = 0.0;
+
+ for (j = 0; j < N; ++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))
+ {
+ continue;
+ }
+
+ 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 * FLOAT_RADIX;
+
+ while (col_norm > g)
+ {
+ f /= FLOAT_RADIX;
+ col_norm /= FLOAT_RADIX_SQ;
+ }
+
+ if ((row_norm + col_norm) < 0.95 * s * f)
+ {
+ not_converged = 1;
+
+ g = 1.0 / f;
+
+ /* apply similarity transformation */
+ for (j = 0; j < N; ++j)
+ {
+ gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) * g);
+ }
+ for (j = 0; j < N; ++j)
+ {
+ gsl_matrix_set(A, j, i, gsl_matrix_get(A, j, i) * f);
+ }
+ }
+ }
+ }
+}
diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.h ./eigen/balance.h
--- /home/palken/tmp2/gsl-1.8/eigen/balance.h 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.h 2006-05-17 10:49:00.000000000 -0600
@@ -0,0 +1,32 @@
+/* eigen/balance.h
+ *
+ * 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.
+ */
+
+#ifndef __GSL_BALANCE_H__
+#define __GSL_BALANCE_H__
+
+#include <gsl/gsl_matrix.h>
+
+/* floating point radix */
+#define FLOAT_RADIX 2.0
+
+#define FLOAT_RADIX_SQ (FLOAT_RADIX * FLOAT_RADIX)
+
+void balance_matrix(gsl_matrix * A);
+
+#endif /* __GSL_BALANCE_H__ */
diff -urN /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h ./eigen/gsl_eigen.h
--- /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h 2005-06-26 07:25:34.000000000 -0600
+++ ./eigen/gsl_eigen.h 2006-05-16 16:38:14.000000000 -0600
@@ -59,6 +59,16 @@
typedef struct {
size_t size;
+ gsl_vector *v; /* temporary vector */
+} gsl_eigen_unsymm_workspace;
+
+gsl_eigen_unsymm_workspace * gsl_eigen_unsymm_alloc (const size_t n);
+void gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w);
+int gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_eigen_unsymm_workspace * w);
+
+typedef struct {
+ size_t size;
double * d;
double * sd;
double * tau;
diff -urN /home/palken/tmp2/gsl-1.8/eigen/Makefile.am ./eigen/Makefile.am
--- /home/palken/tmp2/gsl-1.8/eigen/Makefile.am 2004-09-11 07:45:45.000000000 -0600
+++ ./eigen/Makefile.am 2006-05-18 15:04:48.000000000 -0600
@@ -3,7 +3,7 @@
check_PROGRAMS = test
pkginclude_HEADERS = gsl_eigen.h
-libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c herm.c hermv.c sort.c
+libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c unsymm.c herm.c hermv.c sort.c balance.c
INCLUDES= -I$(top_builddir)
noinst_HEADERS = qrstep.c
diff -urN /home/palken/tmp2/gsl-1.8/eigen/unsymm.c ./eigen/unsymm.c
--- /home/palken/tmp2/gsl-1.8/eigen/unsymm.c 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/unsymm.c 2006-05-18 15:20:58.000000000 -0600
@@ -0,0 +1,400 @@
+/* eigen/unsymm.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_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>
+#include "gsl_eigen.h"
+
+#include "balance.h"
+
+/* maximum iterations before failure is reported */
+#define UNSYMM_MAX_ITERATIONS 100
+
+/*
+ * This module computes the eigenvalues of a real unsymmetric
+ * matrix, using the QR decomposition.
+ *
+ * See Golub & Van Loan, "Matrix Computations" (3rd ed),
+ * algorithm 7.5.2
+ */
+
+inline static void zero_subdiag_small_elements(gsl_matrix * A);
+static inline int francis_qrstep(gsl_matrix * H,
+ gsl_eigen_unsymm_workspace * w);
+static void get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
+ gsl_complex * e2);
+
+/*
+gsl_eigen_unsymm_alloc()
+
+Allocate a workspace for solving the unsymmetric
+eigenvalue/eigenvector problem
+
+Inputs: n - size of matrix
+
+Return: pointer to workspace
+*/
+
+gsl_eigen_unsymm_workspace *
+gsl_eigen_unsymm_alloc(const size_t n)
+{
+ gsl_eigen_unsymm_workspace *w;
+
+ if (n == 0)
+ {
+ GSL_ERROR_NULL ("matrix dimension must be positive integer",
+ GSL_EINVAL);
+ }
+
+ w = ((gsl_eigen_unsymm_workspace *)
+ malloc (sizeof (gsl_eigen_unsymm_workspace)));
+
+ if (w == 0)
+ {
+ GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
+ }
+
+ w->size = n;
+ w->v = gsl_vector_alloc(3);
+
+ return (w);
+}
+
+void
+gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w)
+{
+ gsl_vector_free(w->v);
+
+ free(w);
+}
+
+/*
+gsl_eigen_unsymm()
+
+Solve the unsymmetric eigenvalue problem
+
+A x = \lambda x
+
+for the eigenvalues \lambda using algorithm 7.5.2 of
+Golub & Van Loan, "Matrix Computations" (3rd ed)
+
+Inputs: A - matrix
+ eval - where to store eigenvalues
+ w - workspace
+
+Notes: On output, A contains the upper block triangular Schur
+ decomposition.
+*/
+
+int
+gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_eigen_unsymm_workspace * w)
+{
+ /* check matrix and vector sizes */
+
+ if (A->size1 != A->size2)
+ {
+ GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
+ }
+ else if (eval->size != A->size1)
+ {
+ GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
+ }
+ else
+ {
+ size_t N;
+ size_t p, q;
+ gsl_complex lambda1, lambda2; /* eigenvalues */
+ gsl_matrix_view m;
+ size_t nev; /* number of eigenvalues found so far */
+ size_t nit; /* number of iterations */
+ size_t i,j,k;
+
+ N = A->size1;
+
+ /* special cases */
+ if (N == 1)
+ {
+ GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(A, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, 0, lambda1);
+ return GSL_SUCCESS;
+ }
+
+ if (N == 2)
+ {
+ /*
+ * The 2x2 case is special since the matrix is already
+ * in upper quasi-triangular form so no Schur decomposition
+ * is necessary
+ */
+ get_2b2_eigenvalues(A, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, 0, lambda1);
+ gsl_vector_complex_set(eval, 1, lambda2);
+ return GSL_SUCCESS;
+ }
+
+ nev = 0;
+
+ /* balance the matrix */
+ balance_matrix(A);
+
+ /* compute the Hessenberg reduction of A */
+ gsl_linalg_hessenberg(A);
+
+ /* set small elements on the subdiagonal to 0 */
+ zero_subdiag_small_elements(A);
+
+ m = gsl_matrix_submatrix(A, 0, 0, N, N);
+
+ nit = 0;
+ while ((N > 2) && (nit++ < UNSYMM_MAX_ITERATIONS))
+ {
+ francis_qrstep(&m.matrix, w);
+ zero_subdiag_small_elements(&m.matrix);
+
+ /* Find the largest q such that A_{q,q-1} is 0 */
+ q = N - 1;
+ while ((q > 0) && (gsl_matrix_get(A, q, q - 1) != 0.0))
+ {
+ --q;
+ }
+
+ p = 0;
+ while ((p < (N - 1)) && (gsl_matrix_get(A, p + 1, p) != 0.0))
+ {
+ ++p;
+ }
+
+ if (q == (N - 1))
+ {
+ /*
+ * the (N, N - 1) element of the matrix is 0 -
+ * m_{NN} is a real eigenvalue
+ */
+ GSL_SET_COMPLEX(&lambda1,
+ gsl_matrix_get(&m.matrix, q, q), 0.0);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+
+ --N;
+ m = gsl_matrix_submatrix(A, 0, 0, N, N);
+ }
+ else if (q == (N - 2))
+ {
+ /*
+ * The bottom right 2x2 block of m is an eigenvalue
+ * system
+ */
+
+ m = gsl_matrix_submatrix(A, q, q, 2, 2);
+ get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ gsl_vector_complex_set(eval, nev++, lambda2);
+
+ N -= 2;
+ m = gsl_matrix_submatrix(A, 0, 0, N, N);
+ }
+ }
+
+ if (nit > UNSYMM_MAX_ITERATIONS)
+ {
+ GSL_ERROR("maximum iterations exceeded", GSL_EMAXITER);
+ }
+
+ if (N == 1)
+ {
+ GSL_SET_COMPLEX(&lambda1,
+ gsl_matrix_get(&m.matrix, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ }
+ else if (N == 2)
+ {
+ get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ gsl_vector_complex_set(eval, nev++, lambda2);
+ }
+
+ printf("nit = %d\n", nit);
+
+ return GSL_SUCCESS;
+ }
+} /* gsl_eigen_unsymm() */
+
+/********************************************
+ * INTERNAL ROUTINES *
+ ********************************************/
+
+/*
+zero_subdiag_small_elements()
+ Sets to zero all elements on the subdiaganal of a matrix A
+which satisfy
+
+|A_{i,i-1}| <= eps * (|A_{i,i}| + |A_{i-1,i-1}|)
+*/
+
+inline static void
+zero_subdiag_small_elements(gsl_matrix * A)
+{
+ const size_t N = A->size1;
+ size_t i;
+ double dpel = gsl_matrix_get(A, 0, 0);
+
+ for (i = 1; i < N; ++i)
+ {
+ double sel = gsl_matrix_get(A, i, i - 1);
+ double del = gsl_matrix_get(A, i, i);
+
+ if (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel)))
+ {
+ gsl_matrix_set(A, i, i - 1, 0.0);
+ }
+
+ dpel = del;
+ }
+}
+
+/*
+francis_qrstep()
+ Perform a Francis QR step.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed),
+algorithm 7.5.1
+
+Inputs: H - unreduced upper Hessenberg matrix
+ w - workspace
+*/
+
+static inline int
+francis_qrstep(gsl_matrix * H, gsl_eigen_unsymm_workspace * w)
+{
+ const size_t N = H->size1;
+ double s, t, x, y, z;
+ size_t i;
+ gsl_vector_view v;
+ gsl_matrix_view m;
+ double tau_i;
+ size_t q, r;
+
+ /* s = a_1 + a_2 = H_{mm} + H_{nn} */
+ s = gsl_matrix_get(H, N - 2, N - 2) + gsl_matrix_get(H, N - 1, N - 1);
+
+ /* t = a_1 * a_2 = H_{mm} * H_{nn} - H_{mn} * H_{nm} */
+ t = gsl_matrix_get(H, N - 2, N - 2) * gsl_matrix_get(H, N - 1, N - 1) -
+ gsl_matrix_get(H, N - 2, N - 1) * gsl_matrix_get(H, N - 1, N - 2);
+ x = gsl_matrix_get(H, 0, 0) * gsl_matrix_get(H, 0, 0) +
+ gsl_matrix_get(H, 0, 1) * gsl_matrix_get(H, 1, 0) -
+ s*gsl_matrix_get(H, 0, 0) + t;
+ y = gsl_matrix_get(H, 1, 0) *
+ (gsl_matrix_get(H, 0, 0) + gsl_matrix_get(H, 1, 1) - s);
+ z = gsl_matrix_get(H, 1, 0) * gsl_matrix_get(H, 2, 1);
+
+ v = gsl_vector_subvector(w->v, 0, 3);
+
+ for (i = 0; i < N - 2; ++i)
+ {
+ gsl_vector_set(&v.vector, 0, x);
+ gsl_vector_set(&v.vector, 1, y);
+ gsl_vector_set(&v.vector, 2, z);
+ tau_i = gsl_linalg_householder_transform(&v.vector);
+
+ if (tau_i != 0.0)
+ {
+
+ /* q = max(1, i - 1) */
+ q = (1 > ((int)i - 1)) ? 0 : (i - 1);
+
+ /* apply left householder matrix (I - tau_i v v') to H */
+ m = gsl_matrix_submatrix(H, i, q, 3, N - q);
+ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+ /* r = min(i + 3, N - 1) */
+ r = ((i + 3) < (N - 1)) ? (i + 3) : (N - 1);
+
+ /* apply right householder matrix (I - tau_i v v') to H */
+ m = gsl_matrix_submatrix(H, 0, i, r + 1, 3);
+ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+
+ }
+
+ x = gsl_matrix_get(H, i + 1, i);
+ y = gsl_matrix_get(H, i + 2, i);
+ if (i < (N - 3))
+ {
+ z = gsl_matrix_get(H, i + 3, i);
+ }
+ }
+
+ v = gsl_vector_subvector(w->v, 0, 2);
+ gsl_vector_set(&v.vector, 0, x);
+ gsl_vector_set(&v.vector, 1, y);
+ tau_i = gsl_linalg_householder_transform(&v.vector);
+
+ m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
+ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+ m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
+ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+
+ return GSL_SUCCESS;
+}
+
+/*
+get_2b2_eigenvalues()
+ Compute the eigenvalues of a 2x2 real matrix
+
+Inputs: A - 2x2 matrix
+ e1 - where to store eigenvalue 1
+ e2 - where to store eigenvalue 2
+*/
+
+static void
+get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
+ gsl_complex * e2)
+{
+ double discr; /* discriminant of characteristic poly */
+ double a, b, c, d; /* matrix values */
+
+ a = gsl_matrix_get(A, 0, 0);
+ b = gsl_matrix_get(A, 0, 1);
+ c = gsl_matrix_get(A, 1, 0);
+ d = gsl_matrix_get(A, 1, 1);
+
+ discr = (a + d)*(a + d) - 4.0*(a*d - b*c);
+ if (discr < 0.0)
+ {
+ GSL_SET_REAL(e1, 0.5*(a + d));
+ GSL_SET_REAL(e2, 0.5*(a + d));
+
+ GSL_SET_IMAG(e1, 0.5*sqrt(-discr));
+ GSL_SET_IMAG(e2, -0.5*sqrt(-discr));
+ }
+ else
+ {
+ GSL_SET_REAL(e1, 0.5*(a + d + sqrt(discr)));
+ GSL_SET_REAL(e2, 0.5*(a + d - sqrt(discr)));
+
+ GSL_SET_IMAG(e1, 0.0);
+ GSL_SET_IMAG(e2, 0.0);
+ }
+}
diff -urN /home/palken/tmp2/gsl-1.8/linalg/gsl_linalg.h ./linalg/gsl_linalg.h
--- /home/palken/tmp2/gsl-1.8/linalg/gsl_linalg.h 2005-11-09 14:13:14.000000000 -0700
+++ ./linalg/gsl_linalg.h 2006-05-16 10:36:55.000000000 -0600
@@ -113,6 +113,10 @@
const gsl_vector_complex * v,
gsl_vector_complex * w);
+/* Hessenberg reduction */
+
+int gsl_linalg_hessenberg (gsl_matrix * A);
+
/* Singular Value Decomposition
* exceptions:
diff -urN /home/palken/tmp2/gsl-1.8/linalg/hessenberg.c ./linalg/hessenberg.c
--- /home/palken/tmp2/gsl-1.8/linalg/hessenberg.c 1969-12-31 17:00:00.000000000 -0700
+++ ./linalg/hessenberg.c 2006-05-16 15:11:12.000000000 -0600
@@ -0,0 +1,76 @@
+/* linalg/hessenberg.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 <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+
+/*
+gsl_linalg_hessenberg()
+ Compute the Householder reduction to Hessenberg form of a
+square matrix A.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm
+7.4.2
+*/
+
+int
+gsl_linalg_hessenberg(gsl_matrix * A)
+{
+ if (A->size1 != A->size2)
+ {
+ GSL_ERROR ("Hessenberg reduction requires square matrix",
+ GSL_ENOTSQR);
+ }
+ else
+ {
+ const size_t N = A->size1;
+ size_t i, j;
+ gsl_vector_view v;
+ gsl_vector *v_copy;
+ gsl_matrix_view m;
+ double tau_i; /* beta in algorithm 7.4.2 */
+
+ v_copy = gsl_vector_alloc(N);
+
+ for (i = 0; i < N - 2; ++i)
+ {
+ /* make a copy of A(i + 1:n, i) */
+
+ v = gsl_matrix_column(A, i);
+ gsl_vector_memcpy(v_copy, &v.vector);
+ v = gsl_vector_subvector(v_copy, i + 1, N - (i + 1));
+
+ /* compute householder transformation of A(i+1:n,i) */
+ tau_i = gsl_linalg_householder_transform(&v.vector);
+
+ /* apply left householder matrix (I - tau_i v v') to A */
+ m = gsl_matrix_submatrix(A, i + 1, i, N - (i + 1), N - i);
+ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+ /* apply right householder matrix (I - tau_i v v') to A */
+ m = gsl_matrix_submatrix(A, 0, i + 1, N, N - (i + 1));
+ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+ }
+
+ gsl_vector_free(v_copy);
+
+ return GSL_SUCCESS;
+ }
+}
diff -urN /home/palken/tmp2/gsl-1.8/linalg/householder.c ./linalg/householder.c
--- /home/palken/tmp2/gsl-1.8/linalg/householder.c 2005-06-26 07:25:35.000000000 -0600
+++ ./linalg/householder.c 2006-05-18 14:56:00.000000000 -0600
@@ -41,10 +41,25 @@
else
{
double alpha, beta, tau ;
+ size_t i;
+ double xnorm;
gsl_vector_view x = gsl_vector_subvector (v, 1, n - 1) ;
- double xnorm = gsl_blas_dnrm2 (&x.vector);
+ double vnorm = gsl_blas_dnrm2 (v);
+
+ if (vnorm == 0)
+ {
+ return 0.0;
+ }
+
+ /* v -> v / |v| to prevent overflow/underflow problems */
+ for (i = 0; i < n; ++i)
+ {
+ gsl_vector_set(v, i, gsl_vector_get(v, i) / vnorm);
+ }
+
+ xnorm = gsl_blas_dnrm2 (&x.vector);
if (xnorm == 0)
{
diff -urN /home/palken/tmp2/gsl-1.8/linalg/Makefile.am ./linalg/Makefile.am
--- /home/palken/tmp2/gsl-1.8/linalg/Makefile.am 2004-09-13 12:17:04.000000000 -0600
+++ ./linalg/Makefile.am 2006-05-18 15:06:29.000000000 -0600
@@ -4,7 +4,7 @@
INCLUDES= -I$(top_builddir)
-libgsllinalg_la_SOURCES = multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c qr.c qrpt.c lq.c ptlq.c svd.c householder.c householdercomplex.c cholesky.c symmtd.c hermtd.c bidiag.c balance.c
+libgsllinalg_la_SOURCES = multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c qr.c qrpt.c lq.c ptlq.c svd.c householder.c householdercomplex.c hessenberg.c cholesky.c symmtd.c hermtd.c bidiag.c balance.c
noinst_HEADERS = givens.c apply_givens.c svdstep.c tridiag.h