1 files changed, 285 insertions, 0 deletions
diff --git a/gsl-1.9/eigen/nonsymm.c b/gsl-1.9/eigen/nonsymm.c
new file mode 100644
index 0000000..de9088a
--- /dev/null
+++ b/gsl-1.9/eigen/nonsymm.c
@@ -0,0 +1,285 @@
+/* eigen/nonsymm.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 <config.h>
+#include <stdlib.h>
+#include <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 of a real nonsymmetric
+ * matrix, using the double shift Francis method.
+ *
+ * See the references in francis.c.
+ *
+ * This module gets the matrix ready by balancing it and
+ * reducing it to Hessenberg form before passing it to the
+ * francis module.
+ */
+
+/*
+gsl_eigen_nonsymm_alloc()
+
+Allocate a workspace for solving the nonsymmetric eigenvalue problem.
+The size of this workspace is O(2n)
+
+Inputs: n - size of matrix
+
+Return: pointer to workspace
+*/
+
+gsl_eigen_nonsymm_workspace *
+gsl_eigen_nonsymm_alloc(const size_t n)
+{
+  gsl_eigen_nonsymm_workspace *w;
+
+  if (n == 0)
+    {
+      GSL_ERROR_NULL ("matrix dimension must be positive integer",
+                      GSL_EINVAL);
+    }
+
+  w = (gsl_eigen_nonsymm_workspace *)
+      malloc (sizeof (gsl_eigen_nonsymm_workspace));
+
+  if (w == 0)
+    {
+      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
+    }
+
+  w->size = n;
+  w->Z = NULL;
+  w->do_balance = 0;
+
+  w->diag = gsl_vector_alloc(n);
+
+  if (w->diag == 0)
+    {
+      GSL_ERROR_NULL ("failed to allocate space for balancing vector", GSL_ENOMEM);
+    }
+
+  w->tau = gsl_vector_alloc(n);
+
+  if (w->tau == 0)
+    {
+      GSL_ERROR_NULL ("failed to allocate space for hessenberg coefficients", GSL_ENOMEM);
+    }
+
+  w->francis_workspace_p = gsl_eigen_francis_alloc();
+
+  if (w->francis_workspace_p == 0)
+    {
+      GSL_ERROR_NULL ("failed to allocate space for francis workspace", GSL_ENOMEM);
+    }
+
+  return (w);
+} /* gsl_eigen_nonsymm_alloc() */
+
+/*
+gsl_eigen_nonsymm_free()
+  Free workspace w
+*/
+
+void
+gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w)
+{
+  gsl_vector_free(w->tau);
+
+  gsl_vector_free(w->diag);
+
+  gsl_eigen_francis_free(w->francis_workspace_p);
+
+  free(w);
+} /* gsl_eigen_nonsymm_free() */
+
+/*
+gsl_eigen_nonsymm_params()
+  Set some parameters which define how we solve the eigenvalue
+problem.
+
+Inputs: compute_t - 1 if we want to compute T, 0 if not
+        balance   - 1 if we want to balance the matrix, 0 if not
+        w         - nonsymm workspace
+*/
+
+void
+gsl_eigen_nonsymm_params (const int compute_t, const int balance,
+                          gsl_eigen_nonsymm_workspace *w)
+{
+  gsl_eigen_francis_T(compute_t, w->francis_workspace_p);
+  w->do_balance = balance;
+} /* gsl_eigen_nonsymm_params() */
+
+/*
+gsl_eigen_nonsymm()
+
+Solve the nonsymmetric eigenvalue problem
+
+A x = \lambda x
+
+for the eigenvalues \lambda using the Francis method.
+
+Here we compute the real Schur form
+
+T = Z^t A Z
+
+with the diagonal blocks of T giving us the eigenvalues.
+Z is a matrix of Schur vectors which is not computed by
+this algorithm. See gsl_eigen_nonsymm_Z().
+
+Inputs: A    - general real matrix
+        eval - where to store eigenvalues
+        w    - workspace
+
+Return: success or error
+
+Notes: If T is computed, it is stored in A on output. Otherwise
+       the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
+       blocks.
+*/
+
+int
+gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
+                   gsl_eigen_nonsymm_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
+    {
+      int s;
+
+      if (w->do_balance)
+        {
+          /* balance the matrix */
+          gsl_linalg_balance_matrix(A, w->diag);
+        }
+
+      /* compute the Hessenberg reduction of A */
+      gsl_linalg_hessenberg(A, w->tau);
+
+      if (w->Z)
+        {
+          /*
+           * initialize the matrix Z to U, which is the matrix used
+           * to construct the Hessenberg reduction.
+           */
+
+          /* compute U and store it in Z */
+          gsl_linalg_hessenberg_unpack(A, w->tau, w->Z);
+
+          /* find the eigenvalues and Schur vectors */
+          s = gsl_eigen_francis_Z(A, eval, w->Z, w->francis_workspace_p);
+
+          if (w->do_balance)
+            {
+              /*
+               * The Schur vectors in Z are the vectors for the balanced
+               * matrix. We now must undo the balancing to get the
+               * vectors for the original matrix A.
+               */
+              gsl_linalg_balance_accum(w->Z, w->diag);
+            }
+        }
+      else
+        {
+          /* find the eigenvalues only */
+          s = gsl_eigen_francis(A, eval, w->francis_workspace_p);
+        }
+
+      w->n_evals = w->francis_workspace_p->n_evals;
+
+      return s;
+    }
+} /* gsl_eigen_nonsymm() */
+
+/*
+gsl_eigen_nonsymm_Z()
+
+Solve the nonsymmetric eigenvalue problem
+
+A x = \lambda x
+
+for the eigenvalues \lambda.
+
+Here we compute the real Schur form
+
+T = Z^t A Z
+
+with the diagonal blocks of T giving us the eigenvalues.
+Z is the matrix of Schur vectors.
+
+Inputs: A    - general real matrix
+        eval - where to store eigenvalues
+        Z    - where to store Schur vectors
+        w    - workspace
+
+Return: success or error
+
+Notes: If T is computed, it is stored in A on output. Otherwise
+       the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
+       blocks.
+*/
+
+int
+gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
+                     gsl_matrix * Z, gsl_eigen_nonsymm_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 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_nonsymm(A, eval, w);
+
+      w->Z = NULL;
+
+      return s;
+    }
+} /* gsl_eigen_nonsymm_Z() */