qpms/beyn.h

   1 /** \file beyn.h
   2  * \brief Beyn's algorithm for nonlinear eigenvalue problems.
   3  */
   4 #ifndef BEYN_H
   5 #define BEYN_H
   6
   7 #include <complex.h>
   8 #include <gsl/gsl_matrix.h>
   9 #include <gsl/gsl_complex_math.h>
  10
  11 /// User-supplied function that provides the operator M(z) whose "roots" are to be found.
  12 /** GSL matrix version */
  13 typedef int (*beyn_function_M_gsl_t)(gsl_matrix_complex *target_M, complex double z, void *params);
  14
  15 /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
  16 /** GSL matrix version */
  17 typedef int (*beyn_function_M_inv_Vhat_gsl_t)(gsl_matrix_complex *target_M_inv_Vhat,
  18                 const gsl_matrix_complex *Vhat, complex double z, void *params);
  19
  20 /// User-supplied function that provides the operator M(z) whose "roots" are to be found.
  21 /** Pure C array version */
  22 typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params);
  23
  24 /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
  25 /** Pure C array version */
  26 typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l,
  27                 const complex double *Vhat, complex double z, void *params);
  28
  29 /// Complex plane integration contour structure.
  30 typedef struct beyn_contour_t {
  31         size_t n; ///< Number of discretisation points.
  32         /// "Centre" of the contour.
  33         /**
  34          * This point is used in the rescaling of the \f$ A_1 \f$ matrix as in
  35          * Beyn's Remark 3.2 (b) in order to improve the numerical stability.
  36          * It does not have to be a centre in some strictly defined sense,
  37          * but it should be "somewhere around" where the contour is.
  38          */
  39         complex double centre;
  40         /// Function testing that a point \a z lies inside the contour (optional).
  41         _Bool (*inside_test)(struct beyn_contour_t *, complex double z);
  42         complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
  43 } beyn_contour_t;
  44
  45 /// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes.
  46 /** Free using free(). */
  47 beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints);
  48
  49
  50 typedef enum {
  51         BEYN_CONTOUR_HALFELLIPSE_RE_PLUS = 3,
  52         BEYN_CONTOUR_HALFELLIPSE_RE_MINUS = 1,
  53         BEYN_CONTOUR_HALFELLIPSE_IM_PLUS = 0,
  54         BEYN_CONTOUR_HALFELLIPSE_IM_MINUS = 2,
  55 } beyn_contour_halfellipse_orientation;
  56
  57
  58 /// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes.
  59 /** Free using free(). */
  60 beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints,
  61                 beyn_contour_halfellipse_orientation or);
  62
  63 /// Similar to halfellipse but with rounded corners.
  64 beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im,
  65         double rounding, ///< Must be in interval [0, 0.5)
  66         size_t n, beyn_contour_halfellipse_orientation or);
  67
  68
  69 /// Beyn algorithm result structure (GSL matrix/vector version).
  70 typedef struct beyn_result_gsl_t {
  71         size_t neig; ///< Number of eigenvalues found (a bit redundant?).
  72         gsl_vector_complex *eigval;
  73         gsl_vector_complex *eigval_err;
  74         gsl_vector *residuals;
  75         gsl_matrix_complex *eigvec; // Rows are the eigenvectors
  76         gsl_vector *ranktest_SV;
  77 } beyn_result_gsl_t;
  78
  79 void beyn_result_gsl_free(beyn_result_gsl_t *result);
  80
  81 /// Beyn algorithm result structure (pure C array version).
  82 typedef struct beyn_result_t {
  83         size_t neig; ///< Number of eigenvalues found.
  84         size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
  85         complex double *eigval;
  86         complex double *eigval_err;
  87         double *residuals;
  88         complex double *eigvec; // Rows are the eigenvectors
  89         double *ranktest_SV;
  90
  91         /// Private, we wrap it around beyn_result_gsl_t for now.
  92         beyn_result_gsl_t *gsl;
  93
  94 } beyn_result_t;
  95
  96 /// Converts beyn_result_gsl_t from beyn_result_t.
  97 /** After calling this, use beyn_result_free() on the returned pointer;
  98  *  do NOT run beyn_result_gsl_free() anymore.
  99  */
 100 beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g);
 101
 102 void beyn_result_free(beyn_result_t *result);
 103
 104 beyn_result_gsl_t *beyn_solve_gsl(
 105                 size_t m, ///< Dimension of the matrix \a M.
 106                 size_t l, ///< Number of columns of the random matrix  \f$ \hat V \f$ (larger than the expected number of solutions).
 107                 beyn_function_M_gsl_t M, ///< Function providing the matrix \f$ M(z) \f$.
 108                 beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, ///< Fuction providing the matrix \f$ M^{-1}(z) \hat V \f$ (optional).
 109                 void *params, ///< Parameter pointer passed to M() and M_inv_Vhat().
 110                 const beyn_contour_t *contour, ///< Integration contour.
 111                 double rank_tol, ///< (default: `1e-4`) TODO DOC.
 112                 size_t rank_min_sel, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
 113                 double res_tol ///< (default: `0.0`) TODO DOC.
 114               );
 115
 116 beyn_result_t *beyn_solve(
 117                 size_t m, ///< Dimension of the matrix \a M.
 118                 size_t l, ///< Number of columns of the random matrix  \f$ \hat V \f$ (larger than the expected number of solutions).
 119                 beyn_function_M_t M, ///< Function providing the matrix \f$ M(z) \f$.
 120                 beyn_function_M_inv_Vhat_t M_inv_Vhat, ///< Fuction providing the matrix \f$ M^{-1}(z) \hat V \f$ (optional).
 121                 void *params, ///< Parameter pointer passed to M() and M_inv_Vhat().
 122                 const beyn_contour_t *contour, ///< Integration contour.
 123                 double rank_tol, ///< (default: `1e-4`) TODO DOC.
 124                 size_t rank_min_sel, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
 125                 double res_tol ///< (default: `0.0`) TODO DOC.
 126               );
 127
 128 static inline complex double gsl_complex_tostd(gsl_complex z) { return GSL_REAL(z) + I*GSL_IMAG(z); }
 129 static inline gsl_complex gsl_complex_fromstd(complex double z) { return gsl_complex_rect(creal(z), cimag(z)); }
 130
 131 #endif // BEYN_H