diff --git a/qpms/groups.h b/qpms/groups.h index 6b05332..e78e462 100644 --- a/qpms/groups.h +++ b/qpms/groups.h @@ -81,6 +81,12 @@ static inline qpms_gmi_t qpms_finite_group_inv(const qpms_finite_group_t *G, return G->invi[a]; } +static inline _Bool qpms_iri_is_valid(const qpms_finite_group_t *G, qpms_iri_t iri) { + return (iri > G->nirreps || iri < 0) ? 0 : 1; +} + + + /// NOT IMPLEMENTED Get irrep index by name. qpms_iri_t qpms_finite_group_find_irrep_by_name(qpms_finite_group_t *G, char *name); diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index fed7007..af16a8d 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -280,7 +280,7 @@ cdef class Particle: if isinstance(tgen.holder, CTMatrix): spec = (tgen.holder).spec else: - raise ValueError("bspec argument must be specified separately for str(type(t))") + raise ValueError("bspec argument must be specified separately for %s" % str(type(t))) self.f = TMatrixFunction(tgen, spec) self.p.tmg = self.f.rawpointer() # TODO non-trivial transformations later; if modified, do not forget to update ScatteringSystem constructor @@ -340,6 +340,16 @@ cdef class ScatteringSystem: #cdef list Tmatrices # Here we keep the references to occuring T-matrices cdef EpsMuGenerator medium_holder # Here we keep the reference to medium generator cdef qpms_scatsys_t *s + cdef FinitePointGroup sym + + cdef qpms_iri_t iri_py2c(self, iri, allow_None = True): + if iri is None and allow_None: + return QPMS_NO_IRREP + cdef qpms_iri_t nir = self.nirreps + cdef qpms_iri_t ciri = iri + if ciri < 0 or ciri > nir: + raise ValueError("Invalid irrep index %s (of %d irreps)", str(iri), self.nirreps) + return ciri def check_s(self): # cdef instead? if self.s == NULL: @@ -416,6 +426,7 @@ cdef class ScatteringSystem: self.medium_holder = mediumgen self.s = ss self.tmgobjs = tmgobjs + self.sym = sym pyssw = _ScatteringSystemAtOmega() pyssw.ssw = ssw pyssw.ss_pyref = self @@ -589,7 +600,7 @@ cdef class ScatteringSystem: self.s, qpms_incfield_planewave, &p, 0) return target_np - def find_modes(self, cdouble omega_centre, double omega_rr, double omega_ri, + def find_modes(self, cdouble omega_centre, double omega_rr, double omega_ri, iri = None, size_t contour_points = 20, double rank_tol = 1e-4, size_t rank_min_sel=1, double res_tol = 0): """ @@ -597,6 +608,7 @@ cdef class ScatteringSystem: """ cdef beyn_result_t *res = qpms_scatsys_finite_find_eigenmodes(self.s, + self.iri_py2c(iri), omega_centre, omega_rr, omega_ri, contour_points, rank_tol, rank_min_sel, res_tol) if res == NULL: raise RuntimeError @@ -635,6 +647,7 @@ cdef class ScatteringSystem: 'residuals':residuals, 'eigval_err':eigval_err, 'ranktest_SV':ranktest_SV, + 'iri': iri, } return retdict diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index f408f5a..a492820 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -107,6 +107,7 @@ cdef extern from "qpms_types.h": ctypedef int32_t qpms_ss_pi_t ctypedef int qpms_gmi_t ctypedef int qpms_iri_t + qpms_iri_t QPMS_NO_IRREP ctypedef const char * qpms_permutation_t struct qpms_tmatrix_t: qpms_vswf_set_spec_t *spec @@ -598,9 +599,9 @@ cdef extern from "scatsystem.h": cdouble *qpms_scatsys_scatter_solve(cdouble *target_f, const cdouble *a_inc, qpms_ss_LU ludata) const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) - beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, cdouble omega_centre, - double omega_rr, double omega_ri, size_t contour_npoints, double rank_tol, size_t rank_sel_min, - double res_tol) + beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, qpms_iri_t iri, + cdouble omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, + double rank_tol, size_t rank_sel_min, double res_tol) cdef extern from "ewald.h": struct qpms_csf_result: diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 1f708de..550d8bd 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -309,6 +309,8 @@ typedef int qpms_gmi_t; /// Irreducible representation index. See also groups.h. typedef int qpms_iri_t; +/// Constant labeling that no irrep decomposition is done (full system solved instead). +#define QPMS_NO_IRREP ((qpms_iri_t) -1) /// Permutation representation of a group element. /** For now, it's just a string of the form "(0,1)(3,4,5)" diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 7504d44..7308cc7 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1895,29 +1895,50 @@ complex double *qpms_scatsys_scatter_solve( return f; } +struct qpms_scatsys_finite_eval_Beyn_ImTS_param { + const qpms_scatsys_t *ss; + qpms_iri_t iri; +}; + /// 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); + const struct qpms_scatsys_finite_eval_Beyn_ImTS_param *p = params; + qpms_scatsys_at_omega_t *ssw = qpms_scatsys_at_omega(p->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"); + if (p->iri == QPMS_NO_IRREP) { + QPMS_ASSERT(m == p->ss->fecv_size); + QPMS_ENSURE(NULL != qpms_scatsysw_build_modeproblem_matrix_full( + target, ssw), + "qpms_scatsysw_build_modeproblem_matrix_full() returned NULL"); + } else { + QPMS_ASSERT(m == p->ss->saecv_sizes[p->iri]); + QPMS_ENSURE(NULL != qpms_scatsysw_build_modeproblem_matrix_irrep_packed( + target, ssw, p->iri), + "qpms_scatsysw_build_modeproblem_matrix_irrep_packed() returned NULL"); + } qpms_scatsys_at_omega_free(ssw); return QPMS_SUCCESS; } -beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, +beyn_result_t *qpms_scatsys_finite_find_eigenmodes( + const qpms_scatsys_t * const ss, const qpms_iri_t iri, 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) { + size_t n; // matrix dimension + if (qpms_iri_is_valid(ss->sym, iri)) { + n = ss->saecv_sizes[iri]; + } else if (iri == QPMS_NO_IRREP) { + n = ss->fecv_size; + } else QPMS_WTF; 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, + + struct qpms_scatsys_finite_eval_Beyn_ImTS_param p = {ss, iri}; + beyn_result_t *result = beyn_solve(n, n /* possibly make smaller? */, + qpms_scatsys_finite_eval_Beyn_ImTS, NULL, (void *) &p, contour, rank_tol, rank_sel_min, res_tol); QPMS_ENSURE(result != NULL, "beyn_solve() returned NULL"); diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 54b9852..5c74c1f 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -528,7 +528,9 @@ struct beyn_result_t; // See beyn.h for full definition * rather than the linear algebra. */ struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes( - const qpms_scatsys_t *ss, + const qpms_scatsys_t *ss, + /// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system. + qpms_iri_t iri, 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.