Bug fixed for GSL library

1 files changed, 347 insertions, 0 deletions

diff --git a/thirdparty/includes/GSL/gsl/gsl_eigen.h b/thirdparty/includes/GSL/gsl/gsl_eigen.h
new file mode 100644
index 0000000..6de8dba
--- /dev/null
+++ b/thirdparty/includes/GSL/gsl/gsl_eigen.h
@@ -0,0 +1,347 @@
+/* eigen/gsl_eigen.h
+ * 
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Brian Gough, 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 3 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_EIGEN_H__
+#define __GSL_EIGEN_H__
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+__BEGIN_DECLS
+
+typedef struct {
+  size_t size;
+  double * d;
+  double * sd;
+} gsl_eigen_symm_workspace;
+
+gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
+void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
+int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);
+
+typedef struct {
+  size_t size;
+  double * d;
+  double * sd;
+  double * gc;
+  double * gs;
+} gsl_eigen_symmv_workspace;
+
+gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
+void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
+int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);
+
+typedef struct {
+  size_t size;
+  double * d;
+  double * sd;
+  double * tau;
+} gsl_eigen_herm_workspace;
+
+gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
+void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
+int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
+                         gsl_eigen_herm_workspace * w);
+
+typedef struct {
+  size_t size;
+  double * d;
+  double * sd;
+  double * tau;
+  double * gc;
+  double * gs;
+} gsl_eigen_hermv_workspace;
+
+gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
+void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
+int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval, 
+                           gsl_matrix_complex * evec,
+                           gsl_eigen_hermv_workspace * w);
+
+typedef struct {
+  size_t size;           /* matrix size */
+  size_t max_iterations; /* max iterations since last eigenvalue found */
+  size_t n_iter;         /* number of iterations since last eigenvalue found */
+  size_t n_evals;        /* number of eigenvalues found so far */
+
+  int compute_t;         /* compute Schur form T = Z^t A Z */
+
+  gsl_matrix *H;         /* pointer to Hessenberg matrix */
+  gsl_matrix *Z;         /* pointer to Schur vector matrix */
+} gsl_eigen_francis_workspace;
+
+gsl_eigen_francis_workspace * gsl_eigen_francis_alloc (void);
+void gsl_eigen_francis_free (gsl_eigen_francis_workspace * w);
+void gsl_eigen_francis_T (const int compute_t,
+                          gsl_eigen_francis_workspace * w);
+int gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
+                       gsl_eigen_francis_workspace * w);
+int gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
+                         gsl_matrix * Z,
+                         gsl_eigen_francis_workspace * w);
+
+typedef struct {
+  size_t size;                 /* size of matrices */
+  gsl_vector *diag;            /* diagonal matrix elements from balancing */
+  gsl_vector *tau;             /* Householder coefficients */
+  gsl_matrix *Z;               /* pointer to Z matrix */
+  int do_balance;              /* perform balancing transformation? */
+  size_t n_evals;              /* number of eigenvalues found */
+
+  gsl_eigen_francis_workspace *francis_workspace_p;
+} gsl_eigen_nonsymm_workspace;
+
+gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n);
+void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w);
+void gsl_eigen_nonsymm_params (const int compute_t, const int balance,
+                               gsl_eigen_nonsymm_workspace *w);
+int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
+                       gsl_eigen_nonsymm_workspace * w);
+int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
+                         gsl_matrix * Z, gsl_eigen_nonsymm_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_matrix *Z;               /* pointer to Schur vectors */
+
+  gsl_eigen_nonsymm_workspace *nonsymm_workspace_p;
+} gsl_eigen_nonsymmv_workspace;
+
+gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n);
+void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w);
+void gsl_eigen_nonsymmv_params (const int balance,
+                                gsl_eigen_nonsymmv_workspace *w);
+int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
+                        gsl_matrix_complex * evec,
+                        gsl_eigen_nonsymmv_workspace * w);
+int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
+                          gsl_matrix_complex * evec, gsl_matrix * Z,
+                          gsl_eigen_nonsymmv_workspace * w);
+
+typedef struct {
+  size_t size;            /* size of matrices */
+  gsl_eigen_symm_workspace *symm_workspace_p;
+} gsl_eigen_gensymm_workspace;
+
+gsl_eigen_gensymm_workspace * gsl_eigen_gensymm_alloc (const size_t n);
+void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w);
+int gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B,
+                       gsl_vector * eval, gsl_eigen_gensymm_workspace * w);
+int gsl_eigen_gensymm_standardize (gsl_matrix * A, const gsl_matrix * B);
+
+typedef struct {
+  size_t size;            /* size of matrices */
+  gsl_eigen_symmv_workspace *symmv_workspace_p;
+} gsl_eigen_gensymmv_workspace;
+
+gsl_eigen_gensymmv_workspace * gsl_eigen_gensymmv_alloc (const size_t n);
+void gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w);
+int gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B,
+                        gsl_vector * eval, gsl_matrix * evec,
+                        gsl_eigen_gensymmv_workspace * w);
+
+typedef struct {
+  size_t size;            /* size of matrices */
+  gsl_eigen_herm_workspace *herm_workspace_p;
+} gsl_eigen_genherm_workspace;
+
+gsl_eigen_genherm_workspace * gsl_eigen_genherm_alloc (const size_t n);
+void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w);
+int gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
+                       gsl_vector * eval, gsl_eigen_genherm_workspace * w);
+int gsl_eigen_genherm_standardize (gsl_matrix_complex * A,
+                                   const gsl_matrix_complex * B);
+
+typedef struct {
+  size_t size;            /* size of matrices */
+  gsl_eigen_hermv_workspace *hermv_workspace_p;
+} gsl_eigen_genhermv_workspace;
+
+gsl_eigen_genhermv_workspace * gsl_eigen_genhermv_alloc (const size_t n);
+void gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * w);
+int gsl_eigen_genhermv (gsl_matrix_complex * A, gsl_matrix_complex * B,
+                        gsl_vector * eval, gsl_matrix_complex * evec,
+                        gsl_eigen_genhermv_workspace * w);
+
+typedef struct {
+  size_t size;            /* size of matrices */
+  gsl_vector *work;       /* scratch workspace */
+
+  size_t n_evals;         /* number of eigenvalues found */
+  size_t max_iterations;  /* maximum QZ iterations allowed */
+  size_t n_iter;          /* number of iterations since last eigenvalue found */
+  double eshift;          /* exceptional shift counter */
+
+  int needtop;            /* need to compute top index? */
+
+  double atol;            /* tolerance for splitting A matrix */
+  double btol;            /* tolerance for splitting B matrix */
+
+  double ascale;          /* scaling factor for shifts */
+  double bscale;          /* scaling factor for shifts */
+
+  gsl_matrix *H;          /* pointer to hessenberg matrix */
+  gsl_matrix *R;          /* pointer to upper triangular matrix */
+
+  int compute_s;          /* compute generalized Schur form S */
+  int compute_t;          /* compute generalized Schur form T */
+
+  gsl_matrix *Q;          /* pointer to left Schur vectors */
+  gsl_matrix *Z;          /* pointer to right Schur vectors */
+} gsl_eigen_gen_workspace;
+
+gsl_eigen_gen_workspace * gsl_eigen_gen_alloc (const size_t n);
+void gsl_eigen_gen_free (gsl_eigen_gen_workspace * w);
+void gsl_eigen_gen_params (const int compute_s, const int compute_t,
+                           const int balance, gsl_eigen_gen_workspace * w);
+int gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B,
+                   gsl_vector_complex * alpha, gsl_vector * beta,
+                   gsl_eigen_gen_workspace * w);
+int gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
+                      gsl_vector_complex * alpha, gsl_vector * beta,
+                      gsl_matrix * Q, gsl_matrix * Z,
+                      gsl_eigen_gen_workspace * w);
+
+typedef struct {
+  size_t size;            /* size of matrices */
+
+  gsl_vector *work1;      /* 1-norm of columns of A */
+  gsl_vector *work2;      /* 1-norm of columns of B */
+  gsl_vector *work3;      /* real part of eigenvector */
+  gsl_vector *work4;      /* imag part of eigenvector */
+  gsl_vector *work5;      /* real part of back-transformed eigenvector */
+  gsl_vector *work6;      /* imag part of back-transformed eigenvector */
+
+  gsl_matrix *Q;          /* pointer to left Schur vectors */
+  gsl_matrix *Z;          /* pointer to right Schur vectors */
+
+  gsl_eigen_gen_workspace *gen_workspace_p;
+} gsl_eigen_genv_workspace;
+
+gsl_eigen_genv_workspace * gsl_eigen_genv_alloc (const size_t n);
+void gsl_eigen_genv_free (gsl_eigen_genv_workspace * w);
+int gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B,
+                    gsl_vector_complex * alpha, gsl_vector * beta,
+                    gsl_matrix_complex * evec,
+                    gsl_eigen_genv_workspace * w);
+int gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
+                       gsl_vector_complex * alpha, gsl_vector * beta,
+                       gsl_matrix_complex * evec,
+                       gsl_matrix * Q, gsl_matrix * Z,
+                       gsl_eigen_genv_workspace * w);
+
+
+
+typedef enum {
+  GSL_EIGEN_SORT_VAL_ASC,
+  GSL_EIGEN_SORT_VAL_DESC,
+  GSL_EIGEN_SORT_ABS_ASC,
+  GSL_EIGEN_SORT_ABS_DESC
+}
+gsl_eigen_sort_t;
+
+/* Sort eigensystem results based on eigenvalues.
+ * Sorts in order of increasing value or increasing
+ * absolute value.
+ *
+ * exceptions: GSL_EBADLEN
+ */
+
+int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
+                         gsl_eigen_sort_t sort_type);
+
+int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
+                         gsl_eigen_sort_t sort_type);
+
+int gsl_eigen_nonsymmv_sort(gsl_vector_complex * eval,
+                            gsl_matrix_complex * evec,
+                            gsl_eigen_sort_t sort_type);
+
+int gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec, 
+                             gsl_eigen_sort_t sort_type);
+
+int gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec, 
+                             gsl_eigen_sort_t sort_type);
+
+int gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
+                         gsl_matrix_complex * evec,
+                         gsl_eigen_sort_t sort_type);
+
+/* Prototypes for the schur module */
+
+int gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B,
+                          double *wr1, double *wr2, double *wi,
+                          double *scale1, double *scale2);
+
+int gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
+                             double d1, double d2, const gsl_vector *b,
+                             gsl_vector *x, double *s, double *xnorm,
+                             double smin);
+
+int gsl_schur_solve_equation_z(double ca, const gsl_matrix *A,
+                               gsl_complex *z, double d1, double d2,
+                               const gsl_vector_complex *b,
+                               gsl_vector_complex *x, double *s,
+                               double *xnorm, double smin);
+
+
+/* The following functions are obsolete: */
+
+/* Eigensolve by Jacobi Method
+ *
+ * The data in the matrix input is destroyed.
+ *
+ * exceptions: 
+ */
+int
+gsl_eigen_jacobi(gsl_matrix * matrix,
+                      gsl_vector * eval,
+                      gsl_matrix * evec,
+                      unsigned int max_rot, 
+                      unsigned int * nrot);
+
+
+/* Invert by Jacobi Method
+ *
+ * exceptions: 
+ */
+int
+gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
+                             gsl_matrix * ainv,
+                             unsigned int max_rot);
+
+
+
+__END_DECLS
+
+#endif /* __GSL_EIGEN_H__ */