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 <stddef.h>
   8 #include <complex.h>
   9
  10 /// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found.
  11 /** Pure C array version */
  12 typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params);
  13
  14 /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
  15 /** Pure C array version */
  16 typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l,
  17                 const complex double *Vhat, complex double z, void *params);
  18
  19 /// Complex plane integration contour structure.
  20 typedef struct beyn_contour_t {
  21         size_t n; ///< Number of discretisation points.
  22         /// "Centre" of the contour.
  23         /**
  24          * This point is used in the rescaling of the \f$ A_1 \f$ matrix as in
  25          * Beyn's Remark 3.2 (b) in order to improve the numerical stability.
  26          * It does not have to be a centre in some strictly defined sense,
  27          * but it should be "somewhere around" where the contour is.
  28          */
  29         complex double centre;
  30         /// Function testing that a point \a z lies inside the contour (optional).
  31         _Bool (*inside_test)(struct beyn_contour_t *, complex double z);
  32         complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
  33 } beyn_contour_t;
  34
  35 /// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes.
  36 /** Free using free(). */
  37 beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints);
  38
  39
  40 typedef enum {
  41         BEYN_CONTOUR_HALFELLIPSE_RE_PLUS = 3,
  42         BEYN_CONTOUR_HALFELLIPSE_RE_MINUS = 1,
  43         BEYN_CONTOUR_HALFELLIPSE_IM_PLUS = 0,
  44         BEYN_CONTOUR_HALFELLIPSE_IM_MINUS = 2,
  45 } beyn_contour_halfellipse_orientation;
  46
  47
  48 /// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes.
  49 /** Free using free(). */
  50 beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints,
  51                 beyn_contour_halfellipse_orientation or);
  52
  53 /// Similar to halfellipse but with rounded corners.
  54 beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im,
  55         double rounding, ///< Must be in interval [0, 0.5)
  56         size_t n, beyn_contour_halfellipse_orientation or);
  57
  58
  59 /// Beyn algorithm result structure (pure C array version).
  60 typedef struct beyn_result_t {
  61         size_t neig; ///< Number of eigenvalues found.
  62         size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
  63         complex double *eigval;
  64         complex double *eigval_err;
  65         double *residuals;
  66         complex double *eigvec; // Rows are the eigenvectors
  67         double *ranktest_SV;
  68
  69 } beyn_result_t;
  70
  71 void beyn_result_free(beyn_result_t *result);
  72
  73 /// Solve a non-linear eigenproblem using Beyn's algorithm
  74 beyn_result_t *beyn_solve(
  75                 size_t m, ///< Dimension of the matrix \a M.
  76                 size_t l, ///< Number of columns of the random matrix  \f$ \hat V \f$ (larger than the expected number of solutions).
  77                 beyn_function_M_t M, ///< Function providing the matrix \f$ M(z) \f$.
  78                 beyn_function_M_inv_Vhat_t M_inv_Vhat, ///< Fuction providing the matrix \f$ M^{-1}(z) \hat V \f$ (optional).
  79                 void *params, ///< Parameter pointer passed to M() and M_inv_Vhat().
  80                 const beyn_contour_t *contour, ///< Integration contour.
  81                 double rank_tol, ///< (default: `1e-4`) TODO DOC.
  82                 size_t rank_min_sel, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
  83                 double res_tol ///< (default: `0.0`) TODO DOC.
  84               );
  85
  86 #endif // BEYN_H