include/gsl/gsl_eigen.h

0001 /* eigen/gsl_eigen.h
0002  *
0003  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Brian Gough, Patrick Alken
0004  *
0005  * This program is free software; you can redistribute it and/or modify
0006  * it under the terms of the GNU General Public License as published by
0007  * the Free Software Foundation; either version 3 of the License, or (at
0008  * your option) any later version.
0009  *
0010  * This program is distributed in the hope that it will be useful, but
0011  * WITHOUT ANY WARRANTY; without even the implied warranty of
0012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0013  * General Public License for more details.
0014  *
0015  * You should have received a copy of the GNU General Public License
0016  * along with this program; if not, write to the Free Software
0017  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0018  */
0019
0020 #ifndef __GSL_EIGEN_H__
0021 #define __GSL_EIGEN_H__
0022
0023 #include <gsl/gsl_vector.h>
0024 #include <gsl/gsl_matrix.h>
0025
0026 #undef __BEGIN_DECLS
0027 #undef __END_DECLS
0028 #ifdef __cplusplus
0029 # define __BEGIN_DECLS extern "C" {
0030 # define __END_DECLS }
0031 #else
0032 # define __BEGIN_DECLS /* empty */
0033 # define __END_DECLS /* empty */
0034 #endif
0035
0036 __BEGIN_DECLS
0037
0038 typedef struct {
0039   size_t size;
0040   double * d;
0041   double * sd;
0042 } gsl_eigen_symm_workspace;
0043
0044 gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
0045 void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
0046 int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);
0047
0048 typedef struct {
0049   size_t size;
0050   double * d;
0051   double * sd;
0052   double * gc;
0053   double * gs;
0054 } gsl_eigen_symmv_workspace;
0055
0056 gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
0057 void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
0058 int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);
0059
0060 typedef struct {
0061   size_t size;
0062   double * d;
0063   double * sd;
0064   double * tau;
0065 } gsl_eigen_herm_workspace;
0066
0067 gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
0068 void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
0069 int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
0070                          gsl_eigen_herm_workspace * w);
0071
0072 typedef struct {
0073   size_t size;
0074   double * d;
0075   double * sd;
0076   double * tau;
0077   double * gc;
0078   double * gs;
0079 } gsl_eigen_hermv_workspace;
0080
0081 gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
0082 void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
0083 int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
0084                            gsl_matrix_complex * evec,
0085                            gsl_eigen_hermv_workspace * w);
0086
0087 typedef struct {
0088   size_t size;           /* matrix size */
0089   size_t max_iterations; /* max iterations since last eigenvalue found */
0090   size_t n_iter;         /* number of iterations since last eigenvalue found */
0091   size_t n_evals;        /* number of eigenvalues found so far */
0092
0093   int compute_t;         /* compute Schur form T = Z^t A Z */
0094
0095   gsl_matrix *H;         /* pointer to Hessenberg matrix */
0096   gsl_matrix *Z;         /* pointer to Schur vector matrix */
0097 } gsl_eigen_francis_workspace;
0098
0099 gsl_eigen_francis_workspace * gsl_eigen_francis_alloc (void);
0100 void gsl_eigen_francis_free (gsl_eigen_francis_workspace * w);
0101 void gsl_eigen_francis_T (const int compute_t,
0102                           gsl_eigen_francis_workspace * w);
0103 int gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
0104                        gsl_eigen_francis_workspace * w);
0105 int gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
0106                          gsl_matrix * Z,
0107                          gsl_eigen_francis_workspace * w);
0108
0109 typedef struct {
0110   size_t size;                 /* size of matrices */
0111   gsl_vector *diag;            /* diagonal matrix elements from balancing */
0112   gsl_vector *tau;             /* Householder coefficients */
0113   gsl_matrix *Z;               /* pointer to Z matrix */
0114   int do_balance;              /* perform balancing transformation? */
0115   size_t n_evals;              /* number of eigenvalues found */
0116
0117   gsl_eigen_francis_workspace *francis_workspace_p;
0118 } gsl_eigen_nonsymm_workspace;
0119
0120 gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n);
0121 void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w);
0122 void gsl_eigen_nonsymm_params (const int compute_t, const int balance,
0123                                gsl_eigen_nonsymm_workspace *w);
0124 int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
0125                        gsl_eigen_nonsymm_workspace * w);
0126 int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
0127                          gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w);
0128
0129 typedef struct {
0130   size_t size;                 /* size of matrices */
0131   gsl_vector *work;            /* scratch workspace */
0132   gsl_vector *work2;           /* scratch workspace */
0133   gsl_vector *work3;           /* scratch workspace */
0134
0135   gsl_matrix *Z;               /* pointer to Schur vectors */
0136
0137   gsl_eigen_nonsymm_workspace *nonsymm_workspace_p;
0138 } gsl_eigen_nonsymmv_workspace;
0139
0140 gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n);
0141 void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w);
0142 void gsl_eigen_nonsymmv_params (const int balance,
0143                                 gsl_eigen_nonsymmv_workspace *w);
0144 int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
0145                         gsl_matrix_complex * evec,
0146                         gsl_eigen_nonsymmv_workspace * w);
0147 int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
0148                           gsl_matrix_complex * evec, gsl_matrix * Z,
0149                           gsl_eigen_nonsymmv_workspace * w);
0150
0151 typedef struct {
0152   size_t size;            /* size of matrices */
0153   gsl_eigen_symm_workspace *symm_workspace_p;
0154 } gsl_eigen_gensymm_workspace;
0155
0156 gsl_eigen_gensymm_workspace * gsl_eigen_gensymm_alloc (const size_t n);
0157 void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w);
0158 int gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B,
0159                        gsl_vector * eval, gsl_eigen_gensymm_workspace * w);
0160 int gsl_eigen_gensymm_standardize (gsl_matrix * A, const gsl_matrix * B);
0161
0162 typedef struct {
0163   size_t size;            /* size of matrices */
0164   gsl_eigen_symmv_workspace *symmv_workspace_p;
0165 } gsl_eigen_gensymmv_workspace;
0166
0167 gsl_eigen_gensymmv_workspace * gsl_eigen_gensymmv_alloc (const size_t n);
0168 void gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w);
0169 int gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B,
0170                         gsl_vector * eval, gsl_matrix * evec,
0171                         gsl_eigen_gensymmv_workspace * w);
0172
0173 typedef struct {
0174   size_t size;            /* size of matrices */
0175   gsl_eigen_herm_workspace *herm_workspace_p;
0176 } gsl_eigen_genherm_workspace;
0177
0178 gsl_eigen_genherm_workspace * gsl_eigen_genherm_alloc (const size_t n);
0179 void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w);
0180 int gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
0181                        gsl_vector * eval, gsl_eigen_genherm_workspace * w);
0182 int gsl_eigen_genherm_standardize (gsl_matrix_complex * A,
0183                                    const gsl_matrix_complex * B);
0184
0185 typedef struct {
0186   size_t size;            /* size of matrices */
0187   gsl_eigen_hermv_workspace *hermv_workspace_p;
0188 } gsl_eigen_genhermv_workspace;
0189
0190 gsl_eigen_genhermv_workspace * gsl_eigen_genhermv_alloc (const size_t n);
0191 void gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * w);
0192 int gsl_eigen_genhermv (gsl_matrix_complex * A, gsl_matrix_complex * B,
0193                         gsl_vector * eval, gsl_matrix_complex * evec,
0194                         gsl_eigen_genhermv_workspace * w);
0195
0196 typedef struct {
0197   size_t size;            /* size of matrices */
0198   gsl_vector *work;       /* scratch workspace */
0199
0200   size_t n_evals;         /* number of eigenvalues found */
0201   size_t max_iterations;  /* maximum QZ iterations allowed */
0202   size_t n_iter;          /* number of iterations since last eigenvalue found */
0203   double eshift;          /* exceptional shift counter */
0204
0205   int needtop;            /* need to compute top index? */
0206
0207   double atol;            /* tolerance for splitting A matrix */
0208   double btol;            /* tolerance for splitting B matrix */
0209
0210   double ascale;          /* scaling factor for shifts */
0211   double bscale;          /* scaling factor for shifts */
0212
0213   gsl_matrix *H;          /* pointer to hessenberg matrix */
0214   gsl_matrix *R;          /* pointer to upper triangular matrix */
0215
0216   int compute_s;          /* compute generalized Schur form S */
0217   int compute_t;          /* compute generalized Schur form T */
0218
0219   gsl_matrix *Q;          /* pointer to left Schur vectors */
0220   gsl_matrix *Z;          /* pointer to right Schur vectors */
0221 } gsl_eigen_gen_workspace;
0222
0223 gsl_eigen_gen_workspace * gsl_eigen_gen_alloc (const size_t n);
0224 void gsl_eigen_gen_free (gsl_eigen_gen_workspace * w);
0225 void gsl_eigen_gen_params (const int compute_s, const int compute_t,
0226                            const int balance, gsl_eigen_gen_workspace * w);
0227 int gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B,
0228                    gsl_vector_complex * alpha, gsl_vector * beta,
0229                    gsl_eigen_gen_workspace * w);
0230 int gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
0231                       gsl_vector_complex * alpha, gsl_vector * beta,
0232                       gsl_matrix * Q, gsl_matrix * Z,
0233                       gsl_eigen_gen_workspace * w);
0234
0235 typedef struct {
0236   size_t size;            /* size of matrices */
0237
0238   gsl_vector *work1;      /* 1-norm of columns of A */
0239   gsl_vector *work2;      /* 1-norm of columns of B */
0240   gsl_vector *work3;      /* real part of eigenvector */
0241   gsl_vector *work4;      /* imag part of eigenvector */
0242   gsl_vector *work5;      /* real part of back-transformed eigenvector */
0243   gsl_vector *work6;      /* imag part of back-transformed eigenvector */
0244
0245   gsl_matrix *Q;          /* pointer to left Schur vectors */
0246   gsl_matrix *Z;          /* pointer to right Schur vectors */
0247
0248   gsl_eigen_gen_workspace *gen_workspace_p;
0249 } gsl_eigen_genv_workspace;
0250
0251 gsl_eigen_genv_workspace * gsl_eigen_genv_alloc (const size_t n);
0252 void gsl_eigen_genv_free (gsl_eigen_genv_workspace * w);
0253 int gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B,
0254                     gsl_vector_complex * alpha, gsl_vector * beta,
0255                     gsl_matrix_complex * evec,
0256                     gsl_eigen_genv_workspace * w);
0257 int gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
0258                        gsl_vector_complex * alpha, gsl_vector * beta,
0259                        gsl_matrix_complex * evec,
0260                        gsl_matrix * Q, gsl_matrix * Z,
0261                        gsl_eigen_genv_workspace * w);
0262
0263
0264
0265 typedef enum {
0266   GSL_EIGEN_SORT_VAL_ASC,
0267   GSL_EIGEN_SORT_VAL_DESC,
0268   GSL_EIGEN_SORT_ABS_ASC,
0269   GSL_EIGEN_SORT_ABS_DESC
0270 }
0271 gsl_eigen_sort_t;
0272
0273 /* Sort eigensystem results based on eigenvalues.
0274  * Sorts in order of increasing value or increasing
0275  * absolute value.
0276  *
0277  * exceptions: GSL_EBADLEN
0278  */
0279
0280 int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
0281                          gsl_eigen_sort_t sort_type);
0282
0283 int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
0284                          gsl_eigen_sort_t sort_type);
0285
0286 int gsl_eigen_nonsymmv_sort(gsl_vector_complex * eval,
0287                             gsl_matrix_complex * evec,
0288                             gsl_eigen_sort_t sort_type);
0289
0290 int gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec,
0291                              gsl_eigen_sort_t sort_type);
0292
0293 int gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
0294                              gsl_eigen_sort_t sort_type);
0295
0296 int gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
0297                          gsl_matrix_complex * evec,
0298                          gsl_eigen_sort_t sort_type);
0299
0300 /* Prototypes for the schur module */
0301
0302 int gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B,
0303                           double *wr1, double *wr2, double *wi,
0304                           double *scale1, double *scale2);
0305
0306 int gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
0307                              double d1, double d2, const gsl_vector *b,
0308                              gsl_vector *x, double *s, double *xnorm,
0309                              double smin);
0310
0311 int gsl_schur_solve_equation_z(double ca, const gsl_matrix *A,
0312                                gsl_complex *z, double d1, double d2,
0313                                const gsl_vector_complex *b,
0314                                gsl_vector_complex *x, double *s,
0315                                double *xnorm, double smin);
0316
0317
0318 /* The following functions are obsolete: */
0319
0320 /* Eigensolve by Jacobi Method
0321  *
0322  * The data in the matrix input is destroyed.
0323  *
0324  * exceptions:
0325  */
0326 int
0327 gsl_eigen_jacobi(gsl_matrix * matrix,
0328                       gsl_vector * eval,
0329                       gsl_matrix * evec,
0330                       unsigned int max_rot,
0331                       unsigned int * nrot);
0332
0333
0334 /* Invert by Jacobi Method
0335  *
0336  * exceptions:
0337  */
0338 int
0339 gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
0340                              gsl_matrix * ainv,
0341                              unsigned int max_rot);
0342
0343
0344
0345 __END_DECLS
0346
0347 #endif /* __GSL_EIGEN_H__ */