This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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 
 

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]