htswanalysis/MACS/lib/gsl/gsl-1.11/eigen/gsl_eigen.h

   1 /* eigen/gsl_eigen.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Brian Gough, Patrick Alken
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #ifndef __GSL_EIGEN_H__
  21 #define __GSL_EIGEN_H__
  22
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_matrix.h>
  25
  26 #undef __BEGIN_DECLS
  27 #undef __END_DECLS
  28 #ifdef __cplusplus
  29 # define __BEGIN_DECLS extern "C" {
  30 # define __END_DECLS }
  31 #else
  32 # define __BEGIN_DECLS /* empty */
  33 # define __END_DECLS /* empty */
  34 #endif
  35
  36 __BEGIN_DECLS
  37
  38 typedef struct {
  39   size_t size;
  40   double * d;
  41   double * sd;
  42 } gsl_eigen_symm_workspace;
  43
  44 gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
  45 void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
  46 int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);
  47
  48 typedef struct {
  49   size_t size;
  50   double * d;
  51   double * sd;
  52   double * gc;
  53   double * gs;
  54 } gsl_eigen_symmv_workspace;
  55
  56 gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
  57 void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
  58 int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);
  59
  60 typedef struct {
  61   size_t size;
  62   double * d;
  63   double * sd;
  64   double * tau;
  65 } gsl_eigen_herm_workspace;
  66
  67 gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
  68 void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
  69 int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
  70                          gsl_eigen_herm_workspace * w);
  71
  72 typedef struct {
  73   size_t size;
  74   double * d;
  75   double * sd;
  76   double * tau;
  77   double * gc;
  78   double * gs;
  79 } gsl_eigen_hermv_workspace;
  80
  81 gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
  82 void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
  83 int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
  84                            gsl_matrix_complex * evec,
  85                            gsl_eigen_hermv_workspace * w);
  86
  87 typedef struct {
  88   size_t size;           /* matrix size */
  89   size_t max_iterations; /* max iterations since last eigenvalue found */
  90   size_t n_iter;         /* number of iterations since last eigenvalue found */
  91   size_t n_evals;        /* number of eigenvalues found so far */
  92
  93   int compute_t;         /* compute Schur form T = Z^t A Z */
  94
  95   gsl_matrix *H;         /* pointer to Hessenberg matrix */
  96   gsl_matrix *Z;         /* pointer to Schur vector matrix */
  97 } gsl_eigen_francis_workspace;
  98
  99 gsl_eigen_francis_workspace * gsl_eigen_francis_alloc (void);
 100 void gsl_eigen_francis_free (gsl_eigen_francis_workspace * w);
 101 void gsl_eigen_francis_T (const int compute_t,
 102                           gsl_eigen_francis_workspace * w);
 103 int gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
 104                        gsl_eigen_francis_workspace * w);
 105 int gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
 106                          gsl_matrix * Z,
 107                          gsl_eigen_francis_workspace * w);
 108
 109 typedef struct {
 110   size_t size;                 /* size of matrices */
 111   gsl_vector *diag;            /* diagonal matrix elements from balancing */
 112   gsl_vector *tau;             /* Householder coefficients */
 113   gsl_matrix *Z;               /* pointer to Z matrix */
 114   int do_balance;              /* perform balancing transformation? */
 115   size_t n_evals;              /* number of eigenvalues found */
 116
 117   gsl_eigen_francis_workspace *francis_workspace_p;
 118 } gsl_eigen_nonsymm_workspace;
 119
 120 gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n);
 121 void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w);
 122 void gsl_eigen_nonsymm_params (const int compute_t, const int balance,
 123                                gsl_eigen_nonsymm_workspace *w);
 124 int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
 125                        gsl_eigen_nonsymm_workspace * w);
 126 int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
 127                          gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w);
 128
 129 typedef struct {
 130   size_t size;                 /* size of matrices */
 131   gsl_vector *work;            /* scratch workspace */
 132   gsl_vector *work2;           /* scratch workspace */
 133   gsl_vector *work3;           /* scratch workspace */
 134
 135   gsl_matrix *Z;               /* pointer to Schur vectors */
 136
 137   gsl_eigen_nonsymm_workspace *nonsymm_workspace_p;
 138 } gsl_eigen_nonsymmv_workspace;
 139
 140 gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n);
 141 void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w);
 142 int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
 143                         gsl_matrix_complex * evec,
 144                         gsl_eigen_nonsymmv_workspace * w);
 145 int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
 146                           gsl_matrix_complex * evec, gsl_matrix * Z,
 147                           gsl_eigen_nonsymmv_workspace * w);
 148
 149 typedef struct {
 150   size_t size;            /* size of matrices */
 151   gsl_eigen_symm_workspace *symm_workspace_p;
 152 } gsl_eigen_gensymm_workspace;
 153
 154 gsl_eigen_gensymm_workspace * gsl_eigen_gensymm_alloc (const size_t n);
 155 void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w);
 156 int gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B,
 157                        gsl_vector * eval, gsl_eigen_gensymm_workspace * w);
 158 int gsl_eigen_gensymm_standardize (gsl_matrix * A, const gsl_matrix * B);
 159
 160 typedef struct {
 161   size_t size;            /* size of matrices */
 162   gsl_eigen_symmv_workspace *symmv_workspace_p;
 163 } gsl_eigen_gensymmv_workspace;
 164
 165 gsl_eigen_gensymmv_workspace * gsl_eigen_gensymmv_alloc (const size_t n);
 166 void gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w);
 167 int gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B,
 168                         gsl_vector * eval, gsl_matrix * evec,
 169                         gsl_eigen_gensymmv_workspace * w);
 170
 171 typedef struct {
 172   size_t size;            /* size of matrices */
 173   gsl_eigen_herm_workspace *herm_workspace_p;
 174 } gsl_eigen_genherm_workspace;
 175
 176 gsl_eigen_genherm_workspace * gsl_eigen_genherm_alloc (const size_t n);
 177 void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w);
 178 int gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
 179                        gsl_vector * eval, gsl_eigen_genherm_workspace * w);
 180 int gsl_eigen_genherm_standardize (gsl_matrix_complex * A,
 181                                    const gsl_matrix_complex * B);
 182
 183 typedef struct {
 184   size_t size;            /* size of matrices */
 185   gsl_eigen_hermv_workspace *hermv_workspace_p;
 186 } gsl_eigen_genhermv_workspace;
 187
 188 gsl_eigen_genhermv_workspace * gsl_eigen_genhermv_alloc (const size_t n);
 189 void gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * w);
 190 int gsl_eigen_genhermv (gsl_matrix_complex * A, gsl_matrix_complex * B,
 191                         gsl_vector * eval, gsl_matrix_complex * evec,
 192                         gsl_eigen_genhermv_workspace * w);
 193
 194 typedef struct {
 195   size_t size;            /* size of matrices */
 196   gsl_vector *work;       /* scratch workspace */
 197
 198   size_t n_evals;         /* number of eigenvalues found */
 199   size_t max_iterations;  /* maximum QZ iterations allowed */
 200   size_t n_iter;          /* number of iterations since last eigenvalue found */
 201   double eshift;          /* exceptional shift counter */
 202
 203   int needtop;            /* need to compute top index? */
 204
 205   double atol;            /* tolerance for splitting A matrix */
 206   double btol;            /* tolerance for splitting B matrix */
 207
 208   double ascale;          /* scaling factor for shifts */
 209   double bscale;          /* scaling factor for shifts */
 210
 211   gsl_matrix *H;          /* pointer to hessenberg matrix */
 212   gsl_matrix *R;          /* pointer to upper triangular matrix */
 213
 214   int compute_s;          /* compute generalized Schur form S */
 215   int compute_t;          /* compute generalized Schur form T */
 216
 217   gsl_matrix *Q;          /* pointer to left Schur vectors */
 218   gsl_matrix *Z;          /* pointer to right Schur vectors */
 219 } gsl_eigen_gen_workspace;
 220
 221 gsl_eigen_gen_workspace * gsl_eigen_gen_alloc (const size_t n);
 222 void gsl_eigen_gen_free (gsl_eigen_gen_workspace * w);
 223 void gsl_eigen_gen_params (const int compute_s, const int compute_t,
 224                            const int balance, gsl_eigen_gen_workspace * w);
 225 int gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B,
 226                    gsl_vector_complex * alpha, gsl_vector * beta,
 227                    gsl_eigen_gen_workspace * w);
 228 int gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
 229                       gsl_vector_complex * alpha, gsl_vector * beta,
 230                       gsl_matrix * Q, gsl_matrix * Z,
 231                       gsl_eigen_gen_workspace * w);
 232
 233 typedef struct {
 234   size_t size;            /* size of matrices */
 235
 236   gsl_vector *work1;      /* 1-norm of columns of A */
 237   gsl_vector *work2;      /* 1-norm of columns of B */
 238   gsl_vector *work3;      /* real part of eigenvector */
 239   gsl_vector *work4;      /* imag part of eigenvector */
 240   gsl_vector *work5;      /* real part of back-transformed eigenvector */
 241   gsl_vector *work6;      /* imag part of back-transformed eigenvector */
 242
 243   gsl_matrix *Q;          /* pointer to left Schur vectors */
 244   gsl_matrix *Z;          /* pointer to right Schur vectors */
 245
 246   gsl_eigen_gen_workspace *gen_workspace_p;
 247 } gsl_eigen_genv_workspace;
 248
 249 gsl_eigen_genv_workspace * gsl_eigen_genv_alloc (const size_t n);
 250 void gsl_eigen_genv_free (gsl_eigen_genv_workspace * w);
 251 int gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B,
 252                     gsl_vector_complex * alpha, gsl_vector * beta,
 253                     gsl_matrix_complex * evec,
 254                     gsl_eigen_genv_workspace * w);
 255 int gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
 256                        gsl_vector_complex * alpha, gsl_vector * beta,
 257                        gsl_matrix_complex * evec,
 258                        gsl_matrix * Q, gsl_matrix * Z,
 259                        gsl_eigen_genv_workspace * w);
 260
 261
 262
 263 typedef enum {
 264   GSL_EIGEN_SORT_VAL_ASC,
 265   GSL_EIGEN_SORT_VAL_DESC,
 266   GSL_EIGEN_SORT_ABS_ASC,
 267   GSL_EIGEN_SORT_ABS_DESC
 268 }
 269 gsl_eigen_sort_t;
 270
 271 /* Sort eigensystem results based on eigenvalues.
 272  * Sorts in order of increasing value or increasing
 273  * absolute value.
 274  *
 275  * exceptions: GSL_EBADLEN
 276  */
 277
 278 int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
 279                          gsl_eigen_sort_t sort_type);
 280
 281 int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
 282                          gsl_eigen_sort_t sort_type);
 283
 284 int gsl_eigen_nonsymmv_sort(gsl_vector_complex * eval,
 285                             gsl_matrix_complex * evec,
 286                             gsl_eigen_sort_t sort_type);
 287
 288 int gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec,
 289                              gsl_eigen_sort_t sort_type);
 290
 291 int gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
 292                              gsl_eigen_sort_t sort_type);
 293
 294 int gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
 295                          gsl_matrix_complex * evec,
 296                          gsl_eigen_sort_t sort_type);
 297
 298 /* Prototypes for the schur module */
 299
 300 int gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B,
 301                           double *wr1, double *wr2, double *wi,
 302                           double *scale1, double *scale2);
 303
 304 int gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
 305                              double d1, double d2, const gsl_vector *b,
 306                              gsl_vector *x, double *s, double *xnorm,
 307                              double smin);
 308
 309 int gsl_schur_solve_equation_z(double ca, const gsl_matrix *A,
 310                                gsl_complex *z, double d1, double d2,
 311                                const gsl_vector_complex *b,
 312                                gsl_vector_complex *x, double *s,
 313                                double *xnorm, double smin);
 314
 315
 316 /* The following functions are obsolete: */
 317
 318 /* Eigensolve by Jacobi Method
 319  *
 320  * The data in the matrix input is destroyed.
 321  *
 322  * exceptions:
 323  */
 324 int
 325 gsl_eigen_jacobi(gsl_matrix * matrix,
 326                       gsl_vector * eval,
 327                       gsl_matrix * evec,
 328                       unsigned int max_rot,
 329                       unsigned int * nrot);
 330
 331
 332 /* Invert by Jacobi Method
 333  *
 334  * exceptions:
 335  */
 336 int
 337 gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
 338                              gsl_matrix * ainv,
 339                              unsigned int max_rot);
 340
 341
 342
 343 __END_DECLS
 344
 345 #endif /* __GSL_EIGEN_H__ */