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
38typedef struct {
39 size_t size;
40 double * d;
41 double * sd;
42} gsl_eigen_symm_workspace;
43
44gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
45void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
46int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);
47
48typedef struct {
49 size_t size;
50 double * d;
51 double * sd;
52 double * gc;
53 double * gs;
54} gsl_eigen_symmv_workspace;
55
56gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
57void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
58int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);
59
60typedef struct {
61 size_t size;
62 double * d;
63 double * sd;
64 double * tau;
65} gsl_eigen_herm_workspace;
66
67gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
68void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
69int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
70 gsl_eigen_herm_workspace * w);
71
72typedef 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
81gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
82void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
83int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
84 gsl_matrix_complex * evec,
85 gsl_eigen_hermv_workspace * w);
86
87typedef 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
99gsl_eigen_francis_workspace * gsl_eigen_francis_alloc (void);
100void gsl_eigen_francis_free (gsl_eigen_francis_workspace * w);
101void gsl_eigen_francis_T (const int compute_t,
102 gsl_eigen_francis_workspace * w);
103int gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
104 gsl_eigen_francis_workspace * w);
105int gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
106 gsl_matrix * Z,
107 gsl_eigen_francis_workspace * w);
108
109typedef 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
120gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n);
121void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w);
122void gsl_eigen_nonsymm_params (const int compute_t, const int balance,
123 gsl_eigen_nonsymm_workspace *w);
124int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
125 gsl_eigen_nonsymm_workspace * w);
126int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
127 gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w);
128
129typedef 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
140gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n);
141void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w);
142void gsl_eigen_nonsymmv_params (const int balance,
143 gsl_eigen_nonsymmv_workspace *w);
144int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
145 gsl_matrix_complex * evec,
146 gsl_eigen_nonsymmv_workspace * w);
147int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
148 gsl_matrix_complex * evec, gsl_matrix * Z,
149 gsl_eigen_nonsymmv_workspace * w);
150
151typedef struct {
152 size_t size; /* size of matrices */
153 gsl_eigen_symm_workspace *symm_workspace_p;
154} gsl_eigen_gensymm_workspace;
155
156gsl_eigen_gensymm_workspace * gsl_eigen_gensymm_alloc (const size_t n);
157void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w);
158int gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B,
159 gsl_vector * eval, gsl_eigen_gensymm_workspace * w);
160int gsl_eigen_gensymm_standardize (gsl_matrix * A, const gsl_matrix * B);
161
162typedef struct {
163 size_t size; /* size of matrices */
164 gsl_eigen_symmv_workspace *symmv_workspace_p;
165} gsl_eigen_gensymmv_workspace;
166
167gsl_eigen_gensymmv_workspace * gsl_eigen_gensymmv_alloc (const size_t n);
168void gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w);
169int gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B,
170 gsl_vector * eval, gsl_matrix * evec,
171 gsl_eigen_gensymmv_workspace * w);
172
173typedef struct {
174 size_t size; /* size of matrices */
175 gsl_eigen_herm_workspace *herm_workspace_p;
176} gsl_eigen_genherm_workspace;
177
178gsl_eigen_genherm_workspace * gsl_eigen_genherm_alloc (const size_t n);
179void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w);
180int gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
181 gsl_vector * eval, gsl_eigen_genherm_workspace * w);
182int gsl_eigen_genherm_standardize (gsl_matrix_complex * A,
183 const gsl_matrix_complex * B);
184
185typedef struct {
186 size_t size; /* size of matrices */
187 gsl_eigen_hermv_workspace *hermv_workspace_p;
188} gsl_eigen_genhermv_workspace;
189
190gsl_eigen_genhermv_workspace * gsl_eigen_genhermv_alloc (const size_t n);
191void gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * w);
192int gsl_eigen_genhermv (gsl_matrix_complex * A, gsl_matrix_complex * B,
193 gsl_vector * eval, gsl_matrix_complex * evec,
194 gsl_eigen_genhermv_workspace * w);
195
196typedef struct {
197 size_t size; /* size of matrices */
198 gsl_vector *work; /* scratch workspace */
199
200 size_t n_evals; /* number of eigenvalues found */
201 size_t max_iterations; /* maximum QZ iterations allowed */
202 size_t n_iter; /* number of iterations since last eigenvalue found */
203 double eshift; /* exceptional shift counter */
204
205 int needtop; /* need to compute top index? */
206
207 double atol; /* tolerance for splitting A matrix */
208 double btol; /* tolerance for splitting B matrix */
209
210 double ascale; /* scaling factor for shifts */
211 double bscale; /* scaling factor for shifts */
212
213 gsl_matrix *H; /* pointer to hessenberg matrix */
214 gsl_matrix *R; /* pointer to upper triangular matrix */
215
216 int compute_s; /* compute generalized Schur form S */
217 int compute_t; /* compute generalized Schur form T */
218
219 gsl_matrix *Q; /* pointer to left Schur vectors */
220 gsl_matrix *Z; /* pointer to right Schur vectors */
221} gsl_eigen_gen_workspace;
222
223gsl_eigen_gen_workspace * gsl_eigen_gen_alloc (const size_t n);
224void gsl_eigen_gen_free (gsl_eigen_gen_workspace * w);
225void gsl_eigen_gen_params (const int compute_s, const int compute_t,
226 const int balance, gsl_eigen_gen_workspace * w);
227int gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B,
228 gsl_vector_complex * alpha, gsl_vector * beta,
229 gsl_eigen_gen_workspace * w);
230int gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
231 gsl_vector_complex * alpha, gsl_vector * beta,
232 gsl_matrix * Q, gsl_matrix * Z,
233 gsl_eigen_gen_workspace * w);
234
235typedef struct {
236 size_t size; /* size of matrices */
237
238 gsl_vector *work1; /* 1-norm of columns of A */
239 gsl_vector *work2; /* 1-norm of columns of B */
240 gsl_vector *work3; /* real part of eigenvector */
241 gsl_vector *work4; /* imag part of eigenvector */
242 gsl_vector *work5; /* real part of back-transformed eigenvector */
243 gsl_vector *work6; /* imag part of back-transformed eigenvector */
244
245 gsl_matrix *Q; /* pointer to left Schur vectors */
246 gsl_matrix *Z; /* pointer to right Schur vectors */
247
248 gsl_eigen_gen_workspace *gen_workspace_p;
249} gsl_eigen_genv_workspace;
250
251gsl_eigen_genv_workspace * gsl_eigen_genv_alloc (const size_t n);
252void gsl_eigen_genv_free (gsl_eigen_genv_workspace * w);
253int gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B,
254 gsl_vector_complex * alpha, gsl_vector * beta,
255 gsl_matrix_complex * evec,
256 gsl_eigen_genv_workspace * w);
257int gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
258 gsl_vector_complex * alpha, gsl_vector * beta,
259 gsl_matrix_complex * evec,
260 gsl_matrix * Q, gsl_matrix * Z,
261 gsl_eigen_genv_workspace * w);
262
263
264
265typedef enum {
266 GSL_EIGEN_SORT_VAL_ASC,
267 GSL_EIGEN_SORT_VAL_DESC,
268 GSL_EIGEN_SORT_ABS_ASC,
269 GSL_EIGEN_SORT_ABS_DESC
270}
271gsl_eigen_sort_t;
272
273/* Sort eigensystem results based on eigenvalues.
274 * Sorts in order of increasing value or increasing
275 * absolute value.
276 *
277 * exceptions: GSL_EBADLEN
278 */
279
280int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
281 gsl_eigen_sort_t sort_type);
282
283int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
284 gsl_eigen_sort_t sort_type);
285
286int gsl_eigen_nonsymmv_sort(gsl_vector_complex * eval,
287 gsl_matrix_complex * evec,
288 gsl_eigen_sort_t sort_type);
289
290int gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec,
291 gsl_eigen_sort_t sort_type);
292
293int gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
294 gsl_eigen_sort_t sort_type);
295
296int gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
297 gsl_matrix_complex * evec,
298 gsl_eigen_sort_t sort_type);
299
300/* Prototypes for the schur module */
301
302int gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B,
303 double *wr1, double *wr2, double *wi,
304 double *scale1, double *scale2);
305
306int gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
307 double d1, double d2, const gsl_vector *b,
308 gsl_vector *x, double *s, double *xnorm,
309 double smin);
310
311int gsl_schur_solve_equation_z(double ca, const gsl_matrix *A,
312 gsl_complex *z, double d1, double d2,
313 const gsl_vector_complex *b,
314 gsl_vector_complex *x, double *s,
315 double *xnorm, double smin);
316
317
318/* The following functions are obsolete: */
319
320/* Eigensolve by Jacobi Method
321 *
322 * The data in the matrix input is destroyed.
323 *
324 * exceptions:
325 */
326int
327gsl_eigen_jacobi(gsl_matrix * matrix,
328 gsl_vector * eval,
329 gsl_matrix * evec,
330 unsigned int max_rot,
331 unsigned int * nrot);
332
333
334/* Invert by Jacobi Method
335 *
336 * exceptions:
337 */
338int
339gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
340 gsl_matrix * ainv,
341 unsigned int max_rot);
342
343
344
345__END_DECLS
346
347#endif /* __GSL_EIGEN_H__ */
348