From b4bd1eedacc741d7dbecb9e9be7d37ec720845df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 21 Jan 2020 18:20:22 +0200 Subject: [PATCH] Beyn wrappers for finite system, doxygen Former-commit-id: 236e802bc92c5fb0984f83b00c77555d5c29430d --- qpms/scatsystem.c | 31 +++++++++++++++++++++++++++++++ qpms/scatsystem.h | 43 ++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 71 insertions(+), 3 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 2d645ca..7504d44 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -22,6 +22,7 @@ #include #include "kahansum.h" #include "tolerances.h" +#include "beyn.h" #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS #include "qpmsblas.h" @@ -1894,3 +1895,33 @@ complex double *qpms_scatsys_scatter_solve( return f; } +/// Wrapper for Beyn algorithm (non-periodic system) +static int qpms_scatsys_finite_eval_Beyn_ImTS(complex double *target, + size_t m, complex double omega, void *params) { + const qpms_scatsys_t *ss = params; + qpms_scatsys_at_omega_t *ssw = qpms_scatsys_at_omega(ss, omega); + QPMS_ENSURE(ssw != NULL, "qpms_scatsys_at_omega() returned NULL"); + QPMS_ENSURE(NULL != qpms_scatsysw_build_modeproblem_matrix_full(target, ssw), + "qpms_scatsysw_build_modeproblem_matrix_full() returned NULL"); + qpms_scatsys_at_omega_free(ssw); + return QPMS_SUCCESS; +} + +beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, + complex double omega_centre, double omega_rr, double omega_ri, + size_t contour_npoints, + double rank_tol, size_t rank_sel_min, double res_tol) { + + beyn_contour_t *contour = beyn_contour_ellipse(omega_centre, + omega_rr, omega_ri, contour_npoints); + + beyn_result_t *result = beyn_solve(ss->fecv_size, + ss->fecv_size /* possibly make smaller? */, + qpms_scatsys_finite_eval_Beyn_ImTS, NULL, (void *) ss, + contour, rank_tol, rank_sel_min, res_tol); + QPMS_ENSURE(result != NULL, "beyn_solve() returned NULL"); + + free(contour); + return result; +} + diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 49bc9b8..2aa1a99 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -343,31 +343,43 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( ); /// Creates the full \f$ (I - TS) \f$ matrix of the scattering system. +/** + * \returns \a target on success, NULL on error. + */ complex double *qpms_scatsysw_build_modeproblem_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_at_omega_t *ssw ); /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form. +/** + * \returns \a target on success, NULL on error. + */ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_at_omega_t *ssw, - qpms_iri_t iri + qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym ); /// Alternative implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed(). +/** + * \returns \a target on success, NULL on error. + */ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_at_omega_t *ssw, - qpms_iri_t iri + qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym ); /// Alternative (serial reference) implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed(). +/** + * \returns \a target on success, NULL on error. + */ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_at_omega_t *ssw, - qpms_iri_t iri + qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym ); /// LU factorisation (LAPACKE_zgetrf) result holder. @@ -504,6 +516,31 @@ complex double *qpms_scatsysw_apply_Tmatrices_full( const qpms_scatsys_at_omega_t *ssw ); +struct beyn_result_t; // See beyn.h for full definition + +/// Searches for scattering system's eigenmodes using Beyn's algorithm. +/** + * Currently, elliptical contour is used. + * + * TODO In the future, this will probably support irrep decomposition as well, + * but it does not make much sense for periodic / small systems, as in their + * case the bottleneck is the T-matrix and translation matrix evaluation + * rather than the linear algebra. + */ + +struct beyn_result_t *qpms_scatsys_find_eigenmodes( + const qpms_scatsys_t *ss, + double eta, ///< Ewald sum parameter + const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system. + complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for. + double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for. + double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for. + size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number) + double rank_tol, ///< (default: `1e-4`) TODO DOC. + size_t rank_sel_min, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol. + double res_tol ///< (default: `0.0`) TODO DOC. + ); + #if 0 /// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.